diff --git a/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs
index 090c163..238943e 100644
--- a/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs
@@ -50,7 +50,7 @@ public class GpuPow2ModBenchmarks
return Pow2ModSingleBit(Input.Exponent, Input.Modulus);
}
- private static GpuUInt128 Pow2ModSingleBit(ulong exponent, in GpuUInt128 modulus)
+ private static GpuUInt128 Pow2ModSingleBit(ulong exponent, GpuUInt128 modulus)
{
if (modulus.IsZero || modulus == GpuUInt128.One)
{
diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs
index 48c2050..2504a3e 100644
--- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs
@@ -52,14 +52,31 @@ 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.
///
[Benchmark(Baseline = true)]
- public ulong Baseline()
+ public ulong BaselineCpu()
{
GetData(out var exponents, out var divisors, out _);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
- checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false);
+ checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false);
+ }
+
+ return checksum;
+ }
+
+ ///
+ /// Baseline right-to-left Montgomery ladder; measured ??.?? μs on the small sample set and ??.?? μs on the large one.
+ ///
+ [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);
}
return checksum;
@@ -69,14 +86,31 @@ 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).
///
[Benchmark]
- public ulong BaselineWithKnownCycle()
+ public ulong BaselineWithKnownCycleCpu()
{
GetData(out var exponents, out var divisors, out var cycles);
ulong checksum = 0UL;
for (int i = 0; i < SampleCount; i++)
{
- checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
+ checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]);
+ }
+
+ return checksum;
+ }
+
+ ///
+ /// 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).
+ ///
+ [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]);
}
return checksum;
diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs
index 8c773fb..1342a97 100644
--- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs
@@ -22,320 +22,391 @@ 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);
+ }
+ }
- ///
- /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
- /// large set.
- ///
- // [Benchmark(Baseline = true)]
- public ulong MontgomeryWithoutCycle()
- {
- GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
- ulong checksum = 0UL;
+ ///
+ /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
+ /// large set.
+ ///
+ // [Benchmark(Baseline = true)]
+ public ulong MontgomeryWithoutCycleCpu()
+ {
+ GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+ ulong checksum = 0UL;
- for (int i = 0; i < SampleCount; i++)
- {
- checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false);
- }
+ for (int i = 0; i < SampleCount; i++)
+ {
+ checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false);
+ }
- return checksum;
- }
+ return checksum;
+ }
- ///
- /// 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×).
- ///
- [Benchmark]
- public ulong MontgomeryWithPrecomputedCycle()
- {
- GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
- ulong checksum = 0UL;
+ ///
+ /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
+ /// large set.
+ ///
+ // [Benchmark(Baseline = true)]
+ public ulong MontgomeryWithoutCycleGpu()
+ {
+ GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+ ulong checksum = 0UL;
- for (int i = 0; i < SampleCount; i++)
- {
- checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
- }
+ for (int i = 0; i < SampleCount; i++)
+ {
+ checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false);
+ }
- return checksum;
- }
+ return checksum;
+ }
- ///
- /// 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).
- ///
- // [Benchmark]
- public ulong MontgomeryWithGpuCycleComputation()
- {
- GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
- ulong checksum = 0UL;
+ ///
+ /// 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×).
+ ///
+ [Benchmark]
+ public ulong MontgomeryWithPrecomputedCycleCpu()
+ {
+ GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
+ ulong checksum = 0UL;
- for (int i = 0; i < SampleCount; i++)
- {
- ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
- checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
- }
+ for (int i = 0; i < SampleCount; i++)
+ {
+ checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]);
+ }
- return checksum;
- }
+ return checksum;
+ }
+
+ ///
+ /// 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×).
+ ///
+ [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;
+ }
+
+ ///
+ /// 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).
+ ///
+ // [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;
+ }
+
+ ///
+ /// 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).
+ ///
+ // [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 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 buffer = stackalloc byte[8];
+ _random.NextBytes(buffer);
+ ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
+ return value;
+ }
+
+ private UInt128 NextWideOddModulus()
+ {
+ Span buffer = stackalloc byte[16];
+ _random.NextBytes(buffer);
+ ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
+ Span 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);
+ }
- [Benchmark]
- public ulong MontgomeryWithHeuristicCycleComputation()
- {
- GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
- ulong checksum = 0UL;
+ 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);
+ }
- for (int i = 0; i < SampleCount; i++)
- {
- ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
- checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
- }
+ // When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
- 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 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 buffer = stackalloc byte[8];
- _random.NextBytes(buffer);
- ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
- return value;
- }
-
- private UInt128 NextWideOddModulus()
- {
- Span buffer = stackalloc byte[16];
- _random.NextBytes(buffer);
- ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
- Span 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);
- // }
+ // 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
- PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
- modulus,
- _previousPrimeOrder,
- divisorData,
- PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+ ulong order = PrimeOrderCalculator.Calculate(
+ modulus,
+ _previousPrimeOrder,
+ divisorData,
+ PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- if (orderResult.Order != 0UL)
- {
- _previousPrimeOrder = orderResult.Order;
- return orderResult.Order;
- }
+ if (order != 0UL)
+ {
+ _previousPrimeOrder = order;
+ return 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
- PrimeOrderCalculator.PrimeOrderResultWide orderResult = PrimeOrderCalculator.Calculate(
- modulus,
- _previousWidePrimeOrder,
- PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+ UInt128 order = PrimeOrderCalculator.Calculate(
+ modulus,
+ _previousWidePrimeOrder,
+ PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- if (orderResult.Order != UInt128.Zero)
- {
- _previousWidePrimeOrder = orderResult.Order;
- return orderResult.Order;
- }
+ if (order != UInt128.Zero)
+ {
+ _previousWidePrimeOrder = order;
+ return 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);
+ }
}
diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs
index 440cd96..d3ad626 100644
--- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs
@@ -176,7 +176,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
return (UInt128)result;
}
- private static bool ShouldUseSingleBit(in GpuUInt128 exponent)
+ private static bool ShouldUseSingleBit(GpuUInt128 exponent)
{
if (exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold)
{
@@ -257,7 +257,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
return result;
}
- private static bool IsBitSet(in GpuUInt128 value, int bitIndex)
+ private static bool IsBitSet(GpuUInt128 value, int bitIndex)
{
if (bitIndex >= 64)
{
@@ -267,7 +267,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
return ((value.Low >> bitIndex) & 1UL) != 0UL;
}
- private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount)
+ private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount)
{
if (windowStart != 0)
{
diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs
index 38b3933..3909631 100644
--- a/PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs
+++ b/PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs
@@ -15,9 +15,15 @@ public class GpuContextPoolTests
{
GpuContextPool.DisposeAll();
- using var lease = GpuContextPool.RentPreferred(preferCpu: true);
-
- lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
+ GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
+ try
+ {
+ lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
+ }
+ finally
+ {
+ lease.Dispose();
+ }
GpuContextPool.DisposeAll();
}
@@ -30,17 +36,27 @@ public class GpuContextPoolTests
Context firstContext;
Accelerator firstAccelerator;
- using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
+ GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
+ try
{
firstContext = lease.Context;
firstAccelerator = lease.Accelerator;
}
+ finally
+ {
+ lease.Dispose();
+ }
- using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
+ lease = GpuContextPool.RentPreferred(preferCpu: true);
+ try
{
lease.Context.Should().BeSameAs(firstContext);
lease.Accelerator.Should().BeSameAs(firstAccelerator);
}
+ finally
+ {
+ lease.Dispose();
+ }
GpuContextPool.DisposeAll();
}
@@ -52,17 +68,27 @@ public class GpuContextPoolTests
GpuContextPool.DisposeAll();
Context firstContext;
- using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
+ GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
+ try
{
firstContext = lease.Context;
}
+ finally
+ {
+ lease.Dispose();
+ }
GpuContextPool.DisposeAll();
- using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
+ lease = GpuContextPool.RentPreferred(preferCpu: true);
+ try
{
lease.Context.Should().NotBeSameAs(firstContext);
}
+ finally
+ {
+ lease.Dispose();
+ }
GpuContextPool.DisposeAll();
}
diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs
index 9294c87..dd1d39d 100644
--- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs
+++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs
@@ -199,10 +199,17 @@ public class GpuKernelCompilationTests
[MemberData(nameof(KernelLoaders))]
public void Kernel_compiles_without_internal_compiler_errors(string name, Action compile)
{
- using var lease = GpuContextPool.RentPreferred(preferCpu: true);
- var accelerator = lease.Accelerator;
- Action action = () => compile(accelerator);
- action.Should().NotThrow($"kernel {name} should 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();
+ }
}
private static void CompileDivisorCycleCacheKernel(Accelerator accelerator)
diff --git a/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs b/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs
index 24fb2d0..fa078b0 100644
--- a/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs
+++ b/PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs
@@ -150,29 +150,36 @@ public class NttGpuMathTests
var modulus = new GpuUInt128(0UL, 17UL);
var primitiveRoot = new GpuUInt128(0UL, 3UL);
- 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)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
- var twiddlesInvMont = (MemoryBuffer1D)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++)
+ GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
+ try
{
- ulong prod = (ulong)((((BigInteger)forward[i].Low) * inverse[i].Low) % mod);
- prod.Should().Be(r2);
+ 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)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
+ var twiddlesInvMont = (MemoryBuffer1D)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();
}
}
@@ -338,12 +345,23 @@ public class NttGpuMathTests
values[1] = new GpuUInt128(0UL, 2UL);
var modulus = new GpuUInt128(0UL, 17UL);
var primitiveRoot = new GpuUInt128(0UL, 3UL);
- using (var gpu = GpuContextPool.Rent())
+ GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
+ try
{
var buffer = gpu.Accelerator.Allocate1D(values.Length);
- buffer.View.CopyFromCPU(ref values[0], values.Length);
- NttGpuMath.SquareDevice(gpu.Accelerator, buffer.View, modulus, primitiveRoot);
- buffer.Dispose();
+ try
+ {
+ buffer.View.CopyFromCPU(ref values[0], values.Length);
+ NttGpuMath.SquareDevice(gpu.Accelerator, buffer.View, modulus, primitiveRoot);
+ }
+ finally
+ {
+ buffer.Dispose();
+ }
+ }
+ finally
+ {
+ gpu.Dispose();
}
// Verify that the square cache now holds entries for the accelerator.
diff --git a/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs b/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs
index 84398a1..57fe816 100644
--- a/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs
+++ b/PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs
@@ -33,16 +33,23 @@ public class NttMontgomeryTests
NttGpuMath.Square(cpu, modulus, primitiveRoot);
// Device execution
- using var context = GpuContextPool.Rent();
- var accelerator = context.Accelerator;
- using var buffer = accelerator.Allocate1D(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++)
+ GpuContextPool.GpuContextLease context = GpuContextPool.Rent();
+ try
{
- gpu[i].Should().Be(cpu[i]);
+ var accelerator = context.Accelerator;
+ using var buffer = accelerator.Allocate1D(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
+ {
+ context.Dispose();
}
}
diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs
index 3338cd0..c381c39 100644
--- a/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs
+++ b/PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs
@@ -359,10 +359,9 @@ public class PrimeOrderGpuHeuristicsTests
PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
- PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+ ulong result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
- result.Order.Should().Be(3UL);
+ result.Should().Be(3UL);
}
finally
{
diff --git a/PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs b/PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs
index 1c8fe35..c2bdd5b 100644
--- a/PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs
+++ b/PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs
@@ -11,14 +11,21 @@ public class GpuKernelPoolTests
[Trait("Category", "Fast")]
public void EnsureSmallCyclesOnDevice_returns_same_view_on_subsequent_calls()
{
- using var lease = GpuContextPool.RentPreferred(preferCpu: true);
- var accelerator = lease.Accelerator;
+ GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
+ try
+ {
+ 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);
+ second.Equals(first).Should().BeTrue();
+ second.Length.Should().Be(MersenneDivisorCycles.SmallDivisorsMax + 1);
+ }
+ finally
+ {
+ lease.Dispose();
+ }
}
[Fact]
diff --git a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs
index d7dd522..843131c 100644
--- a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs
+++ b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs
@@ -45,14 +45,14 @@ public class PrimeOrderCalculatorTests
{
ulong expected = ComputeOrderByDoubling(prime);
- PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
+ ulong result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
- PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+ PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
- result.Order.Should().Be(expected);
+ result.Should().Be(expected);
}
[Theory]
@@ -62,21 +62,21 @@ public class PrimeOrderCalculatorTests
{
ulong previousOrder = ComputeOrderByDoubling(previousPrime);
- PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
+ ulong heuristic = PrimeOrderCalculator.Calculate(
prime,
previousOrder,
MontgomeryDivisorData.FromModulus(prime),
- PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+ PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
+ ulong strict = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
- PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault);
+ PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
- strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
- heuristic.Order.Should().Be(strict.Order);
+ heuristic.Should().Be(strict);
}
[Theory]
@@ -85,14 +85,14 @@ public class PrimeOrderCalculatorTests
[InlineData(3UL, 2UL)]
public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder)
{
- PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
+ ulong result = PrimeOrderCalculator.Calculate(
prime,
previousOrder: null,
MontgomeryDivisorData.FromModulus(prime),
- PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+ PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
- result.Order.Should().Be(expectedOrder);
+ result.Should().Be(expectedOrder);
}
[Fact]
@@ -105,14 +105,14 @@ public class PrimeOrderCalculatorTests
maxPowChecks: 8,
mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
- PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
+ ulong result = PrimeOrderCalculator.Calculate(
prime: 13UL,
previousOrder: null,
MontgomeryDivisorData.FromModulus(13UL),
- config);
+ config,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
- result.Order.Should().Be(ComputeOrderByDoubling(13UL));
+ result.Should().Be(ComputeOrderByDoubling(13UL));
}
[Fact]
@@ -126,14 +126,14 @@ public class PrimeOrderCalculatorTests
mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
- PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
+ ulong heuristic = PrimeOrderCalculator.Calculate(
prime: 239UL,
previousOrder: null,
divisorData,
- heuristicConfig);
+ heuristicConfig,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
- heuristic.Order.Should().Be(ComputeOrderByDoubling(239UL));
+ heuristic.Should().Be(ComputeOrderByDoubling(239UL));
var strictConfig = new PrimeOrderCalculator.PrimeOrderSearchConfig(
smallFactorLimit: 2,
@@ -141,14 +141,14 @@ public class PrimeOrderCalculatorTests
maxPowChecks: 1,
mode: PrimeOrderCalculator.PrimeOrderMode.Strict);
- PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
+ ulong strict = PrimeOrderCalculator.Calculate(
prime: 239UL,
previousOrder: null,
divisorData,
- strictConfig);
+ strictConfig,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
- strict.Order.Should().Be(heuristic.Order);
+ strict.Should().Be(heuristic);
}
[Fact]
@@ -157,13 +157,13 @@ public class PrimeOrderCalculatorTests
{
UInt128 prime = UInt128.Parse("18446744073709641691");
- PrimeOrderCalculator.PrimeOrderResultWide result = PrimeOrderCalculator.Calculate(
- prime,
- previousOrder: null,
- PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+ UInt128 result = PrimeOrderCalculator.Calculate(
+ prime,
+ previousOrder: null,
+ PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+ PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
- result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
- result.Order.Should().Be((UInt128)1229782938247309446UL);
+ result.Should().Be((UInt128)1229782938247309446UL);
}
[Fact]
diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
index ec96d17..49baab4 100644
--- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
+++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
@@ -238,20 +238,19 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
}
}
- if (divisorCycle != 0 && divisorCycle == prime)
- {
- processedAll = true;
- return true;
- }
- else if (divisorCycle == 0)
- {
- Console.WriteLine($"Divisor cycle wasn't calculated for ${prime}");
- }
- // if (divisorCycle != 0UL && CheckDivisor(prime, divisorCycle, divisorData) != 0)
- // {
- // processedAll = true;
- // return true;
- // }
+ if (divisorCycle == prime)
+ {
+ // A cycle equal to the tested exponent (which is prime in this path)
+ // guarantees that the candidate divides the corresponding Mersenne
+ // number because the order of 2 modulo the divisor is exactly p.
+ processedAll = true;
+ return true;
+ }
+
+ if (divisorCycle == 0UL)
+ {
+ Console.WriteLine($"Divisor cycle was not calculated for {prime}");
+ }
}
divisor += step;
@@ -269,7 +268,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
{
- ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData);
+ ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
return residue == 1UL ? (byte)1 : (byte)0;
}
diff --git a/PerfectNumbers.Core/DeterministicRandom.cs b/PerfectNumbers.Core/DeterministicRandom.cs
new file mode 100644
index 0000000..b889d31
--- /dev/null
+++ b/PerfectNumbers.Core/DeterministicRandom.cs
@@ -0,0 +1,49 @@
+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;
+ }
+}
diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepper.cs
index 39fb669..0129904 100644
--- a/PerfectNumbers.Core/ExponentRemainderStepper.cs
+++ b/PerfectNumbers.Core/ExponentRemainderStepper.cs
@@ -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.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
+ _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_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.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
+ ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_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.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
+ _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_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.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
+ ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_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.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
+ _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_previousExponent = exponent;
_hasState = true;
isUnity = _currentMontgomery == _montgomeryOne;
diff --git a/PerfectNumbers.Core/FactorEntry.cs b/PerfectNumbers.Core/FactorEntry.cs
new file mode 100644
index 0000000..e4ba3be
--- /dev/null
+++ b/PerfectNumbers.Core/FactorEntry.cs
@@ -0,0 +1,8 @@
+namespace PerfectNumbers.Core
+{
+ public readonly struct FactorEntry(ulong value, int exponent)
+ {
+ public readonly ulong Value = value;
+ public readonly int Exponent = exponent;
+ }
+}
\ No newline at end of file
diff --git a/PerfectNumbers.Core/Gpu/GpuContextPool.cs b/PerfectNumbers.Core/Gpu/GpuContextPool.cs
index 91a6c54..1361e01 100644
--- a/PerfectNumbers.Core/Gpu/GpuContextPool.cs
+++ b/PerfectNumbers.Core/Gpu/GpuContextPool.cs
@@ -126,7 +126,7 @@ public static class GpuContextPool
}
}
- public struct GpuContextLease : IDisposable
+ public struct GpuContextLease
{
private readonly PooledContext _ctx;
private bool disposed;
diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs
index 9000d2b..6b022aa 100644
--- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs
+++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs
@@ -677,7 +677,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
// GC.SuppressFinalize(this);
}
- public readonly struct ExecutionScope : IDisposable
+ public readonly struct ExecutionScope
{
private readonly object? _lock;
diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs
index cee3beb..e3e64f8 100644
--- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs
+++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs
@@ -267,7 +267,7 @@ public struct GpuUInt128 : IComparable, IEquatable
public void Add(in ReadOnlyGpuUInt128 other) => AddInternal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void Add(in GpuUInt128 other) => AddInternal(other.High, other.Low);
+ public void Add(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, IEquatable
public void AddMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void AddMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
+ public void AddMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void AddMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus)
+ public void AddMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
AddModInternal(value.High, value.Low, in readOnlyMod);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void AddMod(in GpuUInt128 value, in GpuUInt128 modulus)
+ public void AddMod(GpuUInt128 value, GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
AddModInternal(value.High, value.Low, in readOnlyMod);
@@ -390,7 +390,7 @@ public struct GpuUInt128 : IComparable, IEquatable
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void AddMod(in GpuUInt128 value, ulong modulus) => AddMod(value.Low, modulus);
+ public void AddMod(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, IEquatable
public void Sub(in ReadOnlyGpuUInt128 other) => SubInternal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void Sub(in GpuUInt128 other) => SubInternal(other.High, other.Low);
+ public void Sub(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, IEquatable
public void Mul64(in ReadOnlyGpuUInt128 other) => Mul64Internal(other.High, other.Low);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void Mul64(in GpuUInt128 other) => Mul64Internal(other.High, other.Low);
+ public void Mul64(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, IEquatable
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static int TrailingZeroCount(in GpuUInt128 value)
+ public static int TrailingZeroCount(GpuUInt128 value)
{
if (value.Low != 0UL)
{
@@ -646,17 +646,17 @@ public struct GpuUInt128 : IComparable, IEquatable
public void SubMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void SubMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
+ public void SubMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void SubMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus)
+ public void SubMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
SubModInternal(value.High, value.Low, in readOnlyMod);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void SubMod(in GpuUInt128 value, in GpuUInt128 modulus)
+ public void SubMod(GpuUInt128 value, GpuUInt128 modulus)
{
ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
SubModInternal(value.High, value.Low, in readOnlyMod);
@@ -727,7 +727,7 @@ public struct GpuUInt128 : IComparable, IEquatable
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void Mul(in GpuUInt128 other)
+ public void Mul(GpuUInt128 other)
{
ReadOnlyGpuUInt128 readOnlyThis = AsReadOnly();
ReadOnlyGpuUInt128 readOnlyOther = other.AsReadOnly();
diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
index 5a41f3f..2220711 100644
--- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
@@ -551,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
}
ulong exponent = exponents[index];
- hits[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+ hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
}
public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -770,7 +770,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
ulong exponent = exponents[index];
- results[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: true);
+ results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
}
public ulong DivisorLimit
diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
index fa4251e..d38c0cc 100644
--- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
+++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs
@@ -415,9 +415,9 @@ internal static partial class PrimeOrderGpuHeuristics
ulong? previousOrder,
PrimeOrderCalculator.PrimeOrderSearchConfig config,
in MontgomeryDivisorData divisorData,
- out PrimeOrderCalculator.PrimeOrderResult result)
+ out ulong order)
{
- result = default;
+ order = 0UL;
var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
var execution = lease.EnterExecutionScope();
@@ -497,18 +497,17 @@ internal static partial class PrimeOrderGpuHeuristics
throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity.");
}
- ulong order = 0UL;
resultBuffer.View.CopyToCPU(ref order, 1);
- PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found
- ? PrimeOrderCalculator.PrimeOrderStatus.Found
- : PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved;
+ if (kernelStatus == PrimeOrderKernelStatus.FactoringFailure)
+ {
+ order = 0UL;
+ }
- result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order);
DisposeResources();
lease.Dispose();
- return true;
+ return order != 0UL;
void DisposeResources()
{
@@ -812,7 +811,7 @@ internal static partial class PrimeOrderGpuHeuristics
}
ulong reduced = phi / factor;
- if (reduced.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL)
+ if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL)
{
return false;
}
@@ -1596,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics
private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
{
- return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL;
+ return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
}
private static ulong CalculateByDoublingKernel(ulong prime)
@@ -1796,7 +1795,7 @@ internal static partial class PrimeOrderGpuHeuristics
return 0UL;
}
- return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false);
+ return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
}
private static Action, MontgomeryDivisorData, ArrayView1D> GetPow2ModKernel(Accelerator accelerator)
@@ -1812,7 +1811,7 @@ internal static partial class PrimeOrderGpuHeuristics
private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders)
{
ulong exponent = exponents[index];
- remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
+ remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false);
}
private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator)
diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs
index ee77260..410978a 100644
--- a/PerfectNumbers.Core/MersenneDivisorCycles.cs
+++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs
@@ -102,7 +102,7 @@ public class MersenneDivisorCycles
return cycleLength == exponent;
}
- return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+ return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
}
public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
@@ -126,7 +126,7 @@ public class MersenneDivisorCycles
}
}
- return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+ return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
}
private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
@@ -234,14 +234,14 @@ public class MersenneDivisorCycles
return (UInt128)cycle;
}
- PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
+ UInt128 wideOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
- if (wideResult.Order != UInt128.Zero)
+ if (wideOrder != UInt128.Zero)
{
- return wideResult.Order;
+ return wideOrder;
}
// Otherwise, find order of 2 mod divisor
@@ -285,15 +285,15 @@ public class MersenneDivisorCycles
return true;
}
- PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
+ ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
- if (orderResult.Order != 0UL)
+ if (computedOrder != 0UL)
{
- cycleLength = orderResult.Order;
+ cycleLength = computedOrder;
return true;
}
@@ -505,7 +505,7 @@ public class MersenneDivisorCycles
}
ulong candidate = order / prime;
- if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+ if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL)
{
order = candidate;
continue;
@@ -940,15 +940,15 @@ public class MersenneDivisorCycles
if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
{
- PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
+ ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor,
previousOrder: null,
divisorData,
PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
- if (orderResult.Order != 0UL)
+ if (computedOrder != 0UL)
{
- cycleLength = orderResult.Order;
+ cycleLength = computedOrder;
return true;
}
}
diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs
index 1080414..f480e12 100644
--- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs
+++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs
@@ -371,10 +371,6 @@ 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)
{
@@ -383,8 +379,8 @@ public static class MersennePrimeFactorTester
return 0UL;
}
- c = (ulong)rng.Next(2, int.MaxValue);
- x = (ulong)rng.Next(2, int.MaxValue);
+ c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+ x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
y = x;
d = 1UL;
@@ -523,9 +519,6 @@ 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)
{
@@ -534,8 +527,8 @@ public static class MersennePrimeFactorTester
return 0;
}
- c = (UInt128)rng.NextInt64(2, long.MaxValue);
- x = (UInt128)rng.NextInt64(2, long.MaxValue);
+ c = (DeterministicRandom.NextUInt128() % (n - 1)) + 1;
+ x = (DeterministicRandom.NextUInt128() % (n - 2)) + 2;
y = x;
d = 1;
diff --git a/PerfectNumbers.Core/PartialFactorResult.cs b/PerfectNumbers.Core/PartialFactorResult.cs
new file mode 100644
index 0000000..dc9c1b0
--- /dev/null
+++ b/PerfectNumbers.Core/PartialFactorResult.cs
@@ -0,0 +1,81 @@
+using System;
+using System.Buffers;
+
+namespace PerfectNumbers.Core
+{
+ public sealed class PartialFactorResult
+ {
+ [ThreadStatic]
+ private static PartialFactorResult? s_poolHead;
+
+ public FactorEntry[]? Factors;
+ private PartialFactorResult? _next;
+ public ulong Cofactor;
+ public bool FullyFactored;
+ public int Count;
+
+ 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;
+ }
+
+ instance._next = null;
+ instance.Factors = factors;
+ instance.Cofactor = cofactor;
+ instance.FullyFactored = fullyFactored;
+ instance.Count = count;
+ return instance;
+ }
+
+ public static readonly PartialFactorResult Empty = new PartialFactorResult
+ {
+ Cofactor = 1UL,
+ FullyFactored = true,
+ };
+
+ public PartialFactorResult WithAdditionalPrime(ulong prime)
+ {
+ FactorEntry[]? source = Factors;
+ if (source is null || Count == 0)
+ {
+ FactorEntry[] local = ArrayPool.Shared.Rent(1);
+ local[0] = new FactorEntry(prime, 1);
+ return Rent(local, 1UL, true, 1);
+ }
+
+ FactorEntry[] extended = ArrayPool.Shared.Rent(Count + 1);
+ Array.Copy(source, 0, extended, 0, Count);
+ extended[Count] = new FactorEntry(prime, 1);
+ Span span = extended.AsSpan(0, Count + 1);
+ span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+ return Rent(extended, 1UL, true, Count + 1);
+ }
+
+ public void Dispose()
+ {
+ if (this != Empty)
+ {
+ FactorEntry[]? factors = Factors;
+ if (factors is not null)
+ {
+ Factors = null;
+ ArrayPool.Shared.Return(factors, clearArray: false);
+ }
+
+ _next = s_poolHead;
+ s_poolHead = this;
+ }
+ }
+ }
+}
diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs
index 4eb6e93..3e37bd3 100644
--- a/PerfectNumbers.Core/PerfectNumberConstants.cs
+++ b/PerfectNumbers.Core/PerfectNumberConstants.cs
@@ -8,6 +8,7 @@ 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
diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs
index 9e1fd27..fc06285 100644
--- a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs
+++ b/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs
@@ -9,30 +9,15 @@ namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
- 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(
+ public static UInt128 Calculate(
UInt128 prime,
UInt128? previousOrder,
PrimeOrderSearchConfig config,
PrimeOrderHeuristicDevice device)
{
- using var scope = UsePow2Mode(device);
+ var scope = UsePow2Mode(device);
MontgomeryDivisorData divisorData;
+ UInt128 result;
if (prime <= ulong.MaxValue)
{
ulong? previous = null;
@@ -51,41 +36,49 @@ internal static partial class PrimeOrderCalculator
ulong prime64 = (ulong)prime;
divisorData = MontgomeryDivisorData.FromModulus(prime64);
- PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device);
- return new PrimeOrderResultWide(result.Status, result.Order);
+ ulong order64 = Calculate(prime64, previous, divisorData, config, device);
+ result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
}
else
{
divisorData = default;
+ result = CalculateWideInternal(prime, previousOrder, divisorData, config);
}
- return CalculateWideInternal(prime, previousOrder, divisorData, config);
+ scope.Dispose();
+ return result;
}
- private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+ private static UInt128 CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
{
if (prime <= UInt128.One)
{
- return new PrimeOrderResultWide(PrimeOrderStatus.Found, UInt128.One);
+ return UInt128.One;
}
if (prime == (UInt128)3UL)
{
- return new PrimeOrderResultWide(PrimeOrderStatus.Found, (UInt128)2UL);
+ return (UInt128)2UL;
}
UInt128 phi = prime - UInt128.One;
- using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
+ var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
+ UInt128 result;
if (phiFactors.Factors is null)
{
- return FinishStrictlyWide(prime, divisorData, config.Mode);
+ result = FinishStrictlyWide(prime, divisorData);
+ }
+ else
+ {
+ result = RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
}
- return RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
+ debugScope.Dispose();
+ return result;
}
- private static PrimeOrderResultWide RunHeuristicPipelineWide(
+ private static UInt128 RunHeuristicPipelineWide(
UInt128 prime,
UInt128? previousOrder,
in MontgomeryDivisorData divisorData,
@@ -95,7 +88,7 @@ internal static partial class PrimeOrderCalculator
{
if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
{
- return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi);
+ return phi;
}
UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData);
@@ -103,27 +96,25 @@ internal static partial class PrimeOrderCalculator
if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config))
{
- return new PrimeOrderResultWide(PrimeOrderStatus.Found, candidateOrder);
+ return candidateOrder;
}
if (config.Mode == PrimeOrderMode.Strict)
{
- return FinishStrictlyWide(prime, divisorData, PrimeOrderMode.Strict);
+ return FinishStrictlyWide(prime, divisorData);
}
if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
{
- return new PrimeOrderResultWide(PrimeOrderStatus.Found, order);
+ return order;
}
- return FinishStrictlyWide(prime, divisorData, config.Mode);
+ return candidateOrder;
}
- private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData, PrimeOrderMode mode)
+ private static UInt128 FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData)
{
- UInt128 strictOrder = CalculateByFactorizationWide(prime, divisorData);
- PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved;
- return new PrimeOrderResultWide(status, strictOrder);
+ return CalculateByFactorizationWide(prime, divisorData);
}
private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs
index f079c42..299cc2c 100644
--- a/PerfectNumbers.Core/PrimeOrderCalculator.cs
+++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs
@@ -1,65 +1,30 @@
-using System;
using System.Buffers;
using System.Buffers.Binary;
-using System.Collections.Generic;
using System.Diagnostics;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
-using System.Security.Cryptography;
using PerfectNumbers.Core.Gpu;
-using System.Threading;
namespace PerfectNumbers.Core;
internal static partial class PrimeOrderCalculator
{
- internal enum PrimeOrderStatus
- {
- Found,
- HeuristicUnresolved,
- }
-
- internal readonly struct PrimeOrderResult
- {
- public PrimeOrderResult(PrimeOrderStatus status, ulong order)
- {
- Status = status;
- Order = order;
- }
-
- public PrimeOrderStatus Status { get; }
-
- public ulong Order { get; }
- }
-
internal enum PrimeOrderMode
{
Heuristic,
Strict,
}
- internal readonly struct PrimeOrderSearchConfig
+ internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode)
{
- public PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode)
- {
- SmallFactorLimit = smallFactorLimit;
- PollardRhoMilliseconds = pollardRhoMilliseconds;
- MaxPowChecks = maxPowChecks;
- Mode = mode;
- }
+ public static PrimeOrderSearchConfig HeuristicDefault => new(smallFactorLimit: 100_000, pollardRhoMilliseconds: 128, maxPowChecks: 24, PrimeOrderMode.Heuristic);
+ public static PrimeOrderSearchConfig StrictDefault => new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict);
- public uint SmallFactorLimit { get; }
-
- public int PollardRhoMilliseconds { get; }
-
- public int MaxPowChecks { get; }
-
- public PrimeOrderMode Mode { get; }
-
- public static PrimeOrderSearchConfig HeuristicDefault => new(100_000, 128, 24, PrimeOrderMode.Heuristic);
-
- public static PrimeOrderSearchConfig StrictDefault => new(1_000_000, 0, 0, PrimeOrderMode.Strict);
+ public readonly uint SmallFactorLimit = smallFactorLimit;
+ public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds;
+ public readonly int MaxPowChecks = maxPowChecks;
+ public readonly PrimeOrderMode Mode = mode;
}
internal enum PrimeOrderHeuristicDevice
@@ -76,15 +41,24 @@ internal static partial class PrimeOrderCalculator
private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
-
[ThreadStatic]
private static PrimeOrderHeuristicDevice s_deviceMode;
[ThreadStatic]
private static bool s_debugLoggingEnabled;
+ [ThreadStatic]
+ private static List? s_partialFactorPending;
+
+ [ThreadStatic]
+ private static Stack? s_partialFactorCompositeStack;
+
+ [ThreadStatic]
+ private static Dictionary? s_partialFactorCounts;
+
private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu;
- private readonly struct Pow2ModeScope : IDisposable
+
+ private readonly struct Pow2ModeScope
{
private readonly bool _previousInitialized;
private readonly bool _previousAllow;
@@ -108,7 +82,7 @@ internal static partial class PrimeOrderCalculator
}
}
- private readonly struct DebugLoggingScope : IDisposable
+ private readonly struct DebugLoggingScope
{
private readonly bool _previous;
@@ -146,93 +120,88 @@ internal static partial class PrimeOrderCalculator
}
}
- public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
- => Calculate(prime, previousOrder, divisorData, config, PrimeOrderHeuristicDevice.Gpu);
-
- public static PrimeOrderResult Calculate(
+ public static ulong Calculate(
ulong prime,
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
- PrimeOrderSearchConfig config,
+ in PrimeOrderSearchConfig config,
PrimeOrderHeuristicDevice device)
{
- using var scope = UsePow2Mode(device);
- return CalculateInternal(prime, previousOrder, divisorData, config);
+ var scope = UsePow2Mode(device);
+ ulong order = CalculateInternal(prime, previousOrder, divisorData, config);
+ scope.Dispose();
+ return order;
}
- private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+ private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
{
if (prime <= 3UL)
{
- return new PrimeOrderResult(PrimeOrderStatus.Found, prime == 3UL ? 2UL : 1UL);
+ return prime == 3UL ? 2UL : 1UL;
}
ulong phi = prime - 1UL;
- if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out PrimeOrderResult gpuResult))
+ if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder))
{
- return gpuResult;
+ return gpuOrder;
}
- using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
- DebugLog("Partial factoring φ(p)");
-
PartialFactorResult phiFactors = PartialFactor(phi, config);
+ ulong result;
if (phiFactors.Factors is null)
{
- DebugLog("No factors found");
- return FinishStrictly(prime, config.Mode);
+ result = FinishStrictlyCpu(prime, divisorData);
+ phiFactors.Dispose();
+ return result;
}
- return RunHeuristicPipeline(prime, previousOrder, config, divisorData, phi, phiFactors);
+ result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors);
+ phiFactors.Dispose();
+ return result;
}
- private static PrimeOrderResult RunHeuristicPipeline(
+ private static ulong RunHeuristicPipelineCpu(
ulong prime,
ulong? previousOrder,
- PrimeOrderSearchConfig config,
+ in PrimeOrderSearchConfig config,
in MontgomeryDivisorData divisorData,
ulong phi,
PartialFactorResult phiFactors)
{
- DebugLog("Trying special max check");
- if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData))
+ if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData))
{
- return new PrimeOrderResult(PrimeOrderStatus.Found, phi);
+ return phi;
}
- DebugLog("Initializing starting order");
- ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData);
- candidateOrder = ExponentLowering(candidateOrder, prime, phiFactors, divisorData);
+ ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData);
+ candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData);
- DebugLog("Trying to confirm order");
- if (TryConfirmOrder(prime, candidateOrder, divisorData, config))
+ if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config))
{
- return new PrimeOrderResult(PrimeOrderStatus.Found, candidateOrder);
+ return candidateOrder;
}
if (config.Mode == PrimeOrderMode.Strict)
{
- return FinishStrictly(prime, PrimeOrderMode.Strict);
+ return FinishStrictlyCpu(prime, divisorData);
}
- if (TryHeuristicFinish(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order))
+ if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order))
{
- return new PrimeOrderResult(PrimeOrderStatus.Found, order);
+ return order;
}
- DebugLog("Heuristic unresolved, finishing strictly");
- return FinishStrictly(prime, config.Mode);
+ return candidateOrder;
}
- private static PrimeOrderResult FinishStrictly(ulong prime, PrimeOrderMode mode)
+ private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData)
{
- ulong strictOrder = CalculateByDoubling(prime);
- return new PrimeOrderResult(mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved, strictOrder);
+ return CalculateByFactorizationCpu(prime, divisorData);
}
- private static bool TrySpecialMax(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
+ private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
{
ReadOnlySpan factorSpan = factors.Factors;
int length = factors.Count;
@@ -240,7 +209,7 @@ internal static partial class PrimeOrderCalculator
{
ulong factor = factorSpan[i].Value;
ulong reduced = phi / factor;
- if (Pow2EqualsOne(reduced, prime, divisorData))
+ if (Pow2EqualsOneCpu(reduced, prime, divisorData))
{
return false;
}
@@ -249,13 +218,13 @@ internal static partial class PrimeOrderCalculator
return true;
}
- private static ulong InitializeStartingOrder(ulong prime, ulong phi, in MontgomeryDivisorData divisorData)
+ private static ulong InitializeStartingOrderCpu(ulong prime, ulong phi, in MontgomeryDivisorData divisorData)
{
ulong order = phi;
if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL)
{
ulong half = phi >> 1;
- if (Pow2EqualsOne(half, prime, divisorData))
+ if (Pow2EqualsOneCpu(half, prime, divisorData))
{
order = half;
}
@@ -264,7 +233,7 @@ internal static partial class PrimeOrderCalculator
return order;
}
- private static ulong ExponentLowering(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
+ private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
{
FactorEntry[]? tempArray = null;
try
@@ -297,7 +266,7 @@ internal static partial class PrimeOrderCalculator
}
ulong reduced = order / primeFactor;
- if (Pow2EqualsOne(reduced, prime, divisorData))
+ if (Pow2EqualsOneCpu(reduced, prime, divisorData))
{
order = reduced;
continue;
@@ -318,7 +287,7 @@ internal static partial class PrimeOrderCalculator
}
}
- private static bool TryConfirmOrder(ulong prime, ulong order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+ private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
{
if (order == 0UL)
{
@@ -326,73 +295,80 @@ internal static partial class PrimeOrderCalculator
}
// Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`.
- DebugLog("Verifying a^order ≡ 1 (mod p)");
- if (!Pow2EqualsOne(order, prime, divisorData))
+ // DebugLog("Verifying a^order ≡ 1 (mod p)");
+ if (!Pow2EqualsOneCpu(order, prime, divisorData))
{
return false;
}
- DebugLog("Partial factoring order");
+ // DebugLog("Partial factoring order");
// TODO: Do we do partial factoring of order multiple times?
PartialFactorResult factorization = PartialFactor(order, config);
- if (factorization.Factors is null)
+ try
{
- return false;
- }
-
- if (!factorization.FullyFactored)
- {
- if (factorization.Cofactor <= 1UL)
- {
- DebugLog("Cofactor <= 1. No factors found");
- return false;
- }
-
- // TODO: Use Open.Numerics.Primality for this final check once it's available.
- DebugLog("Cofactor > 1. Testing primality of cofactor");
- bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
- // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
- if (!isPrime)
+ if (factorization.Factors is null)
{
return false;
}
- DebugLog("Adding cofactor as prime factor");
- factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
- }
-
- ReadOnlySpan span = factorization.Factors;
- DebugLog("Verifying prime-power reductions");
- int length = factorization.Count;
- for (int i = 0; i < length; i++)
- {
- ulong primeFactor = span[i].Value;
- ulong reduced = order;
- for (int iteration = 0; iteration < span[i].Exponent; iteration++)
+ if (!factorization.FullyFactored)
{
- if ((reduced % primeFactor) != 0UL)
- {
- break;
- }
-
- reduced /= primeFactor;
- if (Pow2EqualsOne(reduced, prime, divisorData))
+ if (factorization.Cofactor <= 1UL)
{
return false;
}
- }
- }
- return true;
+ // DebugLog("Cofactor > 1. Testing primality of cofactor");
+ bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
+ // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
+ if (!isPrime)
+ {
+ return false;
+ }
+
+ // DebugLog("Adding cofactor as prime factor");
+ PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
+ factorization.Dispose();
+ factorization = extended;
+ }
+
+ ReadOnlySpan span = factorization.Factors;
+ // DebugLog("Verifying prime-power reductions");
+ int length = factorization.Count;
+ for (int i = 0; i < length; i++)
+ {
+ ulong primeFactor = span[i].Value;
+ ulong reduced = order;
+ for (int iteration = 0; iteration < span[i].Exponent; iteration++)
+ {
+ if ((reduced % primeFactor) != 0UL)
+ {
+ break;
+ }
+
+ reduced /= primeFactor;
+ if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+ {
+ return false;
+ }
+ }
+ }
+
+ return true;
+ }
+ finally
+ {
+ factorization.Dispose();
+ }
}
- private static bool TryHeuristicFinish(
+ private static bool TryHeuristicFinishCpu(
ulong prime,
ulong order,
ulong? previousOrder,
in MontgomeryDivisorData divisorData,
- PrimeOrderSearchConfig config,
+ in PrimeOrderSearchConfig config,
PartialFactorResult phiFactors,
out ulong result)
{
@@ -403,129 +379,142 @@ internal static partial class PrimeOrderCalculator
}
// TODO: Do we do partial factoring of order multiple times?
- DebugLog("Trying heuristic. Partial factoring order");
+ // DebugLog("Trying heuristic. Partial factoring order");
PartialFactorResult orderFactors = PartialFactor(order, config);
- if (orderFactors.Factors is null)
+ try
{
- return false;
- }
-
- if (!orderFactors.FullyFactored)
- {
- if (orderFactors.Cofactor <= 1UL)
+ if (orderFactors.Factors is null)
{
return false;
}
- var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor);
- // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None))
- if (!isPrime)
+ if (!orderFactors.FullyFactored)
{
- return false;
- }
-
- orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
- }
-
- int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4;
- List candidates = new(capacity);
- FactorEntry[] factorArray = orderFactors.Factors!;
- DebugLog("Building candidates list");
- BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
- if (candidates.Count == 0)
- {
- return false;
- }
-
- DebugLog("Sorting candidates");
- SortCandidates(prime, previousOrder, candidates);
-
- int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
- int powUsed = 0;
- int candidateCount = candidates.Count;
- bool allowGpuBatch = true;
- Span candidateSpan = CollectionsMarshal.AsSpan(candidates);
-
- DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
- int index = 0;
- const int MaxGpuBatchSize = 1024;
- const int StackGpuBatchSize = 1024;
- Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize];
- while (index < candidateCount && powUsed < powBudget)
- {
- int remaining = candidateCount - index;
- int budgetRemaining = powBudget - powUsed;
- int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize));
- if (batchSize <= 0)
- {
- break;
- }
-
- ReadOnlySpan batch = candidateSpan.Slice(index, batchSize);
- ulong[]? gpuPool = null;
- Span pooledGpuRemainders = default;
- bool gpuSuccess = false;
- bool gpuStackRemainders = false;
- GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
-
- if (allowGpuBatch && IsGpuPow2Allowed)
- {
- if (batchSize <= StackGpuBatchSize)
+ if (orderFactors.Cofactor <= 1UL)
{
- Span localRemainders = stackGpuRemainders.Slice(0, batchSize);
- status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData);
- if (status == GpuPow2ModStatus.Success)
- {
- gpuSuccess = true;
- gpuStackRemainders = true;
- }
+ return false;
}
- else
+
+ var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor);
+ // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None))
+ if (!isPrime)
{
- gpuPool = ArrayPool.Shared.Rent(batchSize);
- Span pooledRemainders = gpuPool.AsSpan(0, batchSize);
- status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData);
- if (status == GpuPow2ModStatus.Success)
+ return false;
+ }
+
+ PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
+ orderFactors.Dispose();
+ orderFactors = extended;
+ }
+
+ int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4;
+ List candidates = new(capacity);
+ FactorEntry[] factorArray = orderFactors.Factors!;
+ // DebugLog("Building candidates list");
+ BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
+ if (candidates.Count == 0)
+ {
+ return false;
+ }
+
+ // DebugLog("Sorting candidates");
+ SortCandidates(prime, previousOrder, candidates);
+
+ int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
+ int powUsed = 0;
+ int candidateCount = candidates.Count;
+ bool allowGpuBatch = true;
+ Span candidateSpan = CollectionsMarshal.AsSpan(candidates);
+
+ // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
+ int index = 0;
+ const int MaxGpuBatchSize = 1024;
+ const int StackGpuBatchSize = 1024;
+ Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize];
+ while (index < candidateCount && powUsed < powBudget)
+ {
+ int remaining = candidateCount - index;
+ int budgetRemaining = powBudget - powUsed;
+ int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize));
+ if (batchSize <= 0)
+ {
+ break;
+ }
+
+ ReadOnlySpan batch = candidateSpan.Slice(index, batchSize);
+ ulong[]? gpuPool = null;
+ Span pooledGpuRemainders = default;
+ bool gpuSuccess = false;
+ bool gpuStackRemainders = false;
+ GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
+
+ if (allowGpuBatch && IsGpuPow2Allowed)
+ {
+ if (batchSize <= StackGpuBatchSize)
{
- pooledGpuRemainders = pooledRemainders;
- gpuSuccess = true;
+ Span localRemainders = stackGpuRemainders.Slice(0, batchSize);
+ status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData);
+ if (status == GpuPow2ModStatus.Success)
+ {
+ gpuSuccess = true;
+ gpuStackRemainders = true;
+ }
}
else
{
- ArrayPool.Shared.Return(gpuPool, clearArray: false);
- gpuPool = null;
+ gpuPool = ArrayPool.Shared.Rent(batchSize);
+ Span pooledRemainders = gpuPool.AsSpan(0, batchSize);
+ status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData);
+ if (status == GpuPow2ModStatus.Success)
+ {
+ pooledGpuRemainders = pooledRemainders;
+ gpuSuccess = true;
+ }
+ else
+ {
+ ArrayPool.Shared.Return(gpuPool, clearArray: false);
+ gpuPool = null;
+ }
+ }
+
+ if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable))
+ {
+ allowGpuBatch = false;
}
}
- if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable))
+ for (int i = 0; i < batchSize && powUsed < powBudget; i++)
{
- allowGpuBatch = false;
- }
- }
+ ulong candidate = batch[i];
+ powUsed++;
- for (int i = 0; i < batchSize && powUsed < powBudget; i++)
- {
- ulong candidate = batch[i];
- powUsed++;
+ bool equalsOne;
+ if (gpuSuccess)
+ {
+ ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i];
+ equalsOne = remainderValue == 1UL;
+ }
+ else
+ {
+ equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData);
+ }
+ if (!equalsOne)
+ {
+ continue;
+ }
- bool equalsOne;
- if (gpuSuccess)
- {
- ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i];
- equalsOne = remainderValue == 1UL;
- }
- else
- {
- equalsOne = Pow2EqualsOne(candidate, prime, divisorData);
- }
- if (!equalsOne)
- {
- continue;
- }
+ if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget))
+ {
+ continue;
+ }
- if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget))
- {
- continue;
+ if (gpuPool is not null)
+ {
+ ArrayPool.Shared.Return(gpuPool, clearArray: false);
+ }
+
+ result = candidate;
+ return true;
}
if (gpuPool is not null)
@@ -533,20 +522,16 @@ internal static partial class PrimeOrderCalculator
ArrayPool.Shared.Return(gpuPool, clearArray: false);
}
- result = candidate;
- return true;
+ index += batchSize;
}
- if (gpuPool is not null)
- {
- ArrayPool.Shared.Return(gpuPool, clearArray: false);
- }
-
- index += batchSize;
+ // DebugLog("No candidate confirmed");
+ return false;
+ }
+ finally
+ {
+ orderFactors.Dispose();
}
-
- DebugLog("No candidate confirmed");
- return false;
}
private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates)
@@ -716,62 +701,71 @@ internal static partial class PrimeOrderCalculator
}
}
- private static bool TryConfirmCandidate(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, ref int powUsed, int powBudget)
+ private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget)
{
PartialFactorResult factorization = PartialFactor(candidate, config);
- if (factorization.Factors is null)
+ try
{
- return false;
- }
-
- if (!factorization.FullyFactored)
- {
- if (factorization.Cofactor <= 1UL)
+ if (factorization.Factors is null)
{
return false;
}
- bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
- // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
- if (!isPrime)
+ if (!factorization.FullyFactored)
{
- return false;
+ if (factorization.Cofactor <= 1UL)
+ {
+ return false;
+ }
+
+ bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
+ // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
+ if (!isPrime)
+ {
+ return false;
+ }
+
+ PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
+ factorization.Dispose();
+ factorization = extended;
}
- factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
- }
+ ReadOnlySpan span = factorization.Factors;
+ int length = factorization.Count;
+ for (int i = 0; i < length; i++)
+ {
+ ulong primeFactor = span[i].Value;
+ ulong reduced = candidate;
+ for (int iteration = 0; iteration < span[i].Exponent; iteration++)
+ {
+ if ((reduced % primeFactor) != 0UL)
+ {
+ break;
+ }
- ReadOnlySpan span = factorization.Factors;
- int length = factorization.Count;
- for (int i = 0; i < length; i++)
+ reduced /= primeFactor;
+ if (powUsed >= powBudget && powBudget > 0)
+ {
+ return false;
+ }
+
+ powUsed++;
+ if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+ {
+ return false;
+ }
+ }
+ }
+
+ return true;
+ }
+ finally
{
- ulong primeFactor = span[i].Value;
- ulong reduced = candidate;
- for (int iteration = 0; iteration < span[i].Exponent; iteration++)
- {
- if ((reduced % primeFactor) != 0UL)
- {
- break;
- }
-
- reduced /= primeFactor;
- if (powUsed >= powBudget && powBudget > 0)
- {
- return false;
- }
-
- powUsed++;
- if (Pow2EqualsOne(reduced, prime, divisorData))
- {
- return false;
- }
- }
+ factorization.Dispose();
}
-
- return true;
}
- private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
+ private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
{
if (IsGpuPow2Allowed)
{
@@ -782,10 +776,59 @@ internal static partial class PrimeOrderCalculator
}
}
- return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+ return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
}
- private static PartialFactorResult PartialFactor(ulong value, PrimeOrderSearchConfig config)
+ private static List AcquirePendingCompositeList(int capacityHint)
+ {
+ List? list = s_partialFactorPending;
+ if (list is null)
+ {
+ list = new List(capacityHint);
+ s_partialFactorPending = list;
+ return list;
+ }
+
+ list.Clear();
+ if (list.Capacity < capacityHint)
+ {
+ list.Capacity = capacityHint;
+ }
+
+ return list;
+ }
+
+ private static Stack AcquireCompositeStack(int capacityHint)
+ {
+ Stack? stack = s_partialFactorCompositeStack;
+ if (stack is null)
+ {
+ stack = new Stack(capacityHint);
+ s_partialFactorCompositeStack = stack;
+ return stack;
+ }
+
+ stack.Clear();
+ stack.EnsureCapacity(capacityHint);
+ return stack;
+ }
+
+ private static Dictionary AcquireFactorCountDictionary(int capacityHint)
+ {
+ Dictionary? dictionary = s_partialFactorCounts;
+ if (dictionary is null)
+ {
+ dictionary = new Dictionary(capacityHint);
+ s_partialFactorCounts = dictionary;
+ return dictionary;
+ }
+
+ dictionary.Clear();
+ dictionary.EnsureCapacity(capacityHint);
+ return dictionary;
+ }
+
+ private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config)
{
if (value <= 1UL)
{
@@ -803,166 +846,165 @@ internal static partial class PrimeOrderCalculator
bool useDictionary = false;
uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
- ulong remaining;
+ ulong remaining = value;
bool gpuFactored = false;
if (IsGpuHeuristicDevice)
{
gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _);
}
- else
- {
- remaining = value;
- }
- if (!gpuFactored)
+ List pending = AcquirePendingCompositeList(2);
+ Stack? compositeStack = null;
+
+ try
{
- primeSlots.Clear();
- exponentSlots.Clear();
- if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining))
+ if (!gpuFactored)
{
- useDictionary = true;
- counts = new Dictionary(capacity: 8);
- remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
- }
- }
-
- List pending = new();
- if (remaining > 1UL)
- {
- pending.Add(remaining);
- }
-
- if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
- {
- DebugLog("Processing pending composites with Pollard's Rho");
- Stopwatch stopwatch = Stopwatch.StartNew();
- long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks;
- Stack stack = new();
- stack.Push(remaining);
- pending.Clear();
-
- while (stack.Count > 0)
- {
- ulong composite = stack.Pop();
- if (composite <= 1UL)
+ primeSlots.Clear();
+ exponentSlots.Clear();
+ if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining))
{
- continue;
+ useDictionary = true;
+ counts = AcquireFactorCountDictionary(Math.Max(factorCount, 8));
+ remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
}
+ }
+ if (remaining > 1UL)
+ {
+ pending.Add(remaining);
+ }
+
+ if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
+ {
+ Stopwatch stopwatch = Stopwatch.StartNew();
+ long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks;
+ compositeStack = AcquireCompositeStack(Math.Max(pending.Count * 2, 4));
+ compositeStack.Push(remaining);
+ pending.Clear();
+
+ while (compositeStack.Count > 0)
+ {
+ ulong composite = compositeStack.Pop();
+ if (composite <= 1UL)
+ {
+ continue;
+ }
+
+ bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
+ if (isPrime)
+ {
+ AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
+ continue;
+ }
+
+ if (stopwatch.ElapsedTicks > budgetTicks)
+ {
+ pending.Add(composite);
+ continue;
+ }
+
+ if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor))
+ {
+ pending.Add(composite);
+ continue;
+ }
+
+ ulong quotient = composite / factor;
+ compositeStack.Push(factor);
+ compositeStack.Push(quotient);
+ }
+ }
+
+ ulong cofactor = 1UL;
+ int pendingCount = pending.Count;
+ for (int i = 0; i < pendingCount; i++)
+ {
+ ulong composite = pending[i];
bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
if (isPrime)
{
AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
- continue;
}
-
- if (stopwatch.ElapsedTicks > budgetTicks)
+ else if (composite > 1UL)
{
- pending.Add(composite);
- continue;
+ cofactor = checked(cofactor * composite);
}
+ }
- if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor))
+ if (useDictionary)
+ {
+ if ((counts is null || counts.Count == 0) && cofactor == value)
+ {
+ return PartialFactorResult.Rent(null, value, false, 0);
+ }
+
+ if (counts is null || counts.Count == 0)
+ {
+ return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
+ }
+
+ FactorEntry[] factors = ArrayPool.Shared.Rent(counts.Count);
+ int index = 0;
+ foreach (KeyValuePair entry in counts)
+ {
+ factors[index] = new FactorEntry(entry.Key, entry.Value);
+ index++;
+ }
+
+ Span span = factors.AsSpan(0, index);
+ span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+
+ bool fullyFactored = cofactor == 1UL;
+ return PartialFactorResult.Rent(factors, cofactor, fullyFactored, index);
+ }
+
+ int actualCount = 0;
+ for (int i = 0; i < factorCount; i++)
+ {
+ if (primeSlots[i] != 0UL && exponentSlots[i] != 0)
+ {
+ actualCount++;
+ }
+ }
+
+ if (actualCount == 0)
+ {
+ if (cofactor == value)
+ {
+ return PartialFactorResult.Rent(null, value, false, 0);
+ }
+
+ return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
+ }
+
+ FactorEntry[] array = ArrayPool.Shared.Rent(actualCount);
+ int arrayIndex = 0;
+ for (int i = 0; i < factorCount; i++)
+ {
+ ulong primeValue = primeSlots[i];
+ int exponentValue = exponentSlots[i];
+ if (primeValue == 0UL || exponentValue == 0)
{
- pending.Add(composite);
continue;
}
- ulong quotient = composite / factor;
- stack.Push(factor);
- stack.Push(quotient);
+ array[arrayIndex] = new FactorEntry(primeValue, exponentValue);
+ arrayIndex++;
}
- }
- ulong cofactor = 1UL;
- int pendingCount = pending.Count;
- for (int i = 0; i < pendingCount; i++)
+ Span arraySpan = array.AsSpan(0, arrayIndex);
+ arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+ bool fullyFactoredArray = cofactor == 1UL;
+ return PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex);
+ }
+ finally
{
- ulong composite = pending[i];
- bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
- if (isPrime)
- {
- AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1);
- }
- else if (composite > 1UL)
- {
- cofactor = checked(cofactor * composite);
- }
+ pending.Clear();
+ counts?.Clear();
+ compositeStack?.Clear();
}
-
- if (useDictionary)
- {
- if ((counts is null || counts.Count == 0) && cofactor == value)
- {
- DebugLog("cofactor is the same as value, no factors found");
- return new PartialFactorResult(null, value, false, 0);
- }
-
- if (counts is null || counts.Count == 0)
- {
- return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
- }
-
- FactorEntry[] factors = ArrayPool.Shared.Rent(counts.Count);
- int index = 0;
- DebugLog(() => $"Collecting {counts.Count} prime factors");
- foreach (KeyValuePair entry in counts)
- {
- factors[index] = new FactorEntry(entry.Key, entry.Value);
- index++;
- }
-
- Span span = factors.AsSpan(0, index);
- span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
-
- FactorEntry[] resultArray = new FactorEntry[index];
- span.CopyTo(resultArray);
- ArrayPool.Shared.Return(factors, clearArray: false);
-
- bool fullyFactored = cofactor == 1UL;
- return new PartialFactorResult(resultArray, cofactor, fullyFactored, index);
- }
-
- int actualCount = 0;
- for (int i = 0; i < factorCount; i++)
- {
- if (primeSlots[i] != 0UL && exponentSlots[i] != 0)
- {
- actualCount++;
- }
- }
-
- if (actualCount == 0)
- {
- if (cofactor == value)
- {
- DebugLog("cofactor is the same as value, no factors found");
- return new PartialFactorResult(null, value, false, 0);
- }
-
- return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
- }
-
- FactorEntry[] array = new FactorEntry[actualCount];
- int arrayIndex = 0;
- DebugLog(() => $"Collecting {actualCount} prime factors");
- for (int i = 0; i < factorCount; i++)
- {
- ulong primeValue = primeSlots[i];
- int exponentValue = exponentSlots[i];
- if (primeValue == 0UL || exponentValue == 0)
- {
- continue;
- }
-
- array[arrayIndex] = new FactorEntry(primeValue, exponentValue);
- arrayIndex++;
- }
-
- Array.Sort(array, static (a, b) => a.Value.CompareTo(b.Value));
- bool fullyFactoredArray = cofactor == 1UL;
- return new PartialFactorResult(array, cofactor, fullyFactoredArray, array.Length);
}
private const int GpuSmallPrimeFactorSlots = 64;
@@ -1114,8 +1156,6 @@ internal static partial class PrimeOrderCalculator
return true;
}
- Span buffer = stackalloc byte[8];
-
while (true)
{
if (stopwatch.ElapsedTicks > budgetTicks)
@@ -1123,10 +1163,8 @@ internal static partial class PrimeOrderCalculator
return false;
}
- RandomNumberGenerator.Fill(buffer);
- ulong c = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 1UL)) + 1UL;
- RandomNumberGenerator.Fill(buffer);
- ulong x = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 2UL)) + 2UL;
+ ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+ ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
ulong y = x;
ulong d = 1UL;
@@ -1254,7 +1292,7 @@ internal static partial class PrimeOrderCalculator
return;
}
- counts ??= new Dictionary(Math.Max(count, 8));
+ counts ??= AcquireFactorCountDictionary(Math.Max(count, 8));
CopyFactorsToDictionary(primes, exponents, count, counts);
count = 0;
useDictionary = true;
@@ -1280,59 +1318,98 @@ internal static partial class PrimeOrderCalculator
}
}
- private static ulong CalculateByDoubling(ulong prime)
+ private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData)
{
- ulong order = 1UL;
- ulong pow = 2UL;
-
- while (pow != 1UL)
+ ulong phi = prime - 1UL;
+ Dictionary counts = new(capacity: 8);
+ FactorCompletely(phi, counts);
+ if (counts.Count == 0)
{
- pow <<= 1;
- if (pow >= prime)
- {
- pow -= prime;
- }
+ return phi;
+ }
- order++;
+ List> entries = new(counts);
+ entries.Sort(static (a, b) => a.Key.CompareTo(b.Key));
+
+ ulong order = phi;
+ int entryCount = entries.Count;
+ for (int i = 0; i < entryCount; i++)
+ {
+ ulong primeFactor = entries[i].Key;
+ int exponent = entries[i].Value;
+ for (int iteration = 0; iteration < exponent; iteration++)
+ {
+ if ((order % primeFactor) != 0UL)
+ {
+ break;
+ }
+
+ ulong candidate = order / primeFactor;
+ if (Pow2EqualsOneCpu(candidate, prime, divisorData))
+ {
+ order = candidate;
+ continue;
+ }
+
+ break;
+ }
}
return order;
}
+ private static void FactorCompletely(ulong value, Dictionary counts)
+ {
+ if (value <= 1UL)
+ {
+ return;
+ }
+
+ if (Open.Numeric.Primes.Prime.Numbers.IsPrime(value))
+ {
+ AddFactor(counts, value, 1);
+ return;
+ }
+
+ ulong factor = PollardRhoStrict(value);
+ FactorCompletely(factor, counts);
+ FactorCompletely(value / factor, counts);
+ }
+
+ private static ulong PollardRhoStrict(ulong n)
+ {
+ if ((n & 1UL) == 0UL)
+ {
+ return 2UL;
+ }
+
+ while (true)
+ {
+ ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+ ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
+ ulong y = x;
+ ulong d = 1UL;
+
+ while (d == 1UL)
+ {
+ x = AdvancePolynomial(x, c, n);
+ y = AdvancePolynomial(y, c, n);
+ y = AdvancePolynomial(y, c, n);
+ ulong diff = x > y ? x - y : y - x;
+ d = BinaryGcd(diff, n);
+ }
+
+ if (d != n)
+ {
+ return d;
+ }
+ }
+ }
+
private readonly struct CandidateKey(int primary, long secondary, long tertiary)
{
public readonly int Primary = primary;
public readonly long Secondary = secondary;
public readonly long Tertiary = tertiary;
}
-
- private readonly struct FactorEntry(ulong value, int exponent)
- {
- public readonly ulong Value = value;
- public readonly int Exponent = exponent;
- }
-
- private readonly struct PartialFactorResult(FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count)
- {
- public readonly FactorEntry[]? Factors = factors;
- public readonly ulong Cofactor = cofactor;
- public readonly bool FullyFactored = fullyFactored;
- public readonly int Count = count;
- public readonly static PartialFactorResult Empty = new(null, 1UL, true, 0);
-
- public PartialFactorResult WithAdditionalPrime(ulong prime)
- {
- if (Factors is null)
- {
- FactorEntry[] local = [new FactorEntry(prime, 1)];
- return new PartialFactorResult(local, 1UL, true, 1);
- }
-
- FactorEntry[] extended = new FactorEntry[Count + 1];
- Array.Copy(Factors, 0, extended, 0, Count);
- extended[Count] = new FactorEntry(prime, 1);
- Array.Sort(extended, static (a, b) => a.Value.CompareTo(b.Value));
- return new PartialFactorResult(extended, 1UL, true, Count + 1);
- }
- }
}
diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs
index 7f62518..e7bf9c3 100644
--- a/PerfectNumbers.Core/UInt128Extensions.cs
+++ b/PerfectNumbers.Core/UInt128Extensions.cs
@@ -218,7 +218,7 @@ public static class UInt128Extensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- private static bool IsBitSet(in GpuUInt128 value, int bitIndex)
+ private static bool IsBitSet(GpuUInt128 value, int bitIndex)
{
if (bitIndex >= 64)
{
@@ -229,18 +229,16 @@ public static class UInt128Extensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount)
+ private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount)
{
+ ulong mask = (1UL << windowBitCount) - 1UL;
if (windowStart != 0)
{
- GpuUInt128 shifted = exponent;
- shifted.ShiftRight(windowStart);
- ulong mask = (1UL << windowBitCount) - 1UL;
- return shifted.Low & mask;
+ exponent.ShiftRight(windowStart);
+ return exponent.Low & mask;
}
- ulong directMask = (1UL << windowBitCount) - 1UL;
- return exponent.Low & directMask;
+ return exponent.Low & mask;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs
index f46b61c..fa0d030 100644
--- a/PerfectNumbers.Core/ULongExtensions.cs
+++ b/PerfectNumbers.Core/ULongExtensions.cs
@@ -1,4 +1,5 @@
using System;
+using System.Buffers;
using System.Numerics;
using System.Runtime.CompilerServices;
using ILGPU.Algorithms;
@@ -76,7 +77,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;
@@ -326,7 +327,7 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+ public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
{
ulong modulus = divisor.Modulus;
if (exponent == 0UL)
@@ -345,7 +346,78 @@ public static class ULongExtensions
ulong result = divisor.MontgomeryOne;
ulong nPrime = divisor.NPrime;
- ulong[] oddPowers = InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowerCount);
+ ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount);
+
+ 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 (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold)
+ {
+ ArrayPool.Shared.Return(oddPowers);
+ }
+
+ 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);
int index = bitLength - 1;
while (index >= 0)
@@ -461,7 +533,30 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- private static ulong[] InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
+ private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
+ {
+ ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold ? new ulong[oddPowerCount] : ArrayPool.Shared.Rent(oddPowerCount);
+
+ oddPowers[0] = divisor.MontgomeryTwo;
+ if (oddPowerCount == 1)
+ {
+ return oddPowers;
+ }
+
+ ulong previous,
+ square = divisor.MontgomeryTwoSquared;
+
+ for (int i = 1; i < oddPowerCount; i++)
+ {
+ previous = oddPowers[i - 1];
+ oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
+ }
+
+ return oddPowers;
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static ulong[] InitializeMontgomeryOddPowersGpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
{
ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount];
oddPowers[0] = divisor.MontgomeryTwo;
@@ -481,17 +576,31 @@ public static class ULongExtensions
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
+ public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
{
ulong rotationCount = exponent % cycleLength;
- return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false);
+ return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+ public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
{
- return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false);
+ return Pow2MontgomeryModWindowedCpu(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)
{