Compare commits

..
60 changed files with 1317 additions and 4658 deletions

No files matched your search

@@ -15,7 +15,4 @@
<PackageReference Include="BenchmarkDotNet" Version="0.13.12" />
</ItemGroup>
<ItemGroup>
<Compile Remove=".temp.cs" />
</ItemGroup>
</Project>
@@ -50,7 +50,7 @@ public class GpuPow2ModBenchmarks
return Pow2ModSingleBit(Input.Exponent, Input.Modulus);
}
private static GpuUInt128 Pow2ModSingleBit(ulong exponent, GpuUInt128 modulus)
private static GpuUInt128 Pow2ModSingleBit(ulong exponent, in GpuUInt128 modulus)
{
if (modulus.IsZero || modulus == GpuUInt128.One)
{
@@ -1,75 +0,0 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks;
[ShortRunJob(RuntimeMoniker.Net80)]
[MemoryDiagnoser]
public class MutableUInt128AddBenchmarks
{
private static readonly BenchmarkInput[] Inputs =
[
new(Create(0UL, 0x0000_0000_0000_0007UL), Create(0, 0x3UL), "TinyOperands"),
new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), Create(0xFFFF_FFFF, 0x0000_0001UL), "MixedOperands"),
new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), Create(0xFFFF_FFFF, 0xFFFF_FFFBUL), "NearMaxOperands"),
];
[ParamsSource(nameof(GetInputs))]
public BenchmarkInput Input { get; set; }
public static IEnumerable<BenchmarkInput> GetInputs() => Inputs;
private static MutableUInt128 Create(ulong high, ulong low) => new(high, low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void StandardInternal(ref MutableUInt128 source, MutableUInt128 value)
{
ulong temp = source.Low + value.Low;
source.Low = temp;
temp = temp < source.Low ? 1UL : 0UL;
temp += source.High;
temp += value.High;
source.High = temp;
}
/// <summary>
/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.370 ns with mixed operands,
/// 3.361 ns when both limbs and inputs are near their maxima, and 3.368 ns with tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 3.370 ns, NearMaxOperands 3.361 ns, TinyOperands 3.368 ns.
/// </remarks>
[Benchmark(Baseline = true)]
public MutableUInt128 Standard()
{
MutableUInt128 value = Input.InitialValue;
StandardInternal(ref value, Input.Addend);
return value;
}
/// <summary>
/// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 3.175 ns on mixed operands,
/// 3.177 ns on near-max inputs, and 3.178 ns on tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 3.175 ns, NearMaxOperands 3.177 ns, TinyOperands 3.178 ns.
/// Status: Implemented
/// </remarks>
[Benchmark]
public MutableUInt128 Optimized()
{
MutableUInt128 value = Input.InitialValue;
value.Add(Input.Addend);
return value;
}
public readonly record struct BenchmarkInput(MutableUInt128 InitialValue, MutableUInt128 Addend, string Name)
{
public override string ToString()
{
return Name;
}
}
}
@@ -1,100 +0,0 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks;
[ShortRunJob(RuntimeMoniker.Net80)]
[MemoryDiagnoser]
public class MutableUInt128Mul64Benchmarks
{
private static readonly BenchmarkInput[] Inputs =
[
new(Create(0UL, 0x0000_0000_0000_0007UL), 0x3UL, "TinyOperands"),
new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), 0xFFFF_FFFF_0000_0001UL, "MixedOperands"),
new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), 0xFFFF_FFFF_FFFF_FFFBUL, "NearMaxOperands"),
];
[ParamsSource(nameof(GetInputs))]
public BenchmarkInput Input { get; set; }
public static IEnumerable<BenchmarkInput> GetInputs() => Inputs;
private static MutableUInt128 Create(ulong high, ulong low) => new(high, low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong Mul64(ulong left, ulong right, out ulong low)
{
const ulong Mask32 = 0xFFFF_FFFFUL;
ulong leftLow = left & Mask32;
ulong leftHigh = left >> 32;
ulong rightLow = right & Mask32;
ulong rightHigh = right >> 32;
ulong lowLow = leftLow * rightLow;
ulong lowHigh = leftLow * rightHigh;
ulong highLow = leftHigh * rightLow;
ulong highHigh = leftHigh * rightHigh;
ulong carry = (lowLow >> 32) + (lowHigh & Mask32) + (highLow & Mask32);
ulong high = highHigh + (lowHigh >> 32) + (highLow >> 32) + (carry >> 32);
low = ((carry & Mask32) << 32) | (lowLow & Mask32);
return high;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void StandardInternal(ref MutableUInt128 source, ulong value)
{
ulong currentHigh = source.High;
ulong currentLow = source.Low;
ulong carry = Mul64(currentLow, value, out ulong newLow);
Mul64(currentHigh, value, out ulong shiftedHigh);
ulong newHigh = shiftedHigh + carry;
source.Low = newLow;
source.High = newHigh;
}
/// <summary>
/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.113 ns with mixed operands,
/// 3.128 ns when both limbs and inputs are near their maxima, and 3.131 ns with tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 3.113 ns, NearMaxOperands 3.128 ns, TinyOperands 3.131 ns.
/// </remarks>
[Benchmark(Baseline = true)]
public MutableUInt128 Standard()
{
MutableUInt128 value = Input.InitialValue;
StandardInternal(ref value, Input.Addend);
return value;
}
/// <summary>
/// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 2.969 ns on mixed operands,
/// 2.912 ns on near-max inputs, and 2.935 ns on tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 2.969 ns, NearMaxOperands 2.912 ns, TinyOperands 2.935 ns.
/// Status: Implemented
/// </remarks>
[Benchmark]
public MutableUInt128 Optimized()
{
MutableUInt128 value = Input.InitialValue;
value.Mul(Input.Addend);
return value;
}
public readonly record struct BenchmarkInput(MutableUInt128 InitialValue, ulong Addend, string Name)
{
public override string ToString()
{
return Name;
}
}
}
@@ -1,76 +0,0 @@
using System.Collections.Generic;
using System.Numerics;
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks;
[ShortRunJob(RuntimeMoniker.Net80)]
[MemoryDiagnoser]
public class MutableUInt128MultiplyAddBenchmarks
{
private static readonly MultiplyAddInput[] Inputs =
[
new(Create(0UL, 0x0000_0000_0000_0007UL), 0x3UL, 0x5UL, "TinyOperands"),
new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), 0xFFFF_FFFF_0000_0001UL, 0x1357_9BDF_2468_ACEEUL, "MixedOperands"),
new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), 0xFFFF_FFFF_FFFF_FFFBUL, 0xFFFF_FFFF_FFFF_FFF7UL, "NearMaxOperands"),
];
[ParamsSource(nameof(GetInputs))]
public MultiplyAddInput Input { get; set; }
public static IEnumerable<MultiplyAddInput> GetInputs() => Inputs;
private static UInt128 Create(ulong high, ulong low)
{
return ((UInt128)high << 64) | low;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void MultiplyAdd(ref MutableUInt128 source, ulong multiplier, ulong addend)
{
source.Mul(multiplier);
source.Add(addend);
}
/// <summary>
/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.370 ns with mixed operands,
/// 3.361 ns when both limbs and inputs are near their maxima, and 3.368 ns with tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 3.370 ns, NearMaxOperands 3.361 ns, TinyOperands 3.368 ns.
/// </remarks>
[Benchmark(Baseline = true)]
public MutableUInt128 MultiplyAddChained()
{
MutableUInt128 value = new(Input.InitialValue);
MultiplyAdd(ref value, Input.Multiplier, Input.Addend);
return value;
}
/// <summary>
/// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 3.175 ns on mixed operands,
/// 3.177 ns on near-max inputs, and 3.178 ns on tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 3.175 ns, NearMaxOperands 3.177 ns, TinyOperands 3.178 ns.
/// Status: Implemented
/// </remarks>
[Benchmark]
public MutableUInt128 MultiplyAddFullyInlined()
{
MutableUInt128 value = new(Input.InitialValue);
value.MulAdd(Input.Multiplier, Input.Addend);
return value;
}
public readonly record struct MultiplyAddInput(UInt128 InitialValue, ulong Multiplier, ulong Addend, string Name)
{
public override string ToString()
{
return Name;
}
}
}
@@ -1,75 +0,0 @@
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
using PerfectNumbers.Core;
namespace EvenPerfectBitScanner.Benchmarks;
[ShortRunJob(RuntimeMoniker.Net80)]
[MemoryDiagnoser]
public class MutableUInt128SubMutableBenchmarks
{
private static readonly BenchmarkInput[] Inputs =
[
new(Create(0UL, 0x0000_0000_0000_0007UL), Create(0, 0x3UL), "TinyOperands"),
new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), Create(0xFFFF_FFFF, 0x0000_0001UL), "MixedOperands"),
new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), Create(0xFFFF_FFFF, 0xFFFF_FFFBUL), "NearMaxOperands"),
];
[ParamsSource(nameof(GetInputs))]
public BenchmarkInput Input { get; set; }
public static IEnumerable<BenchmarkInput> GetInputs() => Inputs;
private static MutableUInt128 Create(ulong high, ulong low) => new(high, low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void StandardInternal(ref MutableUInt128 source, MutableUInt128 value)
{
ulong low = source.Low;
ulong result = low - value.Low;
ulong borrow = low < value.Low ? 1UL : 0UL;
source.Low = result;
source.High = source.High - value.High - borrow;
}
/// <summary>
/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.370 ns with mixed operands,
/// 3.361 ns when both limbs and inputs are near their maxima, and 3.368 ns with tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 3.370 ns, NearMaxOperands 3.361 ns, TinyOperands 3.368 ns.
/// </remarks>
[Benchmark(Baseline = true)]
public MutableUInt128 Standard()
{
MutableUInt128 value = Input.InitialValue;
StandardInternal(ref value, Input.Addend);
return value;
}
/// <summary>
/// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 3.175 ns on mixed operands,
/// 3.177 ns on near-max inputs, and 3.178 ns on tiny operands.
/// </summary>
/// <remarks>
/// Observed means: MixedOperands 3.175 ns, NearMaxOperands 3.177 ns, TinyOperands 3.178 ns.
/// Status: Implemented
/// </remarks>
[Benchmark]
public MutableUInt128 Optimized()
{
MutableUInt128 value = Input.InitialValue;
value.Sub(Input.Addend);
return value;
}
public readonly record struct BenchmarkInput(MutableUInt128 InitialValue, MutableUInt128 Addend, string Name)
{
public override string ToString()
{
return Name;
}
}
}
@@ -52,31 +52,14 @@ public class Pow2MontgomeryModBenchmarks
/// Baseline right-to-left Montgomery ladder; measured 36.36 μs on the small sample set and 117.32 μs on the large one.
/// </summary>
[Benchmark(Baseline = true)]
public ulong BaselineCpu()
public ulong Baseline()
{
GetData(out var exponents, out var divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false);
}
return checksum;
}
/// <summary>
/// Baseline right-to-left Montgomery ladder; measured ??.?? μs on the small sample set and ??.?? μs on the large one.
/// </summary>
[Benchmark]
public ulong BaselineGpu()
{
GetData(out var exponents, out var divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false);
checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false);
}
return checksum;
@@ -86,31 +69,14 @@ public class Pow2MontgomeryModBenchmarks
/// Uses a known divisor cycle to fold the exponent first, dropping runtimes to 26.04 μs on small inputs (0.72× baseline) and 9.53 μs on large inputs (0.08× baseline).
/// </summary>
[Benchmark]
public ulong BaselineWithKnownCycleCpu()
public ulong BaselineWithKnownCycle()
{
GetData(out var exponents, out var divisors, out var cycles);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]);
}
return checksum;
}
/// <summary>
/// Uses a known divisor cycle to fold the exponent first, dropping runtimes to 26.04 μs on small inputs (0.72× baseline) and 9.53 μs on large inputs (0.08× baseline).
/// </summary>
[Benchmark]
public ulong BaselineWithKnownCycleGpu()
{
GetData(out var exponents, out var divisors, out var cycles);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]);
checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
}
return checksum;
@@ -22,391 +22,320 @@ namespace EvenPerfectBitScanner.Benchmarks;
[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
public class Pow2MontgomeryModCycleComputationBenchmarks
{
private const int SampleCount = 256;
private const ulong OneShiftedLeft63 = 1UL << 63;
private const long OneShiftedLeft60 = 1L << 60;
private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount];
private readonly UInt128[] _wideModuli = new UInt128[SampleCount];
private readonly ulong[] _smallExponents = new ulong[SampleCount];
private readonly ulong[] _largeExponents = new ulong[SampleCount];
private readonly ulong[] _veryLargeExponents = new ulong[SampleCount];
private readonly ulong[] _smallCycles = new ulong[SampleCount];
private readonly ulong[] _largeCycles = new ulong[SampleCount];
private readonly ulong[] _veryLargeCycles = new ulong[SampleCount];
private readonly UInt128[] _wideCycles = new UInt128[SampleCount];
private const int SampleCount = 256;
private const ulong OneShiftedLeft63 = 1UL << 63;
private const long OneShiftedLeft60 = 1L << 60;
private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount];
private readonly UInt128[] _wideModuli = new UInt128[SampleCount];
private readonly ulong[] _smallExponents = new ulong[SampleCount];
private readonly ulong[] _largeExponents = new ulong[SampleCount];
private readonly ulong[] _veryLargeExponents = new ulong[SampleCount];
private readonly ulong[] _smallCycles = new ulong[SampleCount];
private readonly ulong[] _largeCycles = new ulong[SampleCount];
private readonly ulong[] _veryLargeCycles = new ulong[SampleCount];
private readonly UInt128[] _wideCycles = new UInt128[SampleCount];
private readonly Random _random = new(113);
private ulong? _previousPrimeOrder;
private UInt128? _previousWidePrimeOrder;
private readonly Random _random = new(113);
private ulong? _previousPrimeOrder;
private UInt128? _previousWidePrimeOrder;
public enum InputScale
{
Small,
Large,
VeryLarge
}
public enum InputScale
{
Small,
Large,
VeryLarge
}
[ParamsAllValues]
public InputScale Scale { get; set; }
[ParamsAllValues]
public InputScale Scale { get; set; }
[GlobalSetup]
public void Setup()
{
for (int i = 0; i < SampleCount; i++)
{
[GlobalSetup]
public void Setup()
{
for (int i = 0; i < SampleCount; i++)
{
#if DEBUG
Console.WriteLine($"Generating sample {i}");
Console.WriteLine($"Generating sample {i}");
#endif
ulong smallModulus = NextSmallOddModulus();
_smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
_smallExponents[i] = NextSmallExponent();
_smallCycles[i] = CalculateCycleLengthWithHeuristics(smallModulus);
ulong smallModulus = NextSmallOddModulus();
_smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
_smallExponents[i] = NextSmallExponent();
_smallCycles[i] = CalculateCycleLengthWithHeuristics(smallModulus);
(ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
_largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
_largeExponents[i] = NextLargeExponent();
_largeCycles[i] = CalculateCycleLengthWithHeuristics(largeModulus, largeCycle);
(ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
_largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
_largeExponents[i] = NextLargeExponent();
_largeCycles[i] = CalculateCycleLengthWithHeuristics(largeModulus, largeCycle);
ulong veryLargeModulus = NextVeryLargeOddModulus();
_veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus);
_veryLargeExponents[i] = NextVeryLargeExponent();
ulong veryLargeModulus = NextVeryLargeOddModulus();
_veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus);
_veryLargeExponents[i] = NextVeryLargeExponent();
#if DEBUG
Console.WriteLine($"Calculating cycle length {veryLargeModulus}");
Console.WriteLine($"Calculating cycle length {veryLargeModulus}");
#endif
_veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus);
_veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus);
Console.WriteLine($"Generating wide sample {i}");
UInt128 wideModulus = NextWideOddModulus();
_wideModuli[i] = wideModulus;
_wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus);
}
}
Console.WriteLine($"Generating wide sample {i}");
UInt128 wideModulus = NextWideOddModulus();
_wideModuli[i] = wideModulus;
_wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus);
}
}
/// <summary>
/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
/// large set.
/// </summary>
// [Benchmark(Baseline = true)]
public ulong MontgomeryWithoutCycleCpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
/// <summary>
/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
/// large set.
/// </summary>
// [Benchmark(Baseline = true)]
public ulong MontgomeryWithoutCycle()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false);
}
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false);
}
return checksum;
}
return checksum;
}
/// <summary>
/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
/// large set.
/// </summary>
// [Benchmark(Baseline = true)]
public ulong MontgomeryWithoutCycleGpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
/// <summary>
/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
/// on the large one (0.08×).
/// </summary>
[Benchmark]
public ulong MontgomeryWithPrecomputedCycle()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false);
}
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
}
return checksum;
}
return checksum;
}
/// <summary>
/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
/// on the large one (0.08×).
/// </summary>
[Benchmark]
public ulong MontgomeryWithPrecomputedCycleCpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
ulong checksum = 0UL;
/// <summary>
/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
/// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
/// </summary>
// [Benchmark]
public ulong MontgomeryWithGpuCycleComputation()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]);
}
for (int i = 0; i < SampleCount; i++)
{
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
}
return checksum;
}
/// <summary>
/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
/// on the large one (0.08×).
/// </summary>
[Benchmark]
public ulong MontgomeryWithPrecomputedCycleGpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]);
}
return checksum;
}
/// <summary>
/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
/// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
/// </summary>
// [Benchmark]
public ulong MontgomeryWithGpuCycleComputationCpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
}
return checksum;
}
/// <summary>
/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
/// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
/// </summary>
// [Benchmark]
public ulong MontgomeryWithGpuCycleComputationGpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
}
return checksum;
}
[Benchmark]
public ulong MontgomeryWithHeuristicCycleComputationCpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
}
return checksum;
}
[Benchmark]
public ulong MontgomeryWithHeuristicCycleComputationGpu()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
}
return checksum;
}
[Benchmark]
public UInt128 HeuristicWideCycleComputation()
{
UInt128 checksum = UInt128.Zero;
for (int i = 0; i < SampleCount; i++)
{
UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
_wideCycles[i] = cycle;
checksum ^= cycle;
}
return checksum;
}
private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
{
switch (Scale)
{
case InputScale.Small:
exponents = _smallExponents;
divisors = _smallDivisors;
cycles = _smallCycles;
break;
case InputScale.Large:
exponents = _largeExponents;
divisors = _largeDivisors;
cycles = _largeCycles;
break;
case InputScale.VeryLarge:
exponents = _veryLargeExponents;
divisors = _veryLargeDivisors;
cycles = _veryLargeCycles;
break;
default:
throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null);
}
}
private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue);
private ulong NextVeryLargeExponent()
{
Span<byte> buffer = stackalloc byte[8];
_random.NextBytes(buffer);
ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
return value;
}
private ulong NextSmallOddModulus()
{
ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
return value | 1UL;
}
private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
{
int bitLength = _random.Next(48, 64);
ulong modulus = (1UL << bitLength) - 1UL;
return (modulus, (ulong)bitLength);
}
private ulong NextVeryLargeOddModulus()
{
Span<byte> buffer = stackalloc byte[8];
_random.NextBytes(buffer);
ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
return value;
}
private UInt128 NextWideOddModulus()
{
Span<byte> buffer = stackalloc byte[16];
_random.NextBytes(buffer);
ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
Span<byte> highSlice = buffer.Slice(8, 8);
ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
UInt128 value = ((UInt128)high << 64) | low;
return value | UInt128.One;
}
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
{
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL);
}
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
nPrime,
montgomeryOne,
montgomeryTwo,
montgomeryTwoSquared);
}
return checksum;
}
private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
{
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
}
[Benchmark]
public ulong MontgomeryWithHeuristicCycleComputation()
{
GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
ulong checksum = 0UL;
// When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
for (int i = 0; i < SampleCount; i++)
{
ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
}
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(modulus);
// // bool isPrime = PrimeTester.IsPrimeInternal(modulus, CancellationToken.None);
// if (!isPrime)
// {
// return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
// }
return checksum;
}
[Benchmark]
public UInt128 HeuristicWideCycleComputation()
{
UInt128 checksum = UInt128.Zero;
for (int i = 0; i < SampleCount; i++)
{
UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
_wideCycles[i] = cycle;
checksum ^= cycle;
}
return checksum;
}
private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
{
switch (Scale)
{
case InputScale.Small:
exponents = _smallExponents;
divisors = _smallDivisors;
cycles = _smallCycles;
break;
case InputScale.Large:
exponents = _largeExponents;
divisors = _largeDivisors;
cycles = _largeCycles;
break;
case InputScale.VeryLarge:
exponents = _veryLargeExponents;
divisors = _veryLargeDivisors;
cycles = _veryLargeCycles;
break;
default:
throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null);
}
}
private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue);
private ulong NextVeryLargeExponent()
{
Span<byte> buffer = stackalloc byte[8];
_random.NextBytes(buffer);
ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
return value;
}
private ulong NextSmallOddModulus()
{
ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
return value | 1UL;
}
private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
{
int bitLength = _random.Next(48, 64);
ulong modulus = (1UL << bitLength) - 1UL;
return (modulus, (ulong)bitLength);
}
private ulong NextVeryLargeOddModulus()
{
Span<byte> buffer = stackalloc byte[8];
_random.NextBytes(buffer);
ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
return value;
}
private UInt128 NextWideOddModulus()
{
Span<byte> buffer = stackalloc byte[16];
_random.NextBytes(buffer);
ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
Span<byte> highSlice = buffer.Slice(8, 8);
ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
UInt128 value = ((UInt128)high << 64) | low;
return value | UInt128.One;
}
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
{
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL);
}
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
return new MontgomeryDivisorData(
modulus,
nPrime,
montgomeryOne,
montgomeryTwo,
montgomeryTwoSquared);
}
private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
{
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
}
// When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(modulus);
// // bool isPrime = PrimeTester.IsPrimeInternal(modulus, CancellationToken.None);
// if (!isPrime)
// {
// return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
// }
#if DEBUG
Console.WriteLine("Trying heuristic. Prime order calculation");
Console.WriteLine("Trying heuristic. Prime order calculation");
#endif
ulong order = PrimeOrderCalculator.Calculate(
modulus,
_previousPrimeOrder,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
modulus,
_previousPrimeOrder,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
if (order != 0UL)
{
_previousPrimeOrder = order;
return order;
}
if (orderResult.Order != 0UL)
{
_previousPrimeOrder = orderResult.Order;
return orderResult.Order;
}
_previousPrimeOrder = null;
_previousPrimeOrder = null;
Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation");
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
}
Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation");
return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
}
private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
{
if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero)
{
return MersenneDivisorCycles.GetCycle(modulus);
}
private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
{
if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero)
{
return MersenneDivisorCycles.GetCycle(modulus);
}
#if DEBUG
Console.WriteLine("Trying heuristic. Wide prime order calculation");
Console.WriteLine("Trying heuristic. Wide prime order calculation");
#endif
UInt128 order = PrimeOrderCalculator.Calculate(
modulus,
_previousWidePrimeOrder,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderResultWide orderResult = PrimeOrderCalculator.Calculate(
modulus,
_previousWidePrimeOrder,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
if (order != UInt128.Zero)
{
_previousWidePrimeOrder = order;
return order;
}
if (orderResult.Order != UInt128.Zero)
{
_previousWidePrimeOrder = orderResult.Order;
return orderResult.Order;
}
_previousWidePrimeOrder = null;
_previousWidePrimeOrder = null;
Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation");
return MersenneDivisorCycles.GetCycle(modulus);
}
Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation");
return MersenneDivisorCycles.GetCycle(modulus);
}
private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
}
@@ -176,7 +176,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
return (UInt128)result;
}
private static bool ShouldUseSingleBit(GpuUInt128 exponent)
private static bool ShouldUseSingleBit(in GpuUInt128 exponent)
{
if (exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold)
{
@@ -257,7 +257,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
return result;
}
private static bool IsBitSet(GpuUInt128 value, int bitIndex)
private static bool IsBitSet(in GpuUInt128 value, int bitIndex)
{
if (bitIndex >= 64)
{
@@ -267,7 +267,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
return ((value.Low >> bitIndex) & 1UL) != 0UL;
}
private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount)
private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount)
{
if (windowStart != 0)
{
+8 -14
View File
@@ -8,22 +8,16 @@ public static class Program
{
public static void Main(string[] args)
{
if (args.Length > 0)
{
IConfig config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary | ConfigOptions.DisableOptimizationsValidator);
BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config);
return;
}
// var config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary | ConfigOptions.DisableOptimizationsValidator);
// BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config);
BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run();
Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new()
{
Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large
};
// Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new()
// {
// Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large
// };
// GpuPrimeWorkLimiter.SetLimit(1);
// benchmarks.Setup();
GpuPrimeWorkLimiter.SetLimit(1);
benchmarks.Setup();
// Console.WriteLine("Press ENTER to run");
// Console.ReadLine();
@@ -11,7 +11,4 @@
<PackageReference Include="Open.Numeric.Primes" Version="4.0.4" />
</ItemGroup>
<ItemGroup>
<Compile Remove=".temp.cs" />
</ItemGroup>
</Project>
@@ -22,7 +22,4 @@
<ItemGroup>
<ProjectReference Include="../EvenPerfectBitScanner/EvenPerfectBitScanner.csproj" />
</ItemGroup>
<ItemGroup>
<Compile Remove=".temp.cs" />
</ItemGroup>
</Project>
@@ -16,7 +16,4 @@
<PackageReference Include="Open.Numeric" Version="3.0.1" />
<PackageReference Include="Open.Numeric.Primes" Version="4.0.4" />
</ItemGroup>
<ItemGroup>
<Compile Remove=".temp.cs" />
</ItemGroup>
</Project>
+15 -119
View File
@@ -44,7 +44,6 @@ internal static class Program
private static string? _resultsPrefix;
private static readonly Optimized PrimeIterator = new();
private static ulong _byDivisorStartPrime;
private static TimeSpan? _primeTestLimit;
[ThreadStatic]
private static bool _lastCompositeP;
@@ -75,8 +74,6 @@ internal static class Program
bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental)
bool useDivisor = false; // M_p divisibility by specific divisor
bool useByDivisor = false; // Iterative divisor scan across primes UInt128 divisor = UInt128.Zero;
ByDivisorDeltasDevice byDivisorDeltasDevice = ByDivisorDeltasDevice.Cpu;
ByDivisorMontgomeryDevice byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Cpu;
// Device routing
bool useGpuCycles = true;
bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device
@@ -184,51 +181,6 @@ internal static class Program
{
mersenneOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
}
else if (arg.StartsWith("--bydivisor-deltas-device=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> deviceValue = arg.AsSpan(arg.IndexOf('=') + 1);
if (deviceValue.Equals("gpu", StringComparison.OrdinalIgnoreCase))
{
byDivisorDeltasDevice = ByDivisorDeltasDevice.Gpu;
}
else if (deviceValue.Equals("cpu", StringComparison.OrdinalIgnoreCase))
{
byDivisorDeltasDevice = ByDivisorDeltasDevice.Cpu;
}
else
{
Console.WriteLine("Invalid --bydivisor-deltas-device value. Expected cpu or gpu.");
return;
}
}
else if (arg.StartsWith("--bydivisor-montgomery-device=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> deviceValue = arg.AsSpan(arg.IndexOf('=') + 1);
if (deviceValue.Equals("gpu", StringComparison.OrdinalIgnoreCase))
{
byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Gpu;
}
else if (deviceValue.Equals("cpu", StringComparison.OrdinalIgnoreCase))
{
byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Cpu;
}
else
{
Console.WriteLine("Invalid --bydivisor-montgomery-device value. Expected cpu or gpu.");
return;
}
}
else if (arg.StartsWith("--prime-test-limit=", StringComparison.OrdinalIgnoreCase))
{
ReadOnlySpan<char> value = arg.AsSpan(arg.IndexOf('=') + 1);
if (!TryParsePrimeTestLimit(value, out TimeSpan parsedLimitValue))
{
Console.WriteLine("Failed to parse --prime-test-limit value.");
return;
}
_primeTestLimit = parsedLimitValue;
}
else if (arg.Equals("--workaround-mod", StringComparison.OrdinalIgnoreCase))
{
useModuloWorkaround = true;
@@ -425,8 +377,6 @@ internal static class Program
}
}
GpuKernelPool.CpuCycleUsesGpu = orderOnGpu;
if (useByDivisor)
{
_byDivisorStartPrime = startPrimeProvided ? currentP : 0UL;
@@ -567,11 +517,6 @@ internal static class Program
? new MersenneNumberDivisorByDivisorGpuTester()
: new MersenneNumberDivisorByDivisorCpuTester();
_byDivisorTester.BatchSize = scanBatchSize;
if (_byDivisorTester is MersenneNumberDivisorByDivisorCpuTester cpuByDivisorTester)
{
cpuByDivisorTester.DeltasDevice = byDivisorDeltasDevice;
cpuByDivisorTester.MontgomeryDevice = byDivisorMontgomeryDevice;
}
}
// Load RLE blacklist (optional)
@@ -607,18 +552,18 @@ internal static class Program
useDivisor,
useByDivisor,
mersenneOnGpu,
useOrder,
useModuloWorkaround,
_useGcdFilter,
NttGpuMath.GpuTransformBackend,
gpuPrimeThreads,
sliceSize,
scanBatchSize,
_orderWarmupLimitOverride ?? 5_000_000UL,
NttGpuMath.ReductionMode,
mersenneOnGpu ? "gpu" : "cpu",
(GpuContextPool.ForceCpu ? "cpu" : "gpu"),
orderOnGpu ? "gpu" : "cpu");
useOrder,
useModuloWorkaround,
_useGcdFilter,
NttGpuMath.GpuTransformBackend,
gpuPrimeThreads,
sliceSize,
scanBatchSize,
_orderWarmupLimitOverride ?? 5_000_000UL,
NttGpuMath.ReductionMode,
mersenneOnGpu ? "gpu" : "cpu",
(GpuContextPool.ForceCpu ? "cpu" : "gpu"),
orderOnGpu ? "gpu" : "cpu");
if (!string.IsNullOrEmpty(_resultsPrefix))
{
@@ -720,7 +665,7 @@ internal static class Program
if (_useByDivisorMode)
{
RunByDivisorMode(byDivisorCandidates, threadCount, _primeTestLimit);
RunByDivisorMode(byDivisorCandidates, threadCount);
FlushBuffer();
StringBuilderPool.Return(_outputBuilder!);
return;
@@ -861,7 +806,7 @@ internal static class Program
return candidates;
}
private static void RunByDivisorMode(List<ulong> candidates, int threadCount, TimeSpan? primeTestLimit)
private static void RunByDivisorMode(List<ulong> candidates, int threadCount)
{
if (_byDivisorTester is null)
{
@@ -878,8 +823,7 @@ internal static class Program
static () => _lastCompositeP = true,
static () => _lastCompositeP = false,
PrintResult,
threadCount,
primeTestLimit);
threadCount);
}
private static unsafe void LoadResultsFile(string resultsFileName, Action<ulong, bool, bool> lineProcessorAction)
@@ -946,51 +890,6 @@ internal static class Program
}
}
private static bool TryParsePrimeTestLimit(ReadOnlySpan<char> value, out TimeSpan result)
{
result = TimeSpan.Zero;
if (value.IsEmpty)
{
return false;
}
ReadOnlySpan<char> magnitudeSpan = value;
double magnitude;
if (value.Length >= 2 && value.EndsWith("ms", StringComparison.OrdinalIgnoreCase))
{
magnitudeSpan = value[..^2];
if (!double.TryParse(magnitudeSpan, NumberStyles.Float, CultureInfo.InvariantCulture, out magnitude))
{
return false;
}
result = TimeSpan.FromMilliseconds(magnitude);
return true;
}
char suffix = char.ToLowerInvariant(value[^1]);
if (suffix == 's' || suffix == 'm')
{
magnitudeSpan = value[..^1];
if (!double.TryParse(magnitudeSpan, NumberStyles.Float, CultureInfo.InvariantCulture, out magnitude))
{
return false;
}
result = suffix == 's' ? TimeSpan.FromSeconds(magnitude) : TimeSpan.FromMinutes(magnitude);
return true;
}
if (!double.TryParse(value, NumberStyles.Float, CultureInfo.InvariantCulture, out magnitude))
{
return false;
}
result = TimeSpan.FromMilliseconds(magnitude);
return true;
}
private static void PrintHelp()
{
Console.WriteLine("Usage: EvenPerfectBitScanner [options]");
@@ -1024,9 +923,6 @@ internal static class Program
Console.WriteLine(" --divisor-cycles-batch=<value> batch size for cycles generation (default 512)");
Console.WriteLine(" --divisor-cycles-continue continue divisor cycles generation");
Console.WriteLine(" --divisor-cycles-limit=<value> cycle search iterations when --mersenne=divisor");
Console.WriteLine(" --prime-test-limit=<value> Time limit for --mersenne=bydivisor prime checks (e.g. 5s, 500ms)");
Console.WriteLine(" --bydivisor-deltas-device=cpu|gpu Device for --mersenne=bydivisor residue batches (default: cpu)");
Console.WriteLine(" --bydivisor-montgomery-device=cpu|gpu Device for Montgomery data when using --mersenne=bydivisor (default: cpu)");
Console.WriteLine(" --use-order test primality via q order");
Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU");
// mod-automaton removed
@@ -15,15 +15,9 @@ public class GpuContextPoolTests
{
GpuContextPool.DisposeAll();
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
}
finally
{
lease.Dispose();
}
using var lease = GpuContextPool.RentPreferred(preferCpu: true);
lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
GpuContextPool.DisposeAll();
}
@@ -36,27 +30,17 @@ public class GpuContextPoolTests
Context firstContext;
Accelerator firstAccelerator;
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
{
firstContext = lease.Context;
firstAccelerator = lease.Accelerator;
}
finally
{
lease.Dispose();
}
lease = GpuContextPool.RentPreferred(preferCpu: true);
try
using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
{
lease.Context.Should().BeSameAs(firstContext);
lease.Accelerator.Should().BeSameAs(firstAccelerator);
}
finally
{
lease.Dispose();
}
GpuContextPool.DisposeAll();
}
@@ -68,27 +52,17 @@ public class GpuContextPoolTests
GpuContextPool.DisposeAll();
Context firstContext;
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
{
firstContext = lease.Context;
}
finally
{
lease.Dispose();
}
GpuContextPool.DisposeAll();
lease = GpuContextPool.RentPreferred(preferCpu: true);
try
using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
{
lease.Context.Should().NotBeSameAs(firstContext);
}
finally
{
lease.Dispose();
}
GpuContextPool.DisposeAll();
}
@@ -199,17 +199,10 @@ public class GpuKernelCompilationTests
[MemberData(nameof(KernelLoaders))]
public void Kernel_compiles_without_internal_compiler_errors(string name, Action<Accelerator> compile)
{
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
var accelerator = lease.Accelerator;
Action action = () => compile(accelerator);
action.Should().NotThrow($"kernel {name} should compile");
}
finally
{
lease.Dispose();
}
using var lease = GpuContextPool.RentPreferred(preferCpu: true);
var accelerator = lease.Accelerator;
Action action = () => compile(accelerator);
action.Should().NotThrow($"kernel {name} should compile");
}
private static void CompileDivisorCycleCacheKernel(Accelerator accelerator)
@@ -150,36 +150,29 @@ public class NttGpuMathTests
var modulus = new GpuUInt128(0UL, 17UL);
var primitiveRoot = new GpuUInt128(0UL, 3UL);
GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
try
using var gpu = GpuContextPool.Rent();
var accelerator = gpu.Accelerator;
var method = typeof(NttGpuMath).GetMethod("GetSquareCache", BindingFlags.NonPublic | BindingFlags.Static);
method.Should().NotBeNull();
var cache = method!.Invoke(null, new object[] { accelerator, 8, modulus, primitiveRoot })!;
var cacheType = cache.GetType();
var twiddlesMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
var twiddlesInvMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesInvMont")!.GetValue(cache)!;
var r2 = (ulong)cacheType.GetProperty("MontR2Mod64")!.GetValue(cache)!;
var mod = (ulong)cacheType.GetProperty("ModulusLow")!.GetValue(cache)!;
int count = 7; // length 8 -> 7 twiddles
var forward = new GpuUInt128[count];
var inverse = new GpuUInt128[count];
twiddlesMont.CopyToCPU(forward);
twiddlesInvMont.CopyToCPU(inverse);
for (int i = 0; i < count; i++)
{
var accelerator = gpu.Accelerator;
var method = typeof(NttGpuMath).GetMethod("GetSquareCache", BindingFlags.NonPublic | BindingFlags.Static);
method.Should().NotBeNull();
var cache = method!.Invoke(null, new object[] { accelerator, 8, modulus, primitiveRoot })!;
var cacheType = cache.GetType();
var twiddlesMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
var twiddlesInvMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesInvMont")!.GetValue(cache)!;
var r2 = (ulong)cacheType.GetProperty("MontR2Mod64")!.GetValue(cache)!;
var mod = (ulong)cacheType.GetProperty("ModulusLow")!.GetValue(cache)!;
int count = 7; // length 8 -> 7 twiddles
var forward = new GpuUInt128[count];
var inverse = new GpuUInt128[count];
twiddlesMont.CopyToCPU(forward);
twiddlesInvMont.CopyToCPU(inverse);
for (int i = 0; i < count; i++)
{
ulong prod = (ulong)((((BigInteger)forward[i].Low) * inverse[i].Low) % mod);
prod.Should().Be(r2);
}
}
finally
{
gpu.Dispose();
ulong prod = (ulong)((((BigInteger)forward[i].Low) * inverse[i].Low) % mod);
prod.Should().Be(r2);
}
}
@@ -345,23 +338,12 @@ public class NttGpuMathTests
values[1] = new GpuUInt128(0UL, 2UL);
var modulus = new GpuUInt128(0UL, 17UL);
var primitiveRoot = new GpuUInt128(0UL, 3UL);
GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
try
using (var gpu = GpuContextPool.Rent())
{
var buffer = gpu.Accelerator.Allocate1D<GpuUInt128>(values.Length);
try
{
buffer.View.CopyFromCPU(ref values[0], values.Length);
NttGpuMath.SquareDevice(gpu.Accelerator, buffer.View, modulus, primitiveRoot);
}
finally
{
buffer.Dispose();
}
}
finally
{
gpu.Dispose();
buffer.View.CopyFromCPU(ref values[0], values.Length);
NttGpuMath.SquareDevice(gpu.Accelerator, buffer.View, modulus, primitiveRoot);
buffer.Dispose();
}
// Verify that the square cache now holds entries for the accelerator.
@@ -33,23 +33,16 @@ public class NttMontgomeryTests
NttGpuMath.Square(cpu, modulus, primitiveRoot);
// Device execution
GpuContextPool.GpuContextLease context = GpuContextPool.Rent();
try
{
var accelerator = context.Accelerator;
using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
using var context = GpuContextPool.Rent();
var accelerator = context.Accelerator;
using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
for (int i = 0; i < cpu.Length; i++)
{
gpu[i].Should().Be(cpu[i]);
}
}
finally
for (int i = 0; i < cpu.Length; i++)
{
context.Dispose();
gpu[i].Should().Be(cpu[i]);
}
}
@@ -359,9 +359,10 @@ public class PrimeOrderGpuHeuristicsTests
PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
ulong result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
result.Should().Be(3UL);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
result.Order.Should().Be(3UL);
}
finally
{
@@ -11,21 +11,14 @@ public class GpuKernelPoolTests
[Trait("Category", "Fast")]
public void EnsureSmallCyclesOnDevice_returns_same_view_on_subsequent_calls()
{
GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
try
{
var accelerator = lease.Accelerator;
using var lease = GpuContextPool.RentPreferred(preferCpu: true);
var accelerator = lease.Accelerator;
var first = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
var second = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
var first = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
var second = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
second.Equals(first).Should().BeTrue();
second.Length.Should().Be(MersenneDivisorCycles.SmallDivisorsMax + 1);
}
finally
{
lease.Dispose();
}
second.Equals(first).Should().BeTrue();
second.Length.Should().Be(MersenneDivisorCycles.SmallDivisorsMax + 1);
}
[Fact]
@@ -26,8 +26,7 @@ public class MersenneNumberDivisorByDivisorTesterTests
markComposite: () => compositeMarks++,
clearComposite: () => clearedMarks++,
printResult: (prime, searched, detailed, passed) => results.Add((prime, searched, detailed, passed)),
threadCount: 1,
primeTestLimit: null);
threadCount: 1);
tester.ConfiguredMaxPrime.Should().Be(13UL);
compositeMarks.Should().Be(2);
@@ -60,7 +59,7 @@ public class MersenneNumberDivisorByDivisorTesterTests
return DivisorLimit;
}
public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
public bool IsPrime(ulong prime, out bool divisorsExhausted)
{
divisorsExhausted = true;
return prime % 7UL != 0UL && prime % 11UL != 0UL;
@@ -7,88 +7,78 @@ namespace PerfectNumbers.Core.Tests;
public class MersenneNumberDivisorCpuTesterTests
{
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_tester_tracks_divisors_across_primes()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(11UL);
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_tester_tracks_divisors_across_primes()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(11UL);
tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeTrue();
tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeTrue();
tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeTrue();
tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue();
divisorsExhausted.Should().BeTrue();
tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse();
divisorsExhausted.Should().BeTrue();
}
tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse();
divisorsExhausted.Should().BeTrue();
}
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_prime_test_limit_disables_divisor_exhaustion()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(11UL);
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_session_checks_divisors_across_primes()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(13UL);
tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
divisorsExhausted.Should().BeFalse();
}
using var session = tester.CreateDivisorSession();
ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
byte[] hits = new byte[primes.Length];
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_session_checks_divisors_across_primes()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(13UL);
ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
using var session = tester.CreateDivisorSession();
ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
byte[] hits = new byte[primes.Length];
Array.Fill(hits, (byte)0);
ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
}
ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(43UL);
Array.Fill(hits, (byte)0);
ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
}
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(43UL);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
byte[] hits = new byte[exponents.Length];
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(61UL);
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
{
var tester = new MersenneNumberDivisorByDivisorCpuTester();
tester.ConfigureFromMaxPrime(61UL);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
using var session = tester.CreateDivisorSession();
ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
byte[] hits = new byte[exponents.Length];
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
}
}
@@ -86,17 +86,6 @@ public class MersenneNumberDivisorGpuTesterTests
exhausted.Should().BeTrue();
}
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_gpu_prime_test_limit_stops_scan_immediately()
{
var tester = new MersenneNumberDivisorByDivisorGpuTester();
tester.ConfigureFromMaxPrime(43UL);
tester.IsPrime(31UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
divisorsExhausted.Should().BeFalse();
}
[Fact]
[Trait("Category", "Fast")]
public void ByDivisor_tester_requires_configuration()
@@ -22,7 +22,4 @@
<ItemGroup>
<ProjectReference Include="../PerfectNumbers.Core/PerfectNumbers.Core.csproj" />
</ItemGroup>
<ItemGroup>
<Compile Remove=".temp.cs" />
</ItemGroup>
</Project>
@@ -7,9 +7,6 @@ namespace PerfectNumbers.Core.Tests;
public class PrimeOrderCalculatorTests
{
// EvenPerfectBitScanner never exercises the GPU heuristic device or the small-prime guards that the legacy tests cover,
// so those scenarios are commented out because they are unreachable in the production pipeline.
/*
public static IEnumerable<object[]> SmallPrimes()
{
yield return new object[] { 3UL };
@@ -48,14 +45,14 @@ public class PrimeOrderCalculatorTests
{
ulong expected = ComputeOrderByDoubling(prime);
ulong result = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
result.Should().Be(expected);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
result.Order.Should().Be(expected);
}
[Theory]
@@ -65,21 +62,21 @@ public class PrimeOrderCalculatorTests
{
ulong previousOrder = ComputeOrderByDoubling(previousPrime);
ulong heuristic = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
prime,
previousOrder,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
ulong strict = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault);
heuristic.Should().Be(strict);
heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
heuristic.Order.Should().Be(strict.Order);
}
[Theory]
@@ -88,14 +85,14 @@ public class PrimeOrderCalculatorTests
[InlineData(3UL, 2UL)]
public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder)
{
ulong result = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
result.Should().Be(expectedOrder);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
result.Order.Should().Be(expectedOrder);
}
[Fact]
@@ -108,14 +105,14 @@ public class PrimeOrderCalculatorTests
maxPowChecks: 8,
mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
ulong result = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
prime: 13UL,
previousOrder: null,
MontgomeryDivisorData.FromModulus(13UL),
config,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
config);
result.Should().Be(ComputeOrderByDoubling(13UL));
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
result.Order.Should().Be(ComputeOrderByDoubling(13UL));
}
[Fact]
@@ -129,14 +126,14 @@ public class PrimeOrderCalculatorTests
mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
ulong heuristic = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
prime: 239UL,
previousOrder: null,
divisorData,
heuristicConfig,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
heuristicConfig);
heuristic.Should().Be(ComputeOrderByDoubling(239UL));
heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
heuristic.Order.Should().Be(ComputeOrderByDoubling(239UL));
var strictConfig = new PrimeOrderCalculator.PrimeOrderSearchConfig(
smallFactorLimit: 2,
@@ -144,14 +141,14 @@ public class PrimeOrderCalculatorTests
maxPowChecks: 1,
mode: PrimeOrderCalculator.PrimeOrderMode.Strict);
ulong strict = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
prime: 239UL,
previousOrder: null,
divisorData,
strictConfig,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
strictConfig);
strict.Should().Be(heuristic);
strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
strict.Order.Should().Be(heuristic.Order);
}
[Fact]
@@ -160,13 +157,13 @@ public class PrimeOrderCalculatorTests
{
UInt128 prime = UInt128.Parse("18446744073709641691");
UInt128 result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
PrimeOrderCalculator.PrimeOrderResultWide result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
result.Should().Be((UInt128)1229782938247309446UL);
result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
result.Order.Should().Be((UInt128)1229782938247309446UL);
}
[Fact]
@@ -200,5 +197,4 @@ public class PrimeOrderCalculatorTests
return order;
}
*/
}
@@ -1,7 +0,0 @@
namespace PerfectNumbers.Core;
public enum ByDivisorDeltasDevice
{
Cpu,
Gpu
}
@@ -1,7 +0,0 @@
namespace PerfectNumbers.Core;
public enum ByDivisorMontgomeryDevice
{
Cpu = 0,
Gpu = 1,
}
@@ -1,10 +1,8 @@
using System;
using System.Buffers;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Numerics;
using PerfectNumbers.Core;
using PerfectNumbers.Core.Gpu;
namespace PerfectNumbers.Core.Cpu;
@@ -17,49 +15,17 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
private bool _isConfigured;
private int _batchSize = 1_024;
private ulong _configuredMaxPrime;
private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu;
private ByDivisorMontgomeryDevice _montgomeryDevice = ByDivisorMontgomeryDevice.Cpu;
public int BatchSize
{
get => _batchSize;
set
{
int sanitized = Math.Max(1, value);
_batchSize = sanitized;
}
}
public ByDivisorDeltasDevice DeltasDevice
{
get => _deltasDevice;
set
{
lock (_sync)
{
_deltasDevice = value;
}
}
}
public ByDivisorMontgomeryDevice MontgomeryDevice
{
get => _montgomeryDevice;
set
{
lock (_sync)
{
_montgomeryDevice = value;
}
}
set => _batchSize = Math.Max(1, value);
}
public void ConfigureFromMaxPrime(ulong maxPrime)
{
lock (_sync)
{
_configuredMaxPrime = maxPrime;
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
_lastStatusDivisor = 0UL;
_isConfigured = true;
@@ -95,7 +61,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
}
public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
public bool IsPrime(ulong prime, out bool divisorsExhausted)
{
ulong allowedMax;
lock (_sync)
@@ -114,17 +80,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
return true;
}
ulong processedCount;
ulong lastProcessed;
bool processedAll;
ulong processedCount = 0UL;
ulong lastProcessed = 0UL;
bool processedAll = false;
bool composite = CheckDivisors(
prime,
allowedMax,
out lastProcessed,
out processedAll,
out processedCount,
timeLimit);
out processedCount);
if (processedCount > 0UL)
{
@@ -198,38 +163,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
ulong allowedMax,
out ulong lastProcessed,
out bool processedAll,
out ulong processedCount,
TimeSpan? timeLimit)
{
if (_deltasDevice == ByDivisorDeltasDevice.Cpu)
{
return CheckDivisorsCpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit);
}
return CheckDivisorsGpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit);
}
private bool CheckDivisorsCpu(
ulong prime,
ulong allowedMax,
out ulong lastProcessed,
out bool processedAll,
out ulong processedCount,
TimeSpan? timeLimit)
out ulong processedCount)
{
lastProcessed = 0UL;
processedCount = 0UL;
processedAll = false;
PrimeTestTimeLimit limitGuard;
if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard))
{
return false;
}
bool enforceLimit = limitGuard.IsActive;
ulong iterationCounter = 0UL;
if (allowedMax < 3UL)
{
return false;
@@ -273,13 +212,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
while (divisor <= limit)
{
if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired())
{
processedAll = false;
return false;
}
iterationCounter++;
ulong candidate = (ulong)divisor;
processedCount++;
lastProcessed = candidate;
@@ -290,10 +222,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0)
{
// Each divisor appears only once in this monotonically increasing sequence, so
// rebuilding the Montgomery data per candidate avoids caching entries that we
// would never hit again.
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(candidate);
ulong divisorCycle;
if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
@@ -303,9 +232,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
else
{
factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
// EvenPerfectBitScanner reaches this path only for divisors above the precomputed cycle snapshot.
// Those candidates follow q = 2 * prime * k + 1 with prime >= 5, so the downstream order calculator
// can rely on an odd divisor and a prime exponent when deriving the cycle length.
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL)
{
divisorCycle = cycleCache.GetCycleLength(candidate);
@@ -340,248 +266,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
return false;
}
private bool CheckDivisorsGpu(
ulong prime,
ulong allowedMax,
out ulong lastProcessed,
out bool processedAll,
out ulong processedCount,
TimeSpan? timeLimit)
{
lastProcessed = 0UL;
processedCount = 0UL;
processedAll = false;
PrimeTestTimeLimit limitGuard;
if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard))
{
return false;
}
bool enforceLimit = limitGuard.IsActive;
ulong iterationCounter = 0UL;
if (allowedMax < 3UL)
{
return false;
}
GpuUInt128 step128 = new GpuUInt128(prime);
step128.ShiftLeft(1);
if (step128.IsZero)
{
processedAll = true;
return false;
}
GpuUInt128 limit128 = new GpuUInt128(allowedMax);
GpuUInt128 divisor128 = step128;
divisor128.Add(1UL);
if (divisor128.CompareTo(limit128) > 0)
{
processedAll = true;
return false;
}
ulong step = step128.Low;
ulong limit64 = allowedMax;
int batchCapacity = Math.Max(1, _batchSize);
DivisorCycleCache cycleCache = DivisorCycleCache.Shared;
int preferredCycleBatch = cycleCache.PreferredBatchSize;
if (preferredCycleBatch > 0 && preferredCycleBatch < batchCapacity)
{
batchCapacity = preferredCycleBatch;
}
var candidateEvaluator = new MersenneNumberDivisorCandidateGpuEvaluator(batchCapacity);
ulong[] candidateBuffer = ArrayPool<ulong>.Shared.Rent(batchCapacity);
byte[] maskBuffer = ArrayPool<byte>.Shared.Rent(batchCapacity);
Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
const int RemainderTableSize = 6;
byte[] remainderBuffer = ArrayPool<byte>.Shared.Rent(RemainderTableSize);
byte[] stepBuffer = ArrayPool<byte>.Shared.Rent(RemainderTableSize);
Span<byte> remainderSpan = remainderBuffer.AsSpan(0, RemainderTableSize);
Span<byte> stepSpan = stepBuffer.AsSpan(0, RemainderTableSize);
// The remainder and step tables follow the fixed modulus order: 10, 8, 5, 3, 7, 11.
ulong startDivisor64 = divisor128.Low;
remainderSpan[0] = (byte)(startDivisor64 % 10UL);
remainderSpan[1] = (byte)(startDivisor64 % 8UL);
remainderSpan[2] = (byte)(startDivisor64 % 5UL);
remainderSpan[3] = (byte)(startDivisor64 % 3UL);
remainderSpan[4] = (byte)(startDivisor64 % 7UL);
remainderSpan[5] = (byte)(startDivisor64 % 11UL);
stepSpan[0] = (byte)(step % 10UL);
stepSpan[1] = (byte)(step % 8UL);
stepSpan[2] = (byte)(step % 5UL);
stepSpan[3] = (byte)(step % 3UL);
stepSpan[4] = (byte)(step % 7UL);
stepSpan[5] = (byte)(step % 11UL);
bool lastIsSeven = (prime & 3UL) == 3UL;
MersenneNumberDivisorRemainderGpuStepper remainderStepper = new(stepSpan);
MersenneNumberDivisorMontgomeryGpuBuilder? montgomeryBuilder = _montgomeryDevice == ByDivisorMontgomeryDevice.Gpu
? new MersenneNumberDivisorMontgomeryGpuBuilder(batchCapacity)
: null;
MontgomeryDivisorData[]? montgomeryBuffer = montgomeryBuilder != null
? ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity)
: null;
ulong remainingCandidates = ComputeRemainingCandidateCount(startDivisor64, limit64, step);
try
{
while (remainingCandidates > 0UL && divisor128.CompareTo(limit128) <= 0)
{
if (divisor128.High != 0UL)
{
break;
}
ulong currentDivisor = divisor128.Low;
int chunkCount = remainingCandidates > (ulong)batchCapacity
? batchCapacity
: (int)remainingCandidates;
if (chunkCount <= 0)
{
break;
}
Span<ulong> candidateSpan = candidateBuffer.AsSpan(0, chunkCount);
Span<byte> maskSpan = maskBuffer.AsSpan(0, chunkCount);
candidateEvaluator.EvaluateCandidates(
currentDivisor,
step,
limit64,
remainderSpan,
stepSpan,
lastIsSeven,
candidateSpan,
maskSpan);
Span<MontgomeryDivisorData> montgomerySpan = default;
if (montgomeryBuilder != null && montgomeryBuffer != null)
{
montgomerySpan = montgomeryBuffer.AsSpan(0, chunkCount);
// Every candidate divisor is unique across the scan, so build Montgomery data
// for the current batch and drop it immediately instead of maintaining a cache.
montgomeryBuilder.Build(candidateSpan, montgomerySpan);
}
for (int i = 0; i < chunkCount; i++)
{
if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired())
{
processedAll = false;
return false;
}
iterationCounter++;
ulong candidate = candidateSpan[i];
processedCount++;
lastProcessed = candidate;
if (maskSpan[i] == 0)
{
continue;
}
MontgomeryDivisorData divisorData;
if (!montgomerySpan.IsEmpty)
{
divisorData = montgomerySpan[i];
if (divisorData.Modulus != candidate)
{
// The GPU builder might skip entries when masks zero them out; recompute on the CPU
// without caching because divisors never repeat within a session.
divisorData = MontgomeryDivisorData.FromModulus(candidate);
}
}
else
{
// Candidates remain unique even when GPU Montgomery support is disabled, so rebuild
// the divisor data inline instead of persisting unused cache entries.
divisorData = MontgomeryDivisorData.FromModulus(candidate);
}
ulong divisorCycle;
if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
{
divisorCycle = cycleCache.GetCycleLength(candidate);
}
else
{
factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL)
{
divisorCycle = cycleCache.GetCycleLength(candidate);
}
}
if (divisorCycle == prime)
{
processedAll = true;
return true;
}
if (divisorCycle == 0UL)
{
Console.WriteLine($"Divisor cycle was not calculated for {prime}");
}
}
GpuUInt128 increment = step128;
increment.Mul((ulong)chunkCount);
divisor128.Add(increment);
remainderStepper.Advance(chunkCount, remainderSpan);
remainingCandidates -= (ulong)chunkCount;
}
}
finally
{
montgomeryBuilder?.Dispose();
if (montgomeryBuffer != null)
{
ArrayPool<MontgomeryDivisorData>.Shared.Return(montgomeryBuffer, clearArray: false);
}
remainderStepper.Dispose();
ArrayPool<byte>.Shared.Return(remainderBuffer, clearArray: false);
ArrayPool<byte>.Shared.Return(stepBuffer, clearArray: false);
ArrayPool<ulong>.Shared.Return(candidateBuffer, clearArray: false);
ArrayPool<byte>.Shared.Return(maskBuffer, clearArray: false);
candidateEvaluator.Dispose();
}
processedAll = remainingCandidates == 0UL && divisor128.CompareTo(limit128) > 0;
return false;
}
private static ulong ComputeRemainingCandidateCount(ulong currentDivisor, ulong limit, ulong step)
{
if (step == 0UL || currentDivisor > limit)
{
return 0UL;
}
ulong distance = limit - currentDivisor;
ulong stepsAvailable = distance / step;
return stepsAvailable + 1UL;
}
private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
{
ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData);
return residue == 1UL ? (byte)1 : (byte)0;
}
@@ -611,9 +298,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
ulong total = _lastStatusDivisor + processedCount;
// TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
// status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
_lastStatusDivisor = total % interval;
if (total >= interval)
{
do
{
total -= interval;
}
while (total >= interval);
}
_lastStatusDivisor = total;
}
private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
@@ -650,11 +344,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
{
private readonly MersenneNumberDivisorByDivisorCpuTester _owner;
private bool _disposed;
private ulong[]? _primeDeltas;
private ulong[]? _cycleRemainders;
private ulong[]? _residues;
private int _bufferCapacity;
private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator;
internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner)
{
@@ -679,19 +368,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
return;
}
MontgomeryDivisorData effectiveDivisorData = divisorData;
if (effectiveDivisorData.Modulus != divisor)
MontgomeryDivisorData cachedData = divisorData;
if (cachedData.Modulus != divisor)
{
// Each divisor surfaces only once across the monotonic scan, so rebuild the Montgomery data
// instead of retaining a cache entry that would never be reused within future checks.
effectiveDivisorData = MontgomeryDivisorData.FromModulus(divisor);
}
ulong modulus = effectiveDivisorData.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
hits.Clear();
return;
cachedData = MontgomeryDivisorDataCache.Get(divisor);
}
if (divisorCycle == 0UL)
@@ -704,24 +384,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
}
if (_owner._deltasDevice == ByDivisorDeltasDevice.Cpu)
{
CheckDivisorCpu(effectiveDivisorData, divisorCycle, primes, hits);
return;
}
CheckDivisorGpu(effectiveDivisorData, divisorCycle, primes, hits, length);
}
private static void CheckDivisorCpu(
in MontgomeryDivisorData divisorData,
ulong divisorCycle,
ReadOnlySpan<ulong> primes,
Span<byte> hits)
{
// Keep these remainder steppers in place so future updates continue reusing the previously computed residues.
// They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors.
var exponentStepper = new ExponentRemainderStepper(divisorData);
var exponentStepper = new ExponentRemainderStepper(cachedData);
if (!exponentStepper.IsValidModulus)
{
hits.Clear();
@@ -735,7 +400,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0)
: (byte)0;
int length = primes.Length;
for (int i = 1; i < length; i++)
{
remainder = cycleStepper.ComputeNext(primes[i]);
@@ -749,171 +413,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
}
private void CheckDivisorGpu(
in MontgomeryDivisorData divisorData,
ulong divisorCycle,
ReadOnlySpan<ulong> primes,
Span<byte> hits,
int length)
{
EnsureCapacity(length);
Span<ulong> deltaSpan = _primeDeltas!.AsSpan(0, length);
ComputePrimeDeltas(primes, deltaSpan);
Span<ulong> remainderSpan = _cycleRemainders!.AsSpan(0, length);
ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan);
EnsureGpuEvaluator();
Span<ulong> residueSpan = _residues!.AsSpan(0, length);
_gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan);
for (int i = 0; i < length; i++)
{
hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == 1UL ? (byte)1 : (byte)0;
}
}
private void EnsureGpuEvaluator()
{
if (_gpuResidueEvaluator is null)
{
_gpuResidueEvaluator = new MersenneNumberDivisorResidueGpuEvaluator(_owner._batchSize);
}
}
private void ReleaseGpuEvaluator()
{
if (_gpuResidueEvaluator is null)
{
return;
}
_gpuResidueEvaluator.Dispose();
_gpuResidueEvaluator = null;
}
private void EnsureCapacity(int requiredLength)
{
if (requiredLength <= 0)
{
return;
}
if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null)
{
return;
}
int newCapacity = _bufferCapacity;
if (newCapacity == 0)
{
newCapacity = 1;
}
while (newCapacity < requiredLength)
{
newCapacity <<= 1;
}
ReturnBuffers();
_primeDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
_cycleRemainders = ArrayPool<ulong>.Shared.Rent(newCapacity);
_residues = ArrayPool<ulong>.Shared.Rent(newCapacity);
_bufferCapacity = newCapacity;
}
private void ReturnBuffers()
{
if (_primeDeltas is not null)
{
ArrayPool<ulong>.Shared.Return(_primeDeltas, clearArray: false);
_primeDeltas = null;
}
if (_cycleRemainders is not null)
{
ArrayPool<ulong>.Shared.Return(_cycleRemainders, clearArray: false);
_cycleRemainders = null;
}
if (_residues is not null)
{
ArrayPool<ulong>.Shared.Return(_residues, clearArray: false);
_residues = null;
}
_bufferCapacity = 0;
}
private static void ComputePrimeDeltas(ReadOnlySpan<ulong> primes, Span<ulong> deltas)
{
if (primes.Length == 0)
{
return;
}
ulong previous = primes[0];
deltas[0] = previous;
for (int i = 1; i < primes.Length; i++)
{
ulong current = primes[i];
if (current <= previous)
{
throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing.");
}
deltas[i] = current - previous;
previous = current;
}
}
private static void ComputeCycleRemainders(ulong cycleLength, ReadOnlySpan<ulong> primes, ReadOnlySpan<ulong> deltas, Span<ulong> remainders)
{
if (cycleLength == 0UL || primes.Length == 0)
{
remainders.Clear();
return;
}
ulong remainder = primes[0] % cycleLength;
remainders[0] = remainder;
for (int i = 1; i < primes.Length; i++)
{
remainder = AddMod(remainder, deltas[i], cycleLength);
remainders[i] = remainder;
}
}
private static ulong AddMod(ulong value, ulong delta, ulong modulus)
{
if (modulus == 0UL)
{
return 0UL;
}
if (ulong.MaxValue - value < delta)
{
UInt128 extended = (UInt128)value + delta;
return (ulong)(extended % modulus);
}
ulong sum = value + delta;
if (sum >= modulus)
{
sum -= modulus;
if (sum >= modulus)
{
sum %= modulus;
}
}
return sum;
}
public void Dispose()
{
if (_disposed)
@@ -922,7 +421,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
_disposed = true;
ReleaseGpuEvaluator();
_owner.ReturnSession(this);
}
}
@@ -1,49 +0,0 @@
using System;
using System.Threading;
namespace PerfectNumbers.Core;
internal static class DeterministicRandom
{
[ThreadStatic]
private static ulong s_threadState;
private static long s_seedCounter;
public static ulong NextUInt64()
{
ulong state = s_threadState;
if (state == 0UL)
{
state = InitializeThreadState();
}
state ^= state >> 12;
state ^= state << 25;
state ^= state >> 27;
s_threadState = state;
return state * 2685821657736338717UL;
}
public static UInt128 NextUInt128()
{
UInt128 high = (UInt128)NextUInt64();
UInt128 low = (UInt128)NextUInt64();
return (high << 64) | low;
}
private static ulong InitializeThreadState()
{
ulong seed = (ulong)Interlocked.Increment(ref s_seedCounter);
seed += 0x9E3779B97F4A7C15UL;
seed = (seed ^ (seed >> 30)) * 0xBF58476D1CE4E5B9UL;
seed = (seed ^ (seed >> 27)) * 0x94D049BB133111EBUL;
seed ^= seed >> 31;
if (seed == 0UL)
{
seed = 0x9E3779B97F4A7C15UL;
}
return seed;
}
}
+14 -44
View File
@@ -63,9 +63,7 @@ public sealed class DivisorCycleCache
int[]? rentedMissing = null;
Span<int> missingBuffer = length <= StackBufferThreshold
? stackalloc int[length]
: length < PerfectNumberConstants.PooledArrayThreshold
? new Span<int>(new int[length])
: new Span<int>(rentedMissing = ThreadLocalArrayPool<int>.Shared.Rent(length), 0, length);
: new Span<int>(rentedMissing = ArrayPool<int>.Shared.Rent(length), 0, length);
int missingCount = 0;
@@ -108,7 +106,7 @@ public sealed class DivisorCycleCache
if (rentedMissing is not null)
{
ThreadLocalArrayPool<int>.Shared.Return(rentedMissing, clearArray: false);
ArrayPool<int>.Shared.Return(rentedMissing, clearArray: false);
}
}
@@ -130,16 +128,12 @@ public sealed class DivisorCycleCache
ulong[]? rentedDivisors = null;
Span<ulong> workDivisors = count <= StackBufferThreshold
? stackalloc ulong[count]
: count < PerfectNumberConstants.PooledArrayThreshold
? new Span<ulong>(new ulong[count])
: new Span<ulong>(rentedDivisors = ThreadLocalArrayPool<ulong>.Shared.Rent(count), 0, count);
: new Span<ulong>(rentedDivisors = ArrayPool<ulong>.Shared.Rent(count), 0, count);
ulong[]? rentedResults = null;
Span<ulong> workResults = count <= StackBufferThreshold
? stackalloc ulong[count]
: count < PerfectNumberConstants.PooledArrayThreshold
? new Span<ulong>(new ulong[count])
: new Span<ulong>(rentedResults = ThreadLocalArrayPool<ulong>.Shared.Rent(count), 0, count);
: new Span<ulong>(rentedResults = ArrayPool<ulong>.Shared.Rent(count), 0, count);
for (int i = 0; i < count; i++)
{
@@ -155,12 +149,8 @@ public sealed class DivisorCycleCache
if (rentedDivisors is not null)
{
ThreadLocalArrayPool<ulong>.Shared.Return(rentedDivisors, clearArray: false);
}
if (rentedResults is not null)
{
ThreadLocalArrayPool<ulong>.Shared.Return(rentedResults, clearArray: false);
ArrayPool<ulong>.Shared.Return(rentedDivisors, clearArray: false);
ArrayPool<ulong>.Shared.Return(rentedResults!, clearArray: false);
}
}
@@ -182,30 +172,22 @@ public sealed class DivisorCycleCache
ulong[]? rentedPow = null;
Span<ulong> powSpan = length <= StackBufferThreshold
? stackalloc ulong[length]
: length < PerfectNumberConstants.PooledArrayThreshold
? new Span<ulong>(new ulong[length])
: new Span<ulong>(rentedPow = ThreadLocalArrayPool<ulong>.Shared.Rent(length), 0, length);
: new Span<ulong>(rentedPow = ArrayPool<ulong>.Shared.Rent(length), 0, length);
ulong[]? rentedOrder = null;
Span<ulong> orderSpan = length <= StackBufferThreshold
? stackalloc ulong[length]
: length < PerfectNumberConstants.PooledArrayThreshold
? new Span<ulong>(new ulong[length])
: new Span<ulong>(rentedOrder = ThreadLocalArrayPool<ulong>.Shared.Rent(length), 0, length);
: new Span<ulong>(rentedOrder = ArrayPool<ulong>.Shared.Rent(length), 0, length);
ulong[]? rentedResult = null;
Span<ulong> resultSpan = length <= StackBufferThreshold
? stackalloc ulong[length]
: length < PerfectNumberConstants.PooledArrayThreshold
? new Span<ulong>(new ulong[length])
: new Span<ulong>(rentedResult = ThreadLocalArrayPool<ulong>.Shared.Rent(length), 0, length);
: new Span<ulong>(rentedResult = ArrayPool<ulong>.Shared.Rent(length), 0, length);
byte[]? rentedStatus = null;
Span<byte> statusSpan = length <= StackBufferThreshold
? stackalloc byte[length]
: length < PerfectNumberConstants.PooledArrayThreshold
? new Span<byte>(new byte[length])
: new Span<byte>(rentedStatus = ThreadLocalArrayPool<byte>.Shared.Rent(length), 0, length);
: new Span<byte>(rentedStatus = ArrayPool<byte>.Shared.Rent(length), 0, length);
ref ulong divisorRef = ref MemoryMarshal.GetReference(divisors);
divisorBuffer.View.CopyFromCPU(ref divisorRef, length);
@@ -252,22 +234,10 @@ public sealed class DivisorCycleCache
resultSpan.CopyTo(destination);
if (rentedPow is not null)
{
ThreadLocalArrayPool<ulong>.Shared.Return(rentedPow, clearArray: false);
}
if (rentedOrder is not null)
{
ThreadLocalArrayPool<ulong>.Shared.Return(rentedOrder, clearArray: false);
}
if (rentedResult is not null)
{
ThreadLocalArrayPool<ulong>.Shared.Return(rentedResult, clearArray: false);
}
if (rentedStatus is not null)
{
ThreadLocalArrayPool<byte>.Shared.Return(rentedStatus, clearArray: false);
ArrayPool<ulong>.Shared.Return(rentedPow, clearArray: false);
ArrayPool<ulong>.Shared.Return(rentedOrder!, clearArray: false);
ArrayPool<ulong>.Shared.Return(rentedResult!, clearArray: false);
ArrayPool<byte>.Shared.Return(rentedStatus!, clearArray: false);
}
execution.Dispose();
@@ -47,7 +47,7 @@ internal struct ExponentRemainderStepper
// TODO: Route these state resets through the ProcessEightBitWindows helper once the scalar
// Pow2MontgomeryMod implementation adopts it so delta stepping inherits the benchmarked
// 2× gains for large exponents.
_currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
_previousExponent = exponent;
_hasState = true;
return ReduceCurrent();
@@ -59,7 +59,7 @@ internal struct ExponentRemainderStepper
// cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
// TODO: Replace this per-delta powmod with the upcoming windowed ladder so incremental
// updates stop paying the single-bit cost highlighted in GpuPow2ModBenchmarks.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
_previousExponent = exponent;
return ReduceCurrent();
@@ -77,7 +77,7 @@ internal struct ExponentRemainderStepper
{
// TODO: Switch this reload to the shared windowed pow2 helper once available so CPU
// residue checks align with the optimized ProcessEightBitWindows timings.
_currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
_previousExponent = exponent;
_hasState = true;
return _currentMontgomery == _montgomeryOne;
@@ -88,7 +88,7 @@ internal struct ExponentRemainderStepper
// lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor.
// TODO: Use the windowed delta pow2 helper here as well to avoid the single-bit ladder that
// currently lags behind the benchmarked implementation for huge divisor cycles.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
_previousExponent = exponent;
return _currentMontgomery == _montgomeryOne;
@@ -123,7 +123,7 @@ internal struct ExponentRemainderStepper
{
// TODO: Replace this fallback path with the upcoming ProcessEightBitWindows helper so
// fresh Montgomery states also benefit from the faster pow2 ladder measured on CPUs.
_currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
_previousExponent = exponent;
_hasState = true;
isUnity = _currentMontgomery == _montgomeryOne;
-19
View File
@@ -1,19 +0,0 @@
namespace PerfectNumbers.Core
{
public readonly struct FactorEntry
{
public readonly ulong Value;
public readonly int Exponent = 1;
public FactorEntry(ulong value, int exponent)
{
Value = value;
Exponent = exponent;
}
public FactorEntry(ulong value)
{
Value = value;
}
}
}
+132 -135
View File
@@ -6,159 +6,156 @@ namespace PerfectNumbers.Core.Gpu;
public static class GpuContextPool
{
private static readonly bool PoolingEnabled = true;
// TODO: Introduce accelerator-specific warmup so pooled contexts precompile the ProcessEightBitWindows kernels and load
// divisor-cycle data before the first scan begins.
// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
public static bool ForceCpu { get; set; } = false;
private static readonly bool PoolingEnabled = true;
// TODO: Introduce accelerator-specific warmup so pooled contexts precompile the ProcessEightBitWindows kernels and load
// divisor-cycle data before the first scan begins.
// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
public static bool ForceCpu { get; set; } = false;
internal sealed class PooledContext
{
private bool disposed;
public Context Context { get; }
public Accelerator Accelerator { get; }
public bool IsCpu { get; }
public object ExecutionLock { get; } = new();
public PooledContext(bool preferCpu)
internal sealed class PooledContext
{
Context = Context.CreateDefault();
Accelerator = Context.GetPreferredDevice(preferCpu).CreateAccelerator(Context);
IsCpu = Accelerator.AcceleratorType == AcceleratorType.CPU;
// NOTE: Avoid loading/compiling any kernel here to prevent implicit
// CL stream/queue creation during accelerator construction.
// Some OpenCL drivers are fragile when a queue is created immediately
// after device init. Real kernels will be JIT-loaded at first use.
private bool disposed;
public Context Context { get; }
public Accelerator Accelerator { get; }
public bool IsCpu { get; }
public object ExecutionLock { get; } = new();
public PooledContext(bool preferCpu)
{
Context = Context.CreateDefault();
Accelerator = Context.GetPreferredDevice(preferCpu).CreateAccelerator(Context);
IsCpu = Accelerator.AcceleratorType == AcceleratorType.CPU;
// NOTE: Avoid loading/compiling any kernel here to prevent implicit
// CL stream/queue creation during accelerator construction.
// Some OpenCL drivers are fragile when a queue is created immediately
// after device init. Real kernels will be JIT-loaded at first use.
}
public void Dispose()
{
if (disposed)
{
return;
}
// Ensure all cached GPU buffers for this accelerator are released.
NttGpuMath.ClearCaches(Accelerator);
// Release PrimeTester GPU state for this accelerator (kernel/device primes)
// Clear any per-accelerator cached resources to avoid releasing
// after the accelerator is destroyed.
PrimeTester.ClearGpuCaches(Accelerator);
Accelerator.Dispose();
Context.Dispose();
disposed = true;
}
}
private static readonly ConcurrentQueue<PooledContext> CpuPool = new();
private static readonly ConcurrentQueue<PooledContext> GpuPool = new();
private static readonly object CreationLock = new();
public static GpuContextLease Rent()
{
return RentPreferred(ForceCpu);
}
// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
public static GpuContextLease RentPreferred(bool preferCpu)
{
if (PoolingEnabled)
{
if (preferCpu && CpuPool.TryDequeue(out var cpu))
{
return new GpuContextLease(cpu, CreationLock);
}
if (!preferCpu && GpuPool.TryDequeue(out var gpu))
{
return new GpuContextLease(gpu, CreationLock);
}
}
// Serialize accelerator creation to avoid concurrent CL device/stream setup across threads.
lock (CreationLock)
{
return new GpuContextLease(new PooledContext(preferCpu), CreationLock);
}
}
public static void DisposeAll()
{
if (!PoolingEnabled)
{
return;
}
while (CpuPool.TryDequeue(out var cpu))
{
cpu.Dispose();
}
while (GpuPool.TryDequeue(out var gpu))
{
gpu.Dispose();
}
}
private static void Return(PooledContext ctx)
{
if (PoolingEnabled)
{
ctx.Accelerator.Synchronize();
if (ctx.IsCpu)
{
CpuPool.Enqueue(ctx);
}
else
{
GpuPool.Enqueue(ctx);
}
}
else
{
lock (CreationLock)
{
ctx.Dispose();
}
}
}
public void Dispose()
{
if (disposed)
{
return;
}
// Ensure all cached GPU buffers for this accelerator are released.
NttGpuMath.ClearCaches(Accelerator);
// Release PrimeTester GPU state for this accelerator (kernel/device primes)
// Clear any per-accelerator cached resources to avoid releasing
// after the accelerator is destroyed.
PrimeTester.ClearGpuCaches(Accelerator);
Accelerator.Dispose();
Context.Dispose();
disposed = true;
}
}
private static readonly ConcurrentQueue<PooledContext> CpuPool = new();
private static readonly ConcurrentQueue<PooledContext> GpuPool = new();
private static readonly object CreationLock = new();
public static GpuContextLease Rent()
{
return RentPreferred(ForceCpu);
}
// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
public static GpuContextLease RentPreferred(bool preferCpu)
{
if (PoolingEnabled)
{
if (preferCpu && CpuPool.TryDequeue(out var cpu))
{
return new GpuContextLease(cpu, CreationLock);
}
if (!preferCpu && GpuPool.TryDequeue(out var gpu))
{
return new GpuContextLease(gpu, CreationLock);
}
}
// Serialize accelerator creation to avoid concurrent CL device/stream setup across threads.
lock (CreationLock)
{
return new GpuContextLease(new PooledContext(preferCpu), CreationLock);
}
}
public static void DisposeAll()
{
if (!PoolingEnabled)
{
return;
}
while (CpuPool.TryDequeue(out var cpu))
{
cpu.Dispose();
}
while (GpuPool.TryDequeue(out var gpu))
{
gpu.Dispose();
}
}
private static void Return(PooledContext ctx)
{
if (PoolingEnabled)
{
ctx.Accelerator.Synchronize();
if (ctx.IsCpu)
{
CpuPool.Enqueue(ctx);
}
else
{
GpuPool.Enqueue(ctx);
}
}
else
{
lock (CreationLock)
{
ctx.Dispose();
}
}
}
public struct GpuContextLease
public struct GpuContextLease : IDisposable
{
private readonly PooledContext _ctx;
private bool disposed;
private bool disposed;
internal GpuContextLease(PooledContext ctx, object kernelInitLock)
internal GpuContextLease(PooledContext ctx, object kernelInitLock)
{
_ctx = ctx;
KernelInitLock = kernelInitLock;
}
KernelInitLock = kernelInitLock;
}
public Context Context => _ctx.Context;
public Accelerator Accelerator => _ctx.Accelerator;
public object KernelInitLock { get; }
public object KernelInitLock { get; }
public object ExecutionLock => _ctx.ExecutionLock;
public object ExecutionLock => _ctx.ExecutionLock;
internal void ReturnToPool()
{
if (disposed)
{
return;
}
Return(_ctx);
disposed = true;
}
public void Dispose()
{
if (disposed)
{
return;
}
public void Dispose()
{
ReturnToPool();
Return(_ctx);
disposed = true;
}
}
}
+31 -88
View File
@@ -37,7 +37,6 @@ public sealed class KernelContainer
ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? IncrementalOrder;
public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2ModOrder;
public Action<AcceleratorStream, Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int>? MontgomeryOddPowers;
// Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallCycles;
@@ -110,15 +109,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
get
{
var accel = Accelerator; // avoid capturing 'this' in lambda
// lock (_gpu!.Value.KernelInitLock)
// {
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Order, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>>();
});
// }
}
}
}
@@ -127,15 +126,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
get
{
var accel = Accelerator;
// lock (_gpu!.Value.KernelInitLock)
// {
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
// }
}
}
}
@@ -144,15 +143,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
get
{
var accel = Accelerator;
// lock (_gpu!.Value.KernelInitLock)
// {
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Incremental, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
// }
}
}
}
@@ -161,15 +160,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
get
{
var accel = Accelerator;
// lock (_gpu!.Value.KernelInitLock)
// {
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
// }
}
}
}
@@ -178,66 +177,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
get
{
var accel = Accelerator;
// lock (_gpu!.Value.KernelInitLock)
// {
lock (_gpu!.Value.KernelInitLock)
{
return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModOrderKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
});
// }
}
}
public readonly Action<AcceleratorStream, Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int> MontgomeryOddPowerKernel
{
get
{
var accel = Accelerator;
// lock (_gpu!.Value.KernelInitLock)
// {
return KernelContainer.InitOnce(ref _kernels.MontgomeryOddPowers, () =>
{
var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int>(MontgomeryOddPowerKernelScan);
var kernel = KernelUtil.GetKernel(loaded);
return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int>>();
});
// }
}
}
private static void MontgomeryOddPowerKernelScan(
Index1D index,
ArrayView1D<MontgomeryDivisorData, Stride1D.Dense> divisors,
ArrayView1D<int, Stride1D.Dense> counts,
ArrayView1D<ulong, Stride1D.Dense> results,
int stride)
{
int tableIndex = index.X;
int baseOffset = tableIndex * stride;
MontgomeryDivisorData divisor = divisors[tableIndex];
int oddPowerCount = counts[tableIndex];
if (oddPowerCount <= 0)
{
return;
}
ulong montgomeryTwo = divisor.MontgomeryTwo;
results[baseOffset] = montgomeryTwo;
if (oddPowerCount == 1)
{
return;
}
ulong modulus = divisor.Modulus;
ulong nPrime = divisor.NPrime;
ulong square = divisor.MontgomeryTwoSquared;
ulong current = montgomeryTwo;
for (int i = 1; i < oddPowerCount; i++)
{
current = ULongExtensions.MontgomeryMultiply(current, square, modulus, nPrime);
results[baseOffset + i] = current;
}
}
}
@@ -396,10 +344,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
_stream?.Dispose();
_stream = null;
if (_gpu.HasValue)
{
_gpu.Value.ReturnToPool();
}
_gpu?.Dispose();
_gpu = null;
_executionLock = null;
@@ -732,25 +677,25 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
// GC.SuppressFinalize(this);
}
public readonly struct ExecutionScope
public readonly struct ExecutionScope : IDisposable
{
private readonly object? _lock;
public ExecutionScope(object? sync)
{
_lock = sync;
// if (_lock is not null)
// {
// Monitor.Enter(_lock);
// }
if (_lock is not null)
{
Monitor.Enter(_lock);
}
}
public void Dispose()
{
// if (_lock is not null)
// {
// Monitor.Exit(_lock);
// }
if (_lock is not null)
{
Monitor.Exit(_lock);
}
}
}
}
@@ -759,8 +704,6 @@ public class GpuKernelPool
{
private static readonly ConcurrentDictionary<Accelerator, KernelContainer> KernelCache = new(); // TODO: Replace this concurrent map with a simple accelerator-indexed lookup once kernel launchers are prewarmed during startup so we can drop the thread-safe wrapper entirely.
public static bool CpuCycleUsesGpu { get; set; }
private static KernelContainer GetKernels(Accelerator accelerator)
{
return KernelCache.GetOrAdd(accelerator, _ => new KernelContainer());
@@ -777,8 +720,8 @@ public class GpuKernelPool
}
// Ensure single upload per accelerator even if multiple threads race here.
// lock (kernels) // TODO: Remove this lock by pre-uploading the immutable small-cycle snapshot during initialization; once no mutation happens at runtime, the pool must expose a simple reference without synchronization.
// {
lock (kernels) // TODO: Remove this lock by pre-uploading the immutable small-cycle snapshot during initialization; once no mutation happens at runtime, the pool must expose a simple reference without synchronization.
{
if (kernels.SmallCycles is { } existing)
{
return existing.View;
@@ -789,7 +732,7 @@ public class GpuKernelPool
device.View.CopyFromCPU(host);
kernels.SmallCycles = device;
return device.View;
// }
}
}
public static ResiduePrimeViews EnsureSmallPrimesOnDevice(Accelerator accelerator)
@@ -803,8 +746,8 @@ public class GpuKernelPool
return new ResiduePrimeViews(lastOne.View, lastSeven.View, lastOnePow2.View, lastSevenPow2.View);
}
// lock (kernels) // TODO: Inline these small-prime uploads into startup initialization alongside the small-cycle snapshot so we can drop runtime locking and keep the GPU pool free of synchronization.
// {
lock (kernels) // TODO: Inline these small-prime uploads into startup initialization alongside the small-cycle snapshot so we can drop runtime locking and keep the GPU pool free of synchronization.
{
if (kernels.SmallPrimesLastOne is { } existingLastOne &&
kernels.SmallPrimesLastSeven is { } existingLastSeven &&
kernels.SmallPrimesPow2LastOne is { } existingLastOnePow2 &&
@@ -833,7 +776,7 @@ public class GpuKernelPool
kernels.SmallPrimesPow2LastSeven = deviceLastSevenPow2;
return new ResiduePrimeViews(deviceLastOne.View, deviceLastSeven.View, deviceLastOnePow2.View, deviceLastSevenPow2.View);
// }
}
}
public static GpuKernelLease GetKernel(bool useGpuOrder)
+12 -12
View File
@@ -267,7 +267,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void Add(in ReadOnlyGpuUInt128 other) => AddInternal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Add(GpuUInt128 other) => AddInternal(other.High, other.Low);
public void Add(in GpuUInt128 other) => AddInternal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void AddInternal(ulong otherHigh, ulong otherLow)
@@ -294,17 +294,17 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void AddMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void AddMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
public void AddMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void AddMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus)
public void AddMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
AddModInternal(value.High, value.Low, in readOnlyMod);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void AddMod(GpuUInt128 value, GpuUInt128 modulus)
public void AddMod(in GpuUInt128 value, in GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
AddModInternal(value.High, value.Low, in readOnlyMod);
@@ -390,7 +390,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void AddMod(GpuUInt128 value, ulong modulus) => AddMod(value.Low, modulus);
public void AddMod(in GpuUInt128 value, ulong modulus) => AddMod(value.Low, modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void AddMod(ulong value, ulong modulus)
@@ -414,7 +414,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void Sub(in ReadOnlyGpuUInt128 other) => SubInternal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Sub(GpuUInt128 other) => SubInternal(other.High, other.Low);
public void Sub(in GpuUInt128 other) => SubInternal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void SubInternal(ulong otherHigh, ulong otherLow)
@@ -461,7 +461,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void Mul64(in ReadOnlyGpuUInt128 other) => Mul64Internal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Mul64(GpuUInt128 other) => Mul64Internal(other.High, other.Low);
public void Mul64(in GpuUInt128 other) => Mul64Internal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void Mul64Internal(ulong otherHigh, ulong otherLow)
@@ -588,7 +588,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int TrailingZeroCount(GpuUInt128 value)
public static int TrailingZeroCount(in GpuUInt128 value)
{
if (value.Low != 0UL)
{
@@ -646,17 +646,17 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
public void SubMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void SubMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
public void SubMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void SubMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus)
public void SubMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
SubModInternal(value.High, value.Low, in readOnlyMod);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void SubMod(GpuUInt128 value, GpuUInt128 modulus)
public void SubMod(in GpuUInt128 value, in GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
SubModInternal(value.High, value.Low, in readOnlyMod);
@@ -727,7 +727,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Mul(GpuUInt128 other)
public void Mul(in GpuUInt128 other)
{
ReadOnlyGpuUInt128 readOnlyThis = AsReadOnly();
ReadOnlyGpuUInt128 readOnlyOther = other.AsReadOnly();
@@ -7,6 +7,7 @@ using System.Numerics;
using ILGPU;
using ILGPU.Runtime;
using MontgomeryDivisorData = PerfectNumbers.Core.MontgomeryDivisorData;
using MontgomeryDivisorDataCache = PerfectNumbers.Core.MontgomeryDivisorDataCache;
using CycleRemainderStepper = PerfectNumbers.Core.CycleRemainderStepper;
namespace PerfectNumbers.Core.Gpu;
@@ -35,7 +36,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeExponent(Accelerator accelerator) =>
_kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeExponentKernel));
_kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeMontgomeryExponentKernel));
private Action<Index1D, ArrayView<ulong>, ArrayView<byte>, byte> GetRemainderDeltaKernel(Accelerator accelerator) =>
_remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<byte>, byte>(ComputeRemainderDeltasKernel));
@@ -71,7 +72,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
}
public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
public bool IsPrime(ulong prime, out bool divisorsExhausted)
{
ulong allowedMax;
int batchCapacity;
@@ -108,7 +109,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
composite = CheckDivisors(
prime,
allowedMax,
timeLimit,
accelerator,
kernel,
resources.DivisorsBuffer,
@@ -176,7 +176,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
private static bool CheckDivisors(
ulong prime,
ulong allowedMax,
TimeSpan? timeLimit,
Accelerator accelerator,
Action<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>> kernel,
MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> divisorsBuffer,
@@ -204,15 +203,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
bool processedAll = false;
processedCount = 0UL;
lastProcessed = 0UL;
coveredRange = false;
PrimeTestTimeLimit limitGuard;
if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard))
{
return false;
}
bool enforceLimit = limitGuard.IsActive;
Span<MontgomeryDivisorData> divisorDataSpan;
Span<ulong> divisorSpan;
@@ -288,12 +278,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
{
while (currentK <= maxK)
{
if (enforceLimit && limitGuard.HasExpired())
{
processedAll = false;
return false;
}
int chunkCount = Math.Min(chunkCapacity, batchCapacity);
ulong remainingK = maxK - currentK + 1UL;
if ((ulong)chunkCount > remainingK)
@@ -412,8 +396,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ulong divisorValue = filteredDivisorsSpan[i];
divisorSpan[i] = divisorValue;
exponentSpan[i] = prime;
// Divisors never repeat while scanning the monotonically increasing sequence, so regenerate
// the Montgomery data for each entry rather than storing cache lines that no future batch would hit.
divisorDataSpan[i] = MontgomeryDivisorData.FromModulus(divisorValue);
}
@@ -569,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
ulong exponent = exponents[index];
hits[index] = exponent.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0;
hits[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
}
public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -777,7 +759,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
mask[globalIndex] = 1;
}
private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
{
ulong modulus = divisor.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -788,7 +770,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ulong exponent = exponents[index];
results[index] = exponent.Pow2ModBinaryGpu(modulus);
results[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: true);
}
public ulong DivisorLimit
@@ -1,262 +0,0 @@
using System;
using System.Buffers;
using System.Runtime.InteropServices;
using ILGPU;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable
{
private const int RemainderCount = 6;
private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator;
private readonly Action<Index1D, CandidateFilterArgs, ArrayView<ulong>, ArrayView<byte>> _kernel;
private MemoryBuffer1D<ulong, Stride1D.Dense> _candidateBuffer;
private MemoryBuffer1D<byte, Stride1D.Dense> _maskBuffer;
private ulong[] _hostCandidates;
private byte[] _hostMask;
private int _capacity;
private bool _disposed;
internal MersenneNumberDivisorCandidateGpuEvaluator(int batchSize)
{
_lease = GpuContextPool.RentPreferred(preferCpu: false);
_accelerator = _lease.Accelerator;
_kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, CandidateFilterArgs, ArrayView<ulong>, ArrayView<byte>>(EvaluateCandidatesKernel);
_capacity = Math.Max(1, batchSize);
_candidateBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_maskBuffer = _accelerator.Allocate1D<byte>(_capacity);
_hostCandidates = ArrayPool<ulong>.Shared.Rent(_capacity);
_hostMask = ArrayPool<byte>.Shared.Rent(_capacity);
}
public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_candidateBuffer.Dispose();
_maskBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostCandidates, clearArray: false);
ArrayPool<byte>.Shared.Return(_hostMask, clearArray: false);
_lease.Dispose();
}
internal void EvaluateCandidates(
ulong startDivisor,
ulong step,
ulong limit,
ReadOnlySpan<byte> remainders,
ReadOnlySpan<byte> steps,
bool lastIsSeven,
Span<ulong> candidates,
Span<byte> mask)
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(MersenneNumberDivisorCandidateGpuEvaluator));
}
if (remainders.Length != RemainderCount || steps.Length != RemainderCount)
{
throw new ArgumentException("Remainder and step tables must contain exactly six entries.");
}
int count = candidates.Length;
if (count == 0)
{
return;
}
EnsureCapacity(count);
Span<ulong> hostCandidateSpan = _hostCandidates.AsSpan(0, count);
Span<byte> hostMaskSpan = _hostMask.AsSpan(0, count);
ArrayView1D<ulong, Stride1D.Dense> candidateView = _candidateBuffer.View.SubView(0, count);
ArrayView1D<byte, Stride1D.Dense> maskView = _maskBuffer.View.SubView(0, count);
var args = new CandidateFilterArgs(
startDivisor,
step,
limit,
remainders[0],
remainders[1],
remainders[2],
remainders[3],
remainders[4],
remainders[5],
steps[0],
steps[1],
steps[2],
steps[3],
steps[4],
steps[5],
lastIsSeven ? (byte)1 : (byte)0);
_kernel(count, args, candidateView, maskView);
candidateView.CopyToCPU(ref MemoryMarshal.GetReference(hostCandidateSpan), count);
maskView.CopyToCPU(ref MemoryMarshal.GetReference(hostMaskSpan), count);
hostCandidateSpan.CopyTo(candidates);
hostMaskSpan.CopyTo(mask);
}
private void EnsureCapacity(int required)
{
if (required <= _capacity)
{
return;
}
int newCapacity = _capacity == 0 ? 1 : _capacity;
while (newCapacity < required)
{
newCapacity <<= 1;
}
_candidateBuffer.Dispose();
_maskBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostCandidates, clearArray: false);
ArrayPool<byte>.Shared.Return(_hostMask, clearArray: false);
_capacity = newCapacity;
_candidateBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_maskBuffer = _accelerator.Allocate1D<byte>(_capacity);
_hostCandidates = ArrayPool<ulong>.Shared.Rent(_capacity);
_hostMask = ArrayPool<byte>.Shared.Rent(_capacity);
}
private static void EvaluateCandidatesKernel(
Index1D index,
CandidateFilterArgs args,
ArrayView<ulong> candidates,
ArrayView<byte> mask)
{
int idx = index;
if (idx >= mask.Length)
{
return;
}
ulong offset = (ulong)idx;
ulong candidate = args.StartDivisor + args.Step * offset;
candidates[idx] = candidate;
if (candidate > args.Limit || (args.Step != 0UL && candidate < args.StartDivisor))
{
mask[idx] = 0;
return;
}
int remainderValue10 = ComputeRemainder(args.Remainder10, args.Step10, idx, 10);
int remainderValue8 = ComputeRemainder(args.Remainder8, args.Step8, idx, 8);
int remainderValue5 = ComputeRemainder(args.Remainder5, args.Step5, idx, 5);
int remainderValue3 = ComputeRemainder(args.Remainder3, args.Step3, idx, 3);
int remainderValue7 = ComputeRemainder(args.Remainder7, args.Step7, idx, 7);
int remainderValue11 = ComputeRemainder(args.Remainder11, args.Step11, idx, 11);
bool lastIsSeven = args.LastIsSevenFlag != 0;
bool admissible = lastIsSeven
? (remainderValue10 == 3 || remainderValue10 == 7 || remainderValue10 == 9)
: (remainderValue10 == 1 || remainderValue10 == 3 || remainderValue10 == 9);
if (!admissible)
{
mask[idx] = 0;
return;
}
if ((remainderValue8 != 1 && remainderValue8 != 7)
|| remainderValue3 == 0
|| remainderValue5 == 0
|| remainderValue7 == 0
|| remainderValue11 == 0)
{
mask[idx] = 0;
return;
}
mask[idx] = 1;
}
private static int ComputeRemainder(byte remainder, byte step, int index, int modulus)
{
if (modulus == 0)
{
return 0;
}
long value = remainder;
value += (long)step * index;
int result = (int)(value % modulus);
if (result < 0)
{
result += modulus;
}
return result;
}
private readonly struct CandidateFilterArgs
{
public CandidateFilterArgs(
ulong startDivisor,
ulong step,
ulong limit,
byte remainder10,
byte remainder8,
byte remainder5,
byte remainder3,
byte remainder7,
byte remainder11,
byte step10,
byte step8,
byte step5,
byte step3,
byte step7,
byte step11,
byte lastIsSevenFlag)
{
StartDivisor = startDivisor;
Step = step;
Limit = limit;
Remainder10 = remainder10;
Remainder8 = remainder8;
Remainder5 = remainder5;
Remainder3 = remainder3;
Remainder7 = remainder7;
Remainder11 = remainder11;
Step10 = step10;
Step8 = step8;
Step5 = step5;
Step3 = step3;
Step7 = step7;
Step11 = step11;
LastIsSevenFlag = lastIsSevenFlag;
}
public ulong StartDivisor { get; }
public ulong Step { get; }
public ulong Limit { get; }
public byte Remainder10 { get; }
public byte Remainder8 { get; }
public byte Remainder5 { get; }
public byte Remainder3 { get; }
public byte Remainder7 { get; }
public byte Remainder11 { get; }
public byte Step10 { get; }
public byte Step8 { get; }
public byte Step5 { get; }
public byte Step3 { get; }
public byte Step7 { get; }
public byte Step11 { get; }
public byte LastIsSevenFlag { get; }
}
}
@@ -1,179 +0,0 @@
using System;
using System.Buffers;
using System.Runtime.InteropServices;
using PerfectNumbers.Core;
using ILGPU;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
{
private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator;
private readonly Action<Index1D, ArrayView<ulong>, ArrayView<MontgomeryDivisorData>> _kernel;
private MemoryBuffer1D<ulong, Stride1D.Dense> _modulusBuffer;
private MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> _resultBuffer;
private ulong[] _hostModulusBuffer;
private MontgomeryDivisorData[] _hostResultBuffer;
private int _capacity;
private bool _disposed;
internal MersenneNumberDivisorMontgomeryGpuBuilder(int batchSize)
{
_lease = GpuContextPool.RentPreferred(preferCpu: false);
_accelerator = _lease.Accelerator;
_kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<MontgomeryDivisorData>>(BuildKernel);
_capacity = Math.Max(1, batchSize);
_modulusBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultBuffer = _accelerator.Allocate1D<MontgomeryDivisorData>(_capacity);
_hostModulusBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_hostResultBuffer = ArrayPool<MontgomeryDivisorData>.Shared.Rent(_capacity);
}
public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_modulusBuffer.Dispose();
_resultBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostModulusBuffer, clearArray: false);
ArrayPool<MontgomeryDivisorData>.Shared.Return(_hostResultBuffer, clearArray: false);
_lease.Dispose();
}
internal void Build(ReadOnlySpan<ulong> moduli, Span<MontgomeryDivisorData> destination)
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(MersenneNumberDivisorMontgomeryGpuBuilder));
}
int length = moduli.Length;
if (length == 0)
{
return;
}
EnsureCapacity(length);
Span<ulong> hostModulusSpan = _hostModulusBuffer.AsSpan(0, length);
moduli.CopyTo(hostModulusSpan);
ArrayView1D<ulong, Stride1D.Dense> modulusView = _modulusBuffer.View.SubView(0, length);
modulusView.CopyFromCPU(ref MemoryMarshal.GetReference(hostModulusSpan), length);
ArrayView1D<MontgomeryDivisorData, Stride1D.Dense> resultView = _resultBuffer.View.SubView(0, length);
_kernel(length, modulusView, resultView);
Span<MontgomeryDivisorData> hostResultSpan = _hostResultBuffer.AsSpan(0, length);
resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostResultSpan), length);
hostResultSpan.CopyTo(destination);
}
private void EnsureCapacity(int required)
{
if (required <= _capacity)
{
return;
}
int newCapacity = _capacity == 0 ? 1 : _capacity;
while (newCapacity < required)
{
newCapacity <<= 1;
}
_modulusBuffer.Dispose();
_resultBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostModulusBuffer, clearArray: false);
ArrayPool<MontgomeryDivisorData>.Shared.Return(_hostResultBuffer, clearArray: false);
_capacity = newCapacity;
_modulusBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultBuffer = _accelerator.Allocate1D<MontgomeryDivisorData>(_capacity);
_hostModulusBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_hostResultBuffer = ArrayPool<MontgomeryDivisorData>.Shared.Rent(_capacity);
}
private static void BuildKernel(Index1D index, ArrayView<ulong> moduli, ArrayView<MontgomeryDivisorData> destination)
{
int idx = index;
if (idx >= moduli.Length)
{
return;
}
ulong modulus = moduli[idx];
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
destination[idx] = default;
return;
}
ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue64(modulus);
ulong montgomeryTwo = montgomeryOne.MulMod64GpuCompatible(2UL, modulus);
ulong montgomeryTwoSquared = MontgomeryMultiplyGpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
destination[idx] = new MontgomeryDivisorData(
modulus,
nPrime,
montgomeryOne,
montgomeryTwo,
montgomeryTwoSquared);
}
private static ulong ComputeMontgomeryResidue64(ulong modulus)
{
ulong remainder = ulong.MaxValue % modulus;
remainder++;
if (remainder >= modulus)
{
remainder -= modulus;
}
return remainder;
}
private static ulong MontgomeryMultiplyGpu(ulong a, ulong b, ulong modulus, ulong nPrime)
{
GpuUInt128 product = new(a);
product.Mul64(new GpuUInt128(b));
ulong tLow = product.Low;
ulong tHigh = product.High;
ulong m = unchecked(tLow * nPrime);
GpuUInt128 mProduct = new(m);
mProduct.Mul64(new GpuUInt128(modulus));
ulong mLow = mProduct.Low;
ulong mHigh = mProduct.High;
ulong carry = unchecked((tLow + mLow) < tLow ? 1UL : 0UL);
ulong result = tHigh + mHigh + carry;
if (result >= modulus)
{
result -= modulus;
}
return result;
}
private static ulong ComputeMontgomeryNPrime(ulong modulus)
{
ulong inv = modulus;
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
inv *= unchecked(2UL - modulus * inv);
return unchecked(0UL - inv);
}
}
@@ -1,140 +0,0 @@
using System;
using System.Buffers;
using System.Runtime.InteropServices;
using ILGPU;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable
{
private const int RemainderCount = 6;
private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator;
private readonly Action<Index1D, int, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>> _kernel;
private MemoryBuffer1D<byte, Stride1D.Dense> _remainderBuffer;
private MemoryBuffer1D<byte, Stride1D.Dense> _stepBuffer;
private MemoryBuffer1D<byte, Stride1D.Dense> _resultBuffer;
private byte[] _hostRemainders;
private byte[] _hostSteps;
private readonly int _length;
private int _capacity;
private bool _disposed;
internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan<byte> steps)
{
if (steps.Length != RemainderCount)
{
throw new ArgumentException("Step table must contain exactly six entries.");
}
_lease = GpuContextPool.RentPreferred(preferCpu: false);
_accelerator = _lease.Accelerator;
_kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, int, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>>(AdvanceRemaindersKernel);
_length = steps.Length;
_capacity = Math.Max(1, _length);
_remainderBuffer = _accelerator.Allocate1D<byte>(_capacity);
_stepBuffer = _accelerator.Allocate1D<byte>(_capacity);
_resultBuffer = _accelerator.Allocate1D<byte>(_capacity);
_hostRemainders = ArrayPool<byte>.Shared.Rent(_capacity);
_hostSteps = ArrayPool<byte>.Shared.Rent(_capacity);
steps.CopyTo(_hostSteps.AsSpan(0, _length));
ArrayView1D<byte, Stride1D.Dense> stepView = _stepBuffer.View.SubView(0, _length);
stepView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostSteps.AsSpan(0, _length)), _length);
}
public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_remainderBuffer.Dispose();
_stepBuffer.Dispose();
_resultBuffer.Dispose();
ArrayPool<byte>.Shared.Return(_hostRemainders, clearArray: false);
ArrayPool<byte>.Shared.Return(_hostSteps, clearArray: false);
_lease.Dispose();
}
internal void Advance(int advanceCount, Span<byte> remainders)
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(MersenneNumberDivisorRemainderGpuStepper));
}
if (advanceCount <= 0 || remainders.Length == 0)
{
return;
}
if (remainders.Length != _length)
{
throw new ArgumentException("Remainder buffer length must match the configured length.");
}
Span<byte> hostRemainderSpan = _hostRemainders.AsSpan(0, _length);
remainders.CopyTo(hostRemainderSpan);
ArrayView1D<byte, Stride1D.Dense> remainderView = _remainderBuffer.View.SubView(0, _length);
remainderView.CopyFromCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length);
ArrayView1D<byte, Stride1D.Dense> stepView = _stepBuffer.View.SubView(0, _length);
ArrayView1D<byte, Stride1D.Dense> resultView = _resultBuffer.View.SubView(0, _length);
_kernel(_length, advanceCount, remainderView, stepView, resultView);
resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length);
hostRemainderSpan.CopyTo(remainders);
}
private static void AdvanceRemaindersKernel(Index1D index, int advanceCount, ArrayView<byte> remainders, ArrayView<byte> steps, ArrayView<byte> results)
{
int idx = index;
if (idx >= remainders.Length)
{
return;
}
int modulus = GetModulus(idx);
if (modulus == 0)
{
results[idx] = remainders[idx];
return;
}
int remainder = remainders[idx] % modulus;
int step = steps[idx] % modulus;
long delta = (long)step * advanceCount;
int adjusted = (int)((remainder + delta) % modulus);
if (adjusted < 0)
{
adjusted += modulus;
}
results[idx] = (byte)adjusted;
}
private static int GetModulus(int index)
{
return index switch
{
0 => 10,
1 => 8,
2 => 5,
3 => 3,
4 => 7,
5 => 11,
_ => 0,
};
}
}
@@ -1,121 +0,0 @@
using System;
using System.Buffers;
using System.Runtime.InteropServices;
using PerfectNumbers.Core;
using ILGPU;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
{
private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator;
private readonly Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel;
private MemoryBuffer1D<ulong, Stride1D.Dense> _exponentBuffer;
private MemoryBuffer1D<ulong, Stride1D.Dense> _resultBuffer;
private ulong[] _hostBuffer;
private int _capacity;
private bool _disposed;
internal MersenneNumberDivisorResidueGpuEvaluator(int batchSize)
{
_lease = GpuContextPool.RentPreferred(preferCpu: false);
_accelerator = _lease.Accelerator;
_kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeResiduesKernel);
_capacity = Math.Max(1, batchSize);
_exponentBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
}
public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_exponentBuffer.Dispose();
_resultBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
_lease.Dispose();
}
internal void ComputeResidues(ReadOnlySpan<ulong> exponents, in MontgomeryDivisorData divisorData, Span<ulong> destination)
{
if (_disposed)
{
throw new ObjectDisposedException(nameof(MersenneNumberDivisorResidueGpuEvaluator));
}
int length = exponents.Length;
if (length == 0)
{
return;
}
EnsureCapacity(length);
Span<ulong> hostSpan = _hostBuffer.AsSpan(0, length);
exponents.CopyTo(hostSpan);
ArrayView1D<ulong, Stride1D.Dense> exponentView = _exponentBuffer.View.SubView(0, length);
exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length);
ArrayView1D<ulong, Stride1D.Dense> resultView = _resultBuffer.View.SubView(0, length);
_kernel(length, divisorData, exponentView, resultView);
resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostSpan), length);
hostSpan.CopyTo(destination);
}
private void EnsureCapacity(int required)
{
if (required <= _capacity)
{
return;
}
int newCapacity = _capacity;
if (newCapacity == 0)
{
newCapacity = 1;
}
while (newCapacity < required)
{
newCapacity <<= 1;
}
_exponentBuffer.Dispose();
_resultBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
_capacity = newCapacity;
_exponentBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
}
private static void ComputeResiduesKernel(Index1D index, MontgomeryDivisorData divisorData, ArrayView<ulong> exponents, ArrayView<ulong> results)
{
int idx = index;
int length = (int)exponents.Length;
if (idx >= length)
{
return;
}
ulong exponent = exponents[idx];
ulong modulus = divisorData.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{
results[idx] = 0UL;
return;
}
results[idx] = exponent.Pow2ModBinaryGpu(modulus);
}
}
@@ -1,94 +0,0 @@
using System;
using System.Runtime.InteropServices;
using ILGPU.Runtime;
namespace PerfectNumbers.Core.Gpu;
internal static class MontgomeryOddPowerGpu
{
public static bool TryCompute(in MontgomeryDivisorData divisor, int oddPowerCount, Span<ulong> destination)
{
if (oddPowerCount <= 0)
{
return true;
}
Span<MontgomeryDivisorData> divisorSpan = stackalloc MontgomeryDivisorData[1];
divisorSpan[0] = divisor;
Span<int> countSpan = stackalloc int[1];
countSpan[0] = oddPowerCount;
Span<ulong> resultSpan = stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount];
if (!TryComputeBatch(divisorSpan, countSpan, resultSpan, PerfectNumberConstants.MaxOddPowersCount))
{
return false;
}
resultSpan[..oddPowerCount].CopyTo(destination);
return true;
}
public static bool TryComputeBatch(
ReadOnlySpan<MontgomeryDivisorData> divisors,
ReadOnlySpan<int> oddPowerCounts,
Span<ulong> results,
int stride)
{
int divisorCount = divisors.Length;
if (divisorCount == 0)
{
return true;
}
if (oddPowerCounts.Length < divisorCount)
{
throw new ArgumentException(
"oddPowerCounts span must contain at least as many elements as divisors span.",
nameof(oddPowerCounts));
}
if (results.Length < divisorCount * stride)
{
throw new ArgumentException(
"The result span is shorter than the required batched odd power storage.",
nameof(results));
}
if (GpuContextPool.ForceCpu)
{
return false;
}
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
var accelerator = lease.Accelerator;
var stream = lease.Stream;
using var divisorBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(divisorCount);
using var countBuffer = accelerator.Allocate1D<int>(divisorCount);
using var resultBuffer = accelerator.Allocate1D<ulong>(divisorCount * stride);
try
{
divisorBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(divisors), divisorCount);
countBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(oddPowerCounts), divisorCount);
resultBuffer.MemSetToZero();
var kernel = lease.MontgomeryOddPowerKernel;
kernel(stream, divisorCount, divisorBuffer.View, countBuffer.View, resultBuffer.View, stride);
stream.Synchronize();
resultBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), divisorCount * stride);
return true;
}
catch
{
return false;
}
finally
{
execution.Dispose();
lease.Dispose();
}
}
}
@@ -415,9 +415,9 @@ internal static partial class PrimeOrderGpuHeuristics
ulong? previousOrder,
PrimeOrderCalculator.PrimeOrderSearchConfig config,
in MontgomeryDivisorData divisorData,
out ulong order)
out PrimeOrderCalculator.PrimeOrderResult result)
{
order = 0UL;
result = default;
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
@@ -497,17 +497,18 @@ internal static partial class PrimeOrderGpuHeuristics
throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity.");
}
ulong order = 0UL;
resultBuffer.View.CopyToCPU(ref order, 1);
if (kernelStatus == PrimeOrderKernelStatus.FactoringFailure)
{
order = 0UL;
}
PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found
? PrimeOrderCalculator.PrimeOrderStatus.Found
: PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved;
result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order);
DisposeResources();
lease.Dispose();
return order != 0UL;
return true;
void DisposeResources()
{
@@ -811,7 +812,7 @@ internal static partial class PrimeOrderGpuHeuristics
}
ulong reduced = phi / factor;
if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL)
if (reduced.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL)
{
return false;
}
@@ -1595,7 +1596,7 @@ internal static partial class PrimeOrderGpuHeuristics
private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
{
return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL;
return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL;
}
private static ulong CalculateByDoublingKernel(ulong prime)
@@ -1715,10 +1716,10 @@ internal static partial class PrimeOrderGpuHeuristics
GpuUInt128[]? rentedExponents = null;
GpuUInt128[]? rentedResults = null;
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
try
{
var execution = lease.EnterExecutionScope();
Accelerator accelerator = lease.Accelerator;
AcceleratorStream stream = lease.Stream;
var kernel = GetPow2ModWideKernel(accelerator);
@@ -1755,6 +1756,7 @@ internal static partial class PrimeOrderGpuHeuristics
exponentBuffer.Dispose();
remainderBuffer.Dispose();
execution.Dispose();
return true;
}
catch (Exception)
@@ -1774,7 +1776,6 @@ internal static partial class PrimeOrderGpuHeuristics
ArrayPool<GpuUInt128>.Shared.Return(rentedResults, clearArray: false);
}
execution.Dispose();
lease.Dispose();
}
}
@@ -1795,7 +1796,7 @@ internal static partial class PrimeOrderGpuHeuristics
return 0UL;
}
return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false);
}
private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
@@ -1811,57 +1812,7 @@ internal static partial class PrimeOrderGpuHeuristics
private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
{
ulong exponent = exponents[index];
ulong modulus = divisor.Modulus;
// ulong nPrime = divisor.NPrime;
// ulong result = divisor.MontgomeryOne;
// ulong baseValue = divisor.MontgomeryTwo;
// ulong remainingExponent = exponent;
// if (remainingExponent == 0UL)
// {
// remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
// return;
// }
// while (remainingExponent != 0UL)
// {
// if ((remainingExponent & 1UL) != 0UL)
// {
// result = result.MontgomeryMultiply(baseValue, modulus, nPrime);
// }
// remainingExponent >>= 1;
// if (remainingExponent == 0UL)
// {
// break;
// }
// // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term.
// baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime);
// }
// remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
// Plain modular exponentiation without Montgomery reduction kept for reference.
ulong plainResult = 1UL;// % modulus;
ulong plainBase = 2UL;// % modulus;
while (exponent != 0UL)
{
if ((exponent & 1UL) != 0UL)
{
plainResult = (plainResult * plainBase) % modulus;
}
exponent >>= 1;
if (exponent == 0UL)
{
break;
}
plainBase = (plainBase * plainBase) % modulus;
}
remainders[index] = plainResult;
remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
}
private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)
@@ -12,7 +12,7 @@ public interface IMersenneNumberDivisorByDivisorTester
ulong GetAllowedMaxDivisor(ulong prime);
bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null);
bool IsPrime(ulong prime, out bool divisorsExhausted);
IDivisorScanSession CreateDivisorSession();
+94 -258
View File
@@ -1,4 +1,3 @@
using System;
using System.Buffers;
using System.Buffers.Binary;
using System.Collections.Generic;
@@ -19,12 +18,6 @@ public class MersenneDivisorCycles
// Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown
private ulong[]? _smallCycles;
[ThreadStatic]
private static Dictionary<ulong, int>? s_factorCountsPool;
[ThreadStatic]
private static Dictionary<ulong, int>? s_factorScratchPool;
public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles();
public MersenneDivisorCycles()
@@ -109,7 +102,7 @@ public class MersenneDivisorCycles
return cycleLength == exponent;
}
return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
}
public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
@@ -133,7 +126,7 @@ public class MersenneDivisorCycles
}
}
return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
}
private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
@@ -241,14 +234,14 @@ public class MersenneDivisorCycles
return (UInt128)cycle;
}
UInt128 wideOrder = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (wideOrder != UInt128.Zero)
if (wideResult.Order != UInt128.Zero)
{
return wideOrder;
return wideResult.Order;
}
// Otherwise, find order of 2 mod divisor
@@ -280,85 +273,62 @@ public class MersenneDivisorCycles
{
cycleLength = 0UL;
// The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
// so this power-of-two fast-path never triggers on that call path.
// The historical fast-path is left here for reference but commented out because it is unreachable
// on the EvenPerfectBitScanner pipeline.
// if ((divisor & (divisor - 1UL)) == 0UL)
// {
// cycleLength = 1UL;
// return true;
// }
if ((divisor & (divisor - 1UL)) == 0UL)
{
cycleLength = 1UL;
return true;
}
// EvenPerfectBitScanner filters exponents <= 3 before invoking this helper and each
// divisor candidate is at least 2 * exponent + 1, so this guard never matches.
// The fallback calculation is therefore unreachable on the active call path.
// if (divisor <= 3UL || exponent <= 1UL)
// {
// cycleLength = CalculateCycleLength(divisor, divisorData);
// return true;
// }
if (divisor <= 3UL || exponent <= 1UL)
{
cycleLength = CalculateCycleLength(divisor, divisorData);
return true;
}
ulong computedOrder = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (computedOrder != 0UL)
if (orderResult.Order != 0UL)
{
cycleLength = computedOrder;
cycleLength = orderResult.Order;
return true;
}
ulong phi = divisor - 1UL;
// Divisors reach this helper only when they exceed the cached table, which means
// divisor >= 2 * exponent + 1 and phi cannot underflow to zero here.
// The defensive underflow guard is kept for posterity but disabled for the EvenPerfectBitScanner path.
// if (phi == 0UL)
// {
// return false;
// }
if (phi == 0UL)
{
return false;
}
int twoCount = BitOperations.TrailingZeroCount(phi);
ulong reducedPhi = phi >> twoCount;
// Every divisor candidate generated by the by-divisor scan satisfies q ≡ 1 (mod 2 * exponent),
// ensuring reducedPhi keeps a factor of exponent and remains non-zero in this context.
// The divisibility check is therefore unreachable for the current caller and is commented out.
// if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
// {
// return false;
// }
if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
{
return false;
}
ulong k = reducedPhi / exponent;
// reducedPhi always contains at least one factor of exponent, so dividing by exponent leaves k >= 1.
// if (k == 0UL)
// {
// return false;
// }
if (k == 0UL)
{
return false;
}
Dictionary<ulong, int> factorCounts = AcquireFactorCountDictionary(8);
try
Dictionary<ulong, int> factorCounts = new(8);
if (twoCount > 0)
{
// The divisor candidates are odd, so phi is even and this branch always records at least one factor of two.
// if (twoCount > 0)
// {
// factorCounts[2UL] = twoCount;
// }
factorCounts[2UL] = twoCount;
if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
{
return false;
}
cycleLength = ReduceOrder(divisorData, phi, factorCounts, exponent);
return true;
}
finally
if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
{
ReleaseFactorCountDictionary(factorCounts);
return false;
}
cycleLength = ReduceOrder(divisorData, phi, factorCounts);
return true;
}
private static bool AccumulateFactors(
@@ -372,119 +342,61 @@ public class MersenneDivisorCycles
return true;
}
FactorCacheEntry cachedFactorization;
PrimeFactorPartialResult? lease = null;
bool fromCache;
try
if (!TryGetFactorization(value, cache, cacheResult, out FactorCacheEntry factorization))
{
if (!TryGetFactorization(value, cache, cacheResult, out cachedFactorization, out lease, out fromCache))
{
return false;
}
if (fromCache)
{
ulong[] primes = cachedFactorization.Primes;
byte[] exponents = cachedFactorization.Exponents;
int length = cachedFactorization.Count;
for (int i = 0; i < length; i++)
{
AddFactor(counts, primes[i], exponents[i]);
}
return true;
}
if (lease is null || lease.Count == 0)
{
return true;
}
ReadOnlySpan<ulong> primesSpan = lease.Primes;
ReadOnlySpan<byte> exponentsSpan = lease.Exponents;
int spanLength = lease.Count;
for (int i = 0; i < spanLength; i++)
{
AddFactor(counts, primesSpan[i], exponentsSpan[i]);
}
return true;
return false;
}
finally
ulong[] primes = factorization.Primes;
byte[] exponents = factorization.Exponents;
int length = factorization.Count;
for (int i = 0; i < length; i++)
{
lease?.Dispose();
AddFactor(counts, primes[i], exponents[i]);
}
return true;
}
private static bool TryGetFactorization(
ulong value,
Dictionary<ulong, FactorCacheEntry>? cache,
bool cacheResult,
out FactorCacheEntry factorization,
out PrimeFactorPartialResult? lease,
out bool fromCache)
out FactorCacheEntry factorization)
{
factorization = default;
lease = null;
fromCache = false;
if (value <= 1UL)
{
return true;
}
if (cache is not null && cache.TryGetValue(value, out factorization))
{
fromCache = true;
return true;
}
Dictionary<ulong, int> scratch = AcquireFactorScratchDictionary(8);
try
Dictionary<ulong, int> scratch = new(8);
if (!TryFactorIntoCountsInternal(value, scratch))
{
if (!TryFactorIntoCountsInternal(value, scratch))
{
factorization = default;
return false;
}
int count = scratch.Count;
if (count == 0)
{
return true;
}
PrimeFactorPartialResult result = PrimeFactorPartialResult.Rent(count);
lease = result;
Span<ulong> primes = result.PrimeWriteSpan;
Span<byte> exponents = result.ExponentWriteSpan;
int index = 0;
foreach (KeyValuePair<ulong, int> entry in scratch)
{
primes[index] = entry.Key;
exponents[index] = checked((byte)entry.Value);
index++;
}
result.CommitCount(index);
result.Sort();
if (cache is not null && cacheResult)
{
FactorCacheEntry cacheEntry = result.ToCacheEntry();
cache[value] = cacheEntry;
}
return true;
factorization = default!;
return false;
}
finally
int count = scratch.Count;
ulong[] primes = new ulong[count];
byte[] exponents = new byte[count];
int index = 0;
foreach (KeyValuePair<ulong, int> entry in scratch)
{
ReleaseFactorScratchDictionary(scratch);
primes[index] = entry.Key;
exponents[index] = checked((byte)entry.Value);
index++;
}
Array.Sort(primes, exponents);
factorization = new FactorCacheEntry(primes, exponents, count);
if (cache is not null && cacheResult)
{
cache[value] = factorization;
}
return true;
}
private static bool TryFactorIntoCountsInternal(ulong value, Dictionary<ulong, int> counts)
@@ -559,22 +471,15 @@ public class MersenneDivisorCycles
}
}
private static ulong ReduceOrder(
in MontgomeryDivisorData divisorData,
ulong initialOrder,
Dictionary<ulong, int> factorCounts,
ulong targetOrder)
private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
{
if (factorCounts.Count == 0 || initialOrder == targetOrder)
if (factorCounts.Count == 0)
{
return initialOrder;
}
int count = factorCounts.Count;
bool rentFromPool = count >= PerfectNumberConstants.PooledArrayThreshold;
ulong[] primes = rentFromPool
? ThreadLocalArrayPool<ulong>.Shared.Rent(count)
: new ulong[count];
ulong[] primes = ArrayPool<ulong>.Shared.Rent(count);
try
{
int index = 0;
@@ -600,16 +505,9 @@ public class MersenneDivisorCycles
}
ulong candidate = order / prime;
if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL)
if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
{
order = candidate;
if (order == targetOrder)
{
// Early-exit once the order matches the tested exponent;
// further factor processing would not change the outcome for the by-divisor scan.
return order;
}
continue;
}
@@ -621,10 +519,7 @@ public class MersenneDivisorCycles
}
finally
{
if (rentFromPool)
{
ThreadLocalArrayPool<ulong>.Shared.Return(primes, clearArray: false);
}
ArrayPool<ulong>.Shared.Return(primes, clearArray: false);
}
}
@@ -697,53 +592,6 @@ public class MersenneDivisorCycles
}
}
// Thread-local dictionary pools keep the factoring helpers reusable so the GC stays out of the by-divisor CPU loop.
private static Dictionary<ulong, int> AcquireFactorCountDictionary(int capacityHint)
{
Dictionary<ulong, int>? dictionary = s_factorCountsPool;
if (dictionary is null)
{
return new Dictionary<ulong, int>(capacityHint);
}
s_factorCountsPool = null;
dictionary.Clear();
dictionary.EnsureCapacity(capacityHint);
return dictionary;
}
private static void ReleaseFactorCountDictionary(Dictionary<ulong, int> dictionary)
{
dictionary.Clear();
if (s_factorCountsPool is null)
{
s_factorCountsPool = dictionary;
}
}
private static Dictionary<ulong, int> AcquireFactorScratchDictionary(int capacityHint)
{
Dictionary<ulong, int>? dictionary = s_factorScratchPool;
if (dictionary is null)
{
return new Dictionary<ulong, int>(capacityHint);
}
s_factorScratchPool = null;
dictionary.Clear();
dictionary.EnsureCapacity(capacityHint);
return dictionary;
}
private static void ReleaseFactorScratchDictionary(Dictionary<ulong, int> dictionary)
{
dictionary.Clear();
if (s_factorScratchPool is null)
{
s_factorScratchPool = dictionary;
}
}
private static ulong NextRandomUInt64()
{
Span<byte> buffer = stackalloc byte[8];
@@ -810,11 +658,7 @@ public class MersenneDivisorCycles
tasks[taskIndex] = Task.Run(() =>
{
ulong rangeLength = localEnd - localStart + 1UL;
int localLength = checked((int)rangeLength);
bool rentLocalCycles = localLength >= PerfectNumberConstants.PooledArrayThreshold;
(ulong divisor, ulong cycle)[] localCycles = rentLocalCycles
? ThreadLocalArrayPool<(ulong, ulong)>.Shared.Rent(localLength)
: new (ulong divisor, ulong cycle)[localLength];
(ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength));
int localCycleIndex = 0;
try
@@ -862,10 +706,7 @@ public class MersenneDivisorCycles
}
finally
{
if (rentLocalCycles)
{
ThreadLocalArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false);
}
ArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false);
}
});
}
@@ -880,7 +721,7 @@ public class MersenneDivisorCycles
ulong start = 3UL;
var locker = new object();
var pool = ThreadLocalArrayPool<ulong>.Shared;
var pool = ArrayPool<ulong>.Shared;
ulong batchSizeUL = (ulong)batchSize, d, end;
int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx;
ulong[] divisors, outCycles, validDivisors;
@@ -1085,34 +926,29 @@ public class MersenneDivisorCycles
internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength)
{
// The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
// so this power-of-two fast-path never triggers on that call path.
// The early exit is retained for reference but commented out because it is unreachable on the scanner path.
// if ((divisor & (divisor - 1UL)) == 0UL)
// {
// cycleLength = 1UL;
// return true;
// }
if ((divisor & (divisor - 1UL)) == 0UL)
{
cycleLength = 1UL;
return true;
}
// By the divisor-by-divisor pipeline, candidates have the shape 2 * exponent * k + 1 with exponent >= 5,
// so the divisor never drops below 11 on this call path.
// if (divisor <= 3UL)
// {
// cycleLength = CalculateCycleLengthFallback(divisor);
// return true;
// }
if (divisor <= 3UL)
{
cycleLength = CalculateCycleLengthFallback(divisor);
return true;
}
if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
{
ulong computedOrder = PrimeOrderCalculator.Calculate(
PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
if (computedOrder != 0UL)
if (orderResult.Order != 0UL)
{
cycleLength = computedOrder;
cycleLength = orderResult.Order;
return true;
}
}
@@ -18,8 +18,7 @@ public static class MersenneNumberDivisorByDivisorTester
Action markComposite,
Action clearComposite,
Action<ulong, bool, bool, bool> printResult,
int threadCount,
TimeSpan? primeTestLimit)
int threadCount)
{
if (candidates.Count == 0)
{
@@ -113,12 +112,9 @@ public static class MersenneNumberDivisorByDivisorTester
ulong maxPrime = 0UL;
int primeWriteIndex = 0;
// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
// We don't need to additionally check it.
// using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
// bool hasPrime = primeEnumerator.MoveNext();
// ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
bool hasPrime = primeEnumerator.MoveNext();
ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
for (int index = startIndex; index < candidateCount; index++)
{
@@ -131,29 +127,25 @@ public static class MersenneNumberDivisorByDivisorTester
continue;
}
// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
// We don't need to additionally check it.
// while (hasPrime && currentPrime < candidate)
// {
// hasPrime = primeEnumerator.MoveNext();
// if (!hasPrime)
// {
// break;
// }
while (hasPrime && currentPrime < candidate)
{
hasPrime = primeEnumerator.MoveNext();
if (!hasPrime)
{
break;
}
// currentPrime = primeEnumerator.Current;
// }
currentPrime = primeEnumerator.Current;
}
// if (!hasPrime || currentPrime != candidate)
// if (!Number.IsPrime(candidate))
// {
// markComposite();
// printResult(candidate, false, false, false);
// continue;
// }
if (!hasPrime || currentPrime != candidate)
{
markComposite();
printResult(candidate, false, false, false);
continue;
}
primesSpan[primeWriteIndex++] = candidate;
primesSpan[primeWriteIndex++] = candidate;
if (candidate > maxPrime)
{
@@ -244,7 +236,7 @@ public static class MersenneNumberDivisorByDivisorTester
void ProcessPrime(ulong prime)
{
bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted, primeTestLimit);
bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted);
if (!isPrime)
{
@@ -257,52 +249,21 @@ public static class MersenneNumberDivisorByDivisorTester
printResult(prime, true, divisorsExhausted, true);
}
if (workerCount > filteredPrimes.Count)
{
workerCount = filteredPrimes.Count;
}
if (workerCount <= 1)
if (workerCount == 1)
{
foreach (ulong prime in filteredPrimes)
{
ProcessPrime(prime);
}
return;
}
(int Start, int End)[] workerRanges = new (int, int)[workerCount];
int totalPrimeCount = filteredPrimes.Count;
int baseGroupSize = totalPrimeCount / workerCount;
int remaining = totalPrimeCount % workerCount;
int groupStartIndex = 0;
for (int workerIndex = 0; workerIndex < workerCount; workerIndex++)
else
{
int groupSize = baseGroupSize;
if (workerIndex < remaining)
ParallelOptions options = new()
{
groupSize++;
}
MaxDegreeOfParallelism = workerCount
};
startIndex = groupStartIndex;
int endIndex = startIndex + groupSize;
workerRanges[workerIndex] = (startIndex, endIndex);
groupStartIndex = endIndex;
Parallel.ForEach(filteredPrimes, options, ProcessPrime);
}
ParallelOptions parallelOptions = new()
{
MaxDegreeOfParallelism = workerCount
};
Parallel.ForEach(workerRanges, parallelOptions, range =>
{
for (int candidateIndex = range.Start; candidateIndex < range.End; candidateIndex++)
{
ProcessPrime(filteredPrimes[candidateIndex]);
}
});
}
}
@@ -371,6 +371,10 @@ public static class MersennePrimeFactorTester
return 2UL;
}
var rng = new Random(1234567); // TODO: Replace this RNG with the deterministic span-based
// sequence measured faster in PollardRhoBenchmarks so we avoid
// per-call Random allocations.
ulong c, d, x, y, diff;
while (true)
{
@@ -379,8 +383,8 @@ public static class MersennePrimeFactorTester
return 0UL;
}
c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
c = (ulong)rng.Next(2, int.MaxValue);
x = (ulong)rng.Next(2, int.MaxValue);
y = x;
d = 1UL;
@@ -519,6 +523,9 @@ public static class MersennePrimeFactorTester
return 2;
}
var rng = new Random(1234567); // TODO: Replace with the deterministic UInt128 sequence from
// PollardRho128Benchmarks so we eliminate Random allocations.
UInt128 c, d, x, y, diff;
while (true)
{
@@ -527,8 +534,8 @@ public static class MersennePrimeFactorTester
return 0;
}
c = (DeterministicRandom.NextUInt128() % (n - 1)) + 1;
x = (DeterministicRandom.NextUInt128() % (n - 2)) + 2;
c = (UInt128)rng.NextInt64(2, long.MaxValue);
x = (UInt128)rng.NextInt64(2, long.MaxValue);
y = x;
d = 1;
@@ -45,4 +45,32 @@ public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong
montgomeryTwo,
montgomeryTwoSquared);
}
}
internal static class MontgomeryDivisorDataCache
{
[ThreadStatic]
private static bool s_hasCachedValue;
[ThreadStatic]
private static ulong s_cachedModulus;
[ThreadStatic]
private static MontgomeryDivisorData s_cachedData;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static MontgomeryDivisorData Get(ulong modulus)
{
if (s_hasCachedValue && s_cachedModulus == modulus)
{
return s_cachedData;
}
MontgomeryDivisorData computed = MontgomeryDivisorData.FromModulus(modulus);
s_cachedModulus = modulus;
s_cachedData = computed;
s_hasCachedValue = true;
return computed;
}
}
-243
View File
@@ -1,243 +0,0 @@
using System;
using System.Numerics;
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
/// <summary>
/// Mutable helper for 128-bit arithmetic derived from the .NET runtime UInt128 implementation.
/// Portions adapted from https://github.com/dotnet/runtime/blob/v8.0.0/src/libraries/System.Private.CoreLib/src/System/UInt128.cs.
/// </summary>
public struct MutableUInt128
{
public ulong High;
public ulong Low;
public MutableUInt128(ulong value)
{
Low = value;
}
public MutableUInt128(ulong high, ulong low)
{
High = high;
Low = low;
}
public MutableUInt128(UInt128 value)
{
High = (ulong)(value >> 64);
Low = (ulong)value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Set(ulong value)
{
High = 0UL;
Low = value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Set(UInt128 value)
{
High = (ulong)(value >> 64);
Low = (ulong)value;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Set(ulong high, ulong low)
{
High = high;
Low = low;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public readonly UInt128 ToUInt128() => ((UInt128)High << 64) | Low;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Add(ulong value)
{
ulong low = Low + value;
if (low < value)
{
High++;
}
Low = low;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Add(MutableUInt128 value)
{
ulong temp = Low + value.Low;
Low = temp;
temp = temp < Low ? 1UL : 0UL;
temp += High;
temp += value.High;
High = temp;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Sub(ulong value)
{
ulong temp = Low;
ulong result = temp - value;
temp = result > temp ? 1UL : 0UL;
Low = result;
High -= temp;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Sub(MutableUInt128 value)
{
ulong temp = Low;
ulong result = temp - value.Low;
temp = temp < value.Low ? 1UL : 0UL;
Low = result;
High = High - value.High - temp;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Mul(ulong value)
{
ulong currentHigh = High;
ulong currentLow = Low;
ulong temp = Mul64(currentLow, value, out currentLow);
Mul64(currentHigh, value, out currentHigh);
temp += currentHigh;
Low = currentLow;
High = temp;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Mul64(ulong left, ulong right, out ulong low)
{
const ulong Mask32 = 0xFFFF_FFFFUL;
ulong leftHigh = left >> 32;
ulong rightLow = right & Mask32;
ulong leftLow = left & Mask32;
ulong lowLow = leftLow * rightLow;
ulong temp = right >> 32;
leftLow *= temp;
rightLow *= leftHigh;
leftHigh *= temp;
temp = (lowLow >> 32) + (leftLow & Mask32) + (rightLow & Mask32);
leftHigh += (leftLow >> 32) + (rightLow >> 32) + (temp >> 32);
low = ((temp & Mask32) << 32) | (lowLow & Mask32);
return leftHigh;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void MulAdd(ulong multiplier, ulong addend)
{
const ulong Mask32 = 0xFFFF_FFFFUL;
ulong multiplierLow = multiplier & Mask32;
ulong multiplierHigh = multiplier >> 32;
ulong product11 = Low;
ulong product01 = product11 & Mask32;
product11 >>= 32;
ulong product00 = product01 * multiplierLow;
product01 *= multiplierHigh;
ulong product10 = product11 * multiplierLow;
product11 *= multiplierHigh;
ulong high = (product00 >> 32) + (product01 & Mask32) + (product10 & Mask32);
product11 = product11 + (product01 >> 32) + (product10 >> 32) + (high >> 32);
product10 = ((high & Mask32) << 32) | (product00 & Mask32);
high = High;
product00 = high & Mask32;
high >>= 32;
product01 = product00 * multiplierLow;
product00 *= multiplierHigh;
high *= multiplierLow;
high = (product01 >> 32) + (product00 & Mask32) + (high & Mask32);
product01 = ((high & Mask32) << 32) | (product01 & Mask32);
// Calculate low. We're assigning addend to local variable for a tiny bit better performance.
product00 = addend;
high = product10 + product00;
Low = high;
// Calculate high including carry
high = high < product00 ? 1UL : 0UL;
high += product11 + product01;
// Discard overflow above 128 bits to match Multiply's truncation semantics.
High = high;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public readonly ulong Mod(ulong modulus)
{
// This case will never happen in production code.
// if (modulus == 0UL)
// {
// throw new DivideByZeroException();
// }
ulong high = High;
if (high == 0UL)
{
return Low % modulus;
}
MutableUInt128 remainder = default;
int highBitLength = 64 - BitOperations.LeadingZeroCount(high);
if (highBitLength > 0)
{
ProcessWord(ref remainder, high, highBitLength, modulus);
}
ProcessWord(ref remainder, Low, 64, modulus);
return remainder.Low;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void ProcessWord(ref MutableUInt128 remainder, ulong word, int bitCount, ulong modulus)
{
for (int bit = bitCount - 1; bit >= 0; bit--)
{
remainder.ShiftLeft();
if (((word >> bit) & 1UL) != 0)
{
remainder.Add(1UL);
}
remainder.Reduce(modulus);
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void ShiftLeft()
{
High = (High << 1) | (Low >> 63);
Low <<= 1;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void Reduce(ulong modulus)
{
if (High != 0UL || Low >= modulus)
{
Sub(modulus);
if (High != 0UL || Low >= modulus)
{
Sub(modulus);
}
}
}
}
@@ -1,94 +0,0 @@
using System;
namespace PerfectNumbers.Core
{
public sealed class PartialFactorResult
{
[ThreadStatic]
private static PartialFactorResult? s_poolHead;
public static readonly PartialFactorResult Empty = new()
{
Cofactor = 1UL,
FullyFactored = true,
};
public FactorEntry[]? Factors;
public ulong Cofactor;
public bool FullyFactored;
public int Count;
private PartialFactorResult? _next;
private bool _factorsFromPool;
private PartialFactorResult()
{
}
public static PartialFactorResult Rent(in FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count)
{
PartialFactorResult? instance = s_poolHead;
if (instance is null)
{
instance = new PartialFactorResult();
}
else
{
s_poolHead = instance._next;
}
FactorEntry[]? factorArray = factors;
instance._next = null;
instance.Factors = factorArray;
instance.Cofactor = cofactor;
instance.FullyFactored = fullyFactored;
instance.Count = count;
instance._factorsFromPool = factorArray is not null && factorArray.Length >= PerfectNumberConstants.PooledArrayThreshold;
return instance;
}
public PartialFactorResult WithAdditionalPrime(ulong prime)
{
FactorEntry[]? source = Factors;
if (source is null || Count == 0)
{
FactorEntry[] single = new FactorEntry[1];
single[0] = new FactorEntry(prime);
return Rent(single, 1UL, true, 1);
}
int newCount = Count + 1;
FactorEntry[] extended = newCount >= PerfectNumberConstants.PooledArrayThreshold
? ThreadLocalArrayPool<FactorEntry>.Shared.Rent(newCount)
: new FactorEntry[newCount];
Array.Copy(source, 0, extended, 0, Count);
extended[Count] = new FactorEntry(prime);
Span<FactorEntry> span = extended.AsSpan(0, newCount);
span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
return Rent(extended, 1UL, true, newCount);
}
public void Dispose()
{
if (this == Empty)
{
return;
}
FactorEntry[]? factors = Factors;
if (factors is not null)
{
Factors = null;
if (_factorsFromPool)
{
ThreadLocalArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
}
}
_factorsFromPool = false;
_next = s_poolHead;
s_poolHead = this;
}
}
}
@@ -8,7 +8,6 @@ public static class PerfectNumberConstants
public const int MaxQForDivisorCycles = 4_000_000;
public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000;
public const int MaxOddPowersCount = 128;
public const int PooledArrayThreshold = 64;
// TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal
// divisor-cycle datasets we generate offline.
// TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune
@@ -22,7 +22,4 @@
<InternalsVisibleTo Include="EvenPerfectBitScanner.Benchmarks" />
</ItemGroup>
<ItemGroup>
<Compile Remove=".temp.cs" />
</ItemGroup>
</Project>
@@ -1,150 +0,0 @@
using System;
using System.Buffers;
namespace PerfectNumbers.Core
{
internal sealed class PrimeFactorPartialResult : IDisposable
{
[ThreadStatic]
private static PrimeFactorPartialResult? s_poolHead;
private PrimeFactorPartialResult? _next;
private ulong[]? _primes;
private byte[]? _exponents;
private bool _arraysFromPool;
private int _capacity;
private PrimeFactorPartialResult()
{
}
public static PrimeFactorPartialResult Rent(int capacity)
{
PrimeFactorPartialResult? instance = s_poolHead;
if (instance is null)
{
instance = new PrimeFactorPartialResult();
}
else
{
s_poolHead = instance._next;
}
instance._next = null;
instance.EnsureCapacity(capacity);
instance.Count = 0;
return instance;
}
public int Count { get; private set; }
internal ReadOnlySpan<ulong> Primes => _primes is null ? ReadOnlySpan<ulong>.Empty : _primes.AsSpan(0, Count);
internal ReadOnlySpan<byte> Exponents => _exponents is null ? ReadOnlySpan<byte>.Empty : _exponents.AsSpan(0, Count);
internal Span<ulong> PrimeWriteSpan => _primes is null ? Span<ulong>.Empty : _primes.AsSpan(0, _capacity);
internal Span<byte> ExponentWriteSpan => _exponents is null ? Span<byte>.Empty : _exponents.AsSpan(0, _capacity);
internal void CommitCount(int count)
{
Count = count;
}
internal void Sort()
{
if (Count == 0 || _primes is null || _exponents is null)
{
return;
}
Array.Sort(_primes, _exponents, 0, Count);
}
internal MersenneDivisorCycles.FactorCacheEntry ToCacheEntry()
{
if (Count == 0 || _primes is null || _exponents is null)
{
return default;
}
ulong[] primes = new ulong[Count];
byte[] exponents = new byte[Count];
Array.Copy(_primes, 0, primes, 0, Count);
Array.Copy(_exponents, 0, exponents, 0, Count);
return new MersenneDivisorCycles.FactorCacheEntry(primes, exponents, Count);
}
public void Dispose()
{
if (_arraysFromPool)
{
if (_primes is not null)
{
ThreadLocalArrayPool<ulong>.Shared.Return(_primes, clearArray: false);
}
if (_exponents is not null)
{
ThreadLocalArrayPool<byte>.Shared.Return(_exponents, clearArray: false);
}
_arraysFromPool = false;
_primes = null;
_exponents = null;
}
_capacity = 0;
Count = 0;
_next = s_poolHead;
s_poolHead = this;
}
private void EnsureCapacity(int capacity)
{
if (capacity <= 0)
{
_capacity = 0;
return;
}
if (_primes is not null && _primes.Length >= capacity)
{
_capacity = capacity;
return;
}
if (_arraysFromPool)
{
if (_primes is not null)
{
ThreadLocalArrayPool<ulong>.Shared.Return(_primes, clearArray: false);
}
if (_exponents is not null)
{
ThreadLocalArrayPool<byte>.Shared.Return(_exponents, clearArray: false);
}
_arraysFromPool = false;
_primes = null;
_exponents = null;
}
if (capacity >= PerfectNumberConstants.PooledArrayThreshold)
{
_primes = ThreadLocalArrayPool<ulong>.Shared.Rent(capacity);
_exponents = ThreadLocalArrayPool<byte>.Shared.Rent(capacity);
_arraysFromPool = true;
}
else
{
_primes = new ulong[capacity];
_exponents = new byte[capacity];
}
_capacity = capacity;
}
}
}
@@ -9,15 +9,30 @@ namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
public static UInt128 Calculate(
internal readonly struct PrimeOrderResultWide
{
public PrimeOrderResultWide(PrimeOrderStatus status, UInt128 order)
{
Status = status;
Order = order;
}
public PrimeOrderStatus Status { get; }
public UInt128 Order { get; }
}
public static PrimeOrderResultWide Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
=> Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
public static PrimeOrderResultWide Calculate(
UInt128 prime,
UInt128? previousOrder,
PrimeOrderSearchConfig config,
PrimeOrderHeuristicDevice device)
{
var scope = UsePow2Mode(device);
using var scope = UsePow2Mode(device);
MontgomeryDivisorData divisorData;
UInt128 result;
if (prime <= ulong.MaxValue)
{
ulong? previous = null;
@@ -36,49 +51,41 @@ internal static partial class PrimeOrderCalculator
ulong prime64 = (ulong)prime;
divisorData = MontgomeryDivisorData.FromModulus(prime64);
ulong order64 = Calculate(prime64, previous, divisorData, config, device);
result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device);
return new PrimeOrderResultWide(result.Status, result.Order);
}
else
{
divisorData = default;
result = CalculateWideInternal(prime, previousOrder, divisorData, config);
}
scope.Dispose();
return result;
return CalculateWideInternal(prime, previousOrder, divisorData, config);
}
private static UInt128 CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
{
if (prime <= UInt128.One)
{
return UInt128.One;
return new PrimeOrderResultWide(PrimeOrderStatus.Found, UInt128.One);
}
if (prime == (UInt128)3UL)
{
return (UInt128)2UL;
return new PrimeOrderResultWide(PrimeOrderStatus.Found, (UInt128)2UL);
}
UInt128 phi = prime - UInt128.One;
var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
UInt128 result;
if (phiFactors.Factors is null)
{
result = FinishStrictlyWide(prime, divisorData);
}
else
{
result = RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
return FinishStrictlyWide(prime, divisorData, config.Mode);
}
debugScope.Dispose();
return result;
return RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
}
private static UInt128 RunHeuristicPipelineWide(
private static PrimeOrderResultWide RunHeuristicPipelineWide(
UInt128 prime,
UInt128? previousOrder,
in MontgomeryDivisorData divisorData,
@@ -88,7 +95,7 @@ internal static partial class PrimeOrderCalculator
{
if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
{
return phi;
return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi);
}
UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData);
@@ -96,25 +103,27 @@ internal static partial class PrimeOrderCalculator
if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config))
{
return candidateOrder;
return new PrimeOrderResultWide(PrimeOrderStatus.Found, candidateOrder);
}
if (config.Mode == PrimeOrderMode.Strict)
{
return FinishStrictlyWide(prime, divisorData);
return FinishStrictlyWide(prime, divisorData, PrimeOrderMode.Strict);
}
if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
{
return order;
return new PrimeOrderResultWide(PrimeOrderStatus.Found, order);
}
return candidateOrder;
return FinishStrictlyWide(prime, divisorData, config.Mode);
}
private static UInt128 FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData)
private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData, PrimeOrderMode mode)
{
return CalculateByFactorizationWide(prime, divisorData);
UInt128 strictOrder = CalculateByFactorizationWide(prime, divisorData);
PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved;
return new PrimeOrderResultWide(status, strictOrder);
}
private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
@@ -162,7 +171,7 @@ internal static partial class PrimeOrderCalculator
return order;
}
tempArray = ThreadLocalArrayPool<FactorEntry128>.Shared.Rent(length + 1);
tempArray = ArrayPool<FactorEntry128>.Shared.Rent(length + 1);
Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
factorSpan.CopyTo(buffer);
@@ -200,7 +209,7 @@ internal static partial class PrimeOrderCalculator
{
if (tempArray is not null)
{
ThreadLocalArrayPool<FactorEntry128>.Shared.Return(tempArray, clearArray: false);
ArrayPool<FactorEntry128>.Shared.Return(tempArray, clearArray: false);
}
}
}
@@ -664,7 +673,7 @@ internal static partial class PrimeOrderCalculator
return new PartialFactorResult128(null, value, false, 0);
}
FactorEntry128[]? rented = counts.Count > 0 ? ThreadLocalArrayPool<FactorEntry128>.Shared.Rent(counts.Count) : null;
FactorEntry128[]? rented = counts.Count > 0 ? ArrayPool<FactorEntry128>.Shared.Rent(counts.Count) : null;
int index = 0;
if (rented is not null)
{
@@ -682,7 +691,7 @@ internal static partial class PrimeOrderCalculator
span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
resultArray = new FactorEntry128[index];
span.CopyTo(resultArray);
ThreadLocalArrayPool<FactorEntry128>.Shared.Return(rented, clearArray: false);
ArrayPool<FactorEntry128>.Shared.Return(rented, clearArray: false);
}
bool fullyFactored = cofactor == UInt128.One;
@@ -1062,20 +1071,15 @@ internal static partial class PrimeOrderCalculator
private readonly struct FactorEntry128
{
public FactorEntry128(UInt128 value)
{
Value = value;
}
public FactorEntry128(UInt128 value, int exponent)
{
Value = value;
Exponent = exponent;
}
public FactorEntry128(UInt128 value, int exponent)
{
Value = value;
Exponent = exponent;
}
public UInt128 Value { get; }
public readonly UInt128 Value;
public readonly int Exponent = 1;
public int Exponent { get; }
}
private readonly struct PartialFactorResult128
@@ -1088,13 +1092,13 @@ internal static partial class PrimeOrderCalculator
Count = count;
}
public readonly FactorEntry128[]? Factors;
public FactorEntry128[]? Factors { get; }
public readonly UInt128 Cofactor;
public UInt128 Cofactor { get; }
public readonly bool FullyFactored;
public bool FullyFactored { get; }
public readonly int Count;
public int Count { get; }
public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0);
@@ -1102,13 +1106,14 @@ internal static partial class PrimeOrderCalculator
{
if (Factors is null)
{
FactorEntry128[] local = [new FactorEntry128(prime)];
return new PartialFactorResult128(local, UInt128.One, true, 1);
FactorEntry128[] local = new FactorEntry128[1];
local[0] = new FactorEntry128(prime, 1);
return new PartialFactorResult128(local, UInt128.One, true, 1);
}
FactorEntry128[] extended = new FactorEntry128[Count + 1];
Array.Copy(Factors, extended, Count);
extended[Count] = new FactorEntry128(prime);
extended[Count] = new FactorEntry128(prime, 1);
return new PartialFactorResult128(extended, UInt128.One, true, Count + 1);
}
}
File diff suppressed because it is too large. Load diff
-72
View File
@@ -1,72 +0,0 @@
using System;
using System.Diagnostics;
namespace PerfectNumbers.Core;
internal readonly struct PrimeTestTimeLimit
{
private readonly long _deadlineTimestamp;
private PrimeTestTimeLimit(bool isActive, long deadlineTimestamp)
{
IsActive = isActive;
_deadlineTimestamp = deadlineTimestamp;
}
internal bool IsActive { get; }
internal static bool TryCreate(TimeSpan? timeLimit, out PrimeTestTimeLimit limit)
{
if (!timeLimit.HasValue)
{
limit = new PrimeTestTimeLimit(isActive: false, deadlineTimestamp: long.MaxValue);
return true;
}
TimeSpan limitValue = timeLimit.GetValueOrDefault();
if (limitValue <= TimeSpan.Zero)
{
limit = new PrimeTestTimeLimit(isActive: true, deadlineTimestamp: long.MinValue);
return false;
}
double seconds = limitValue.TotalSeconds;
if (double.IsInfinity(seconds))
{
limit = new PrimeTestTimeLimit(isActive: true, deadlineTimestamp: long.MaxValue);
return true;
}
double scaled = seconds * Stopwatch.Frequency;
long limitTicks = scaled >= long.MaxValue ? long.MaxValue : (long)scaled;
long startTimestamp = Stopwatch.GetTimestamp();
long deadline = startTimestamp <= long.MaxValue - limitTicks ? startTimestamp + limitTicks : long.MaxValue;
if (deadline < startTimestamp)
{
deadline = long.MaxValue;
}
limit = new PrimeTestTimeLimit(isActive: true, deadlineTimestamp: deadline);
return true;
}
internal bool HasExpired()
{
if (!IsActive)
{
return false;
}
if (_deadlineTimestamp == long.MaxValue)
{
return false;
}
if (_deadlineTimestamp == long.MinValue)
{
return true;
}
return Stopwatch.GetTimestamp() >= _deadlineTimestamp;
}
}
@@ -1,27 +0,0 @@
using System;
using System.Buffers;
using System.Runtime.CompilerServices;
namespace PerfectNumbers.Core;
internal static class ThreadLocalArrayPool<T>
{
[ThreadStatic]
private static ArrayPool<T>? s_threadPool;
public static ArrayPool<T> Shared
{
[MethodImpl(MethodImplOptions.AggressiveInlining)]
get
{
ArrayPool<T>? pool = s_threadPool;
if (pool is null)
{
pool = ArrayPool<T>.Create();
s_threadPool = pool;
}
return pool;
}
}
}
+8 -6
View File
@@ -218,7 +218,7 @@ public static class UInt128Extensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool IsBitSet(GpuUInt128 value, int bitIndex)
private static bool IsBitSet(in GpuUInt128 value, int bitIndex)
{
if (bitIndex >= 64)
{
@@ -229,16 +229,18 @@ public static class UInt128Extensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount)
private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount)
{
ulong mask = (1UL << windowBitCount) - 1UL;
if (windowStart != 0)
{
exponent.ShiftRight(windowStart);
return exponent.Low & mask;
GpuUInt128 shifted = exponent;
shifted.ShiftRight(windowStart);
ulong mask = (1UL << windowBitCount) - 1UL;
return shifted.Low & mask;
}
return exponent.Low & mask;
ulong directMask = (1UL << windowBitCount) - 1UL;
return exponent.Low & directMask;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
+8 -171
View File
@@ -1,5 +1,4 @@
using System;
using System.Buffers;
using System.Numerics;
using System.Runtime.CompilerServices;
using ILGPU.Algorithms;
@@ -77,7 +76,7 @@ public static class ULongExtensions
public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift)
{
ulong high = value >> shift;
ulong mask = (1UL << shift) - 1UL;
ulong mask = (1UL << shift) - 1UL;
ulong low = value & mask;
UInt128 highContribution = (UInt128)high * multiplier;
@@ -310,37 +309,6 @@ public static class ULongExtensions
return state.MulModWithNativeModulo(b, modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus)
{
if (modulus <= 1UL)
{
return 0UL;
}
ulong remainingExponent = exponent;
ulong result = 1UL % modulus;
ulong baseValue = 2UL % modulus;
while (remainingExponent != 0UL)
{
if ((remainingExponent & 1UL) != 0UL)
{
result = result.MulMod64GpuCompatible(baseValue, modulus);
}
remainingExponent >>= 1;
if (remainingExponent == 0UL)
{
break;
}
baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus);
}
return result;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime)
{
@@ -358,7 +326,7 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
{
ulong modulus = divisor.Modulus;
if (exponent == 0UL)
@@ -377,93 +345,7 @@ public static class ULongExtensions
ulong result = divisor.MontgomeryOne;
ulong nPrime = divisor.NPrime;
Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : default;
ulong[]? pooledOddPowers = null;
if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold)
{
oddPowerStorage = new ulong[oddPowerCount];
}
else
{
pooledOddPowers = ArrayPool<ulong>.Shared.Rent(oddPowerCount);
oddPowerStorage = pooledOddPowers;
}
Span<ulong> oddPowers = oddPowerStorage[..oddPowerCount];
InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers);
int index = bitLength - 1;
while (index >= 0)
{
if (((exponent >> index) & 1UL) == 0UL)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
index--;
continue;
}
int windowStart = index - windowSize + 1;
if (windowStart < 0)
{
windowStart = 0;
}
while (((exponent >> windowStart) & 1UL) == 0UL)
{
windowStart++;
}
int windowLength = index - windowStart + 1;
for (int square = 0; square < windowLength; square++)
{
result = result.MontgomeryMultiply(result, modulus, nPrime);
}
ulong mask = (1UL << windowLength) - 1UL;
ulong windowValue = (exponent >> windowStart) & mask;
int tableIndex = (int)((windowValue - 1UL) >> 1);
ulong multiplier = oddPowers[tableIndex];
result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
index = windowStart - 1;
}
if (pooledOddPowers is not null)
{
ArrayPool<ulong>.Shared.Return(pooledOddPowers);
}
if (keepMontgomery)
{
return result;
}
return result.MontgomeryMultiply(1UL, modulus, nPrime);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
{
ulong modulus = divisor.Modulus;
if (exponent == 0UL)
{
return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
}
if (exponent <= Pow2WindowFallbackThreshold)
{
return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
}
int bitLength = GetPortableBitLength(exponent);
int windowSize = GetWindowSize(bitLength);
int oddPowerCount = 1 << (windowSize - 1);
ulong result = divisor.MontgomeryOne;
ulong nPrime = divisor.NPrime;
ulong[] oddPowers = InitializeMontgomeryOddPowersGpu(divisor, modulus, nPrime, oddPowerCount);
ulong[] oddPowers = InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowerCount);
int index = bitLength - 1;
while (index >= 0)
@@ -579,38 +461,7 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> destination)
{
int oddPowerCount = destination.Length;
if (oddPowerCount == 0)
{
return;
}
bool computedOnGpu = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
&& MontgomeryOddPowerGpu.TryCompute(divisor, oddPowerCount, destination);
if (computedOnGpu)
{
return;
}
destination[0] = divisor.MontgomeryTwo;
if (oddPowerCount == 1)
{
return;
}
ulong square = divisor.MontgomeryTwoSquared;
ulong current = destination[0];
for (int i = 1; i < oddPowerCount; i++)
{
current = current.MontgomeryMultiply(square, modulus, nPrime);
destination[i] = current;
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong[] InitializeMontgomeryOddPowersGpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
private static ulong[] InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
{
ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount];
oddPowers[0] = divisor.MontgomeryTwo;
@@ -630,31 +481,17 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
{
ulong rotationCount = exponent % cycleLength;
return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false);
return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
{
return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false);
return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
{
ulong rotationCount = exponent % cycleLength;
return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
{
return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
{
-18
View File
@@ -1,18 +0,0 @@
# GPU Modular Arithmetic in the `--mersenne=bydivisor` Path
The GPU implementation now routes single-use exponent checks through standard modular arithmetic while keeping Montgomery reduction only where repeated multiplications reuse the same modulus.
## Kernels using standard modular arithmetic
- `CheckKernel` in `MersenneNumberDivisorByDivisorGpuTester` evaluates each divisor candidate with `Pow2ModBinaryGpu` and compares the result directly against 1. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L561-L645】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
- `ComputeExponentKernel` in the same tester precomputes batched residues with the plain helper so CPU lookups receive standard-domain results. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L780-L791】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
- `ComputeResiduesKernel` inside `MersenneNumberDivisorResidueGpuEvaluator` now mirrors that logic when building residue tables for mixed CPU/GPU scans. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs†L102-L114】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
- The order-search helpers in `PrimeOrderGpuHeuristics` (`TrySpecialMaxKernel`, `ExponentLoweringKernel`, and `Pow2EqualsOneKernel`) also rely on `Pow2ModBinaryGpu` for their equality checks. 【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L798-L866】【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L1596-L1598】
## Kernels still using Montgomery arithmetic
- `MontgomeryOddPowerKernelScan` in `GpuKernelPool` continues to build windowed odd-power tables with `MontgomeryMultiply`, which amortizes the Montgomery setup over many multiplications per modulus. 【F:PerfectNumbers.Core/Gpu/GpuKernelPool.cs†L200-L237】【F:PerfectNumbers.Core/ULongExtensions.cs†L532-L559】
## Shared helper
- `Pow2ModBinaryGpu` in `ULongExtensions` implements the GPU-friendly binary exponentiation loop that powers the standard modular path across these kernels. 【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
@@ -8,7 +8,4 @@
<ItemGroup>
<ProjectReference Include="../../PerfectNumbers.Core/PerfectNumbers.Core.csproj" />
</ItemGroup>
<ItemGroup>
<Compile Remove=".temp.cs" />
</ItemGroup>
</Project>
</Project>