TODO added

Remove default prime order overloads

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)
         {

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.
     /// </summary>
     [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;
+    }
+
+    /// <summary>
+    /// Baseline right-to-left Montgomery ladder; measured ??.?? μs on the small sample set and ??.?? μs on the large one.
+    /// </summary>
+    [Benchmark]
+    public ulong BaselineGpu()
+    {
+        GetData(out var exponents, out var divisors, out _);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false);
         }
 
         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).
     /// </summary>
     [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;
+    }
+
+    /// <summary>
+    /// Uses a known divisor cycle to fold the exponent first, dropping runtimes to 26.04 μs on small inputs (0.72× baseline) and 9.53 μs on large inputs (0.08× baseline).
+    /// </summary>
+    [Benchmark]
+    public ulong BaselineWithKnownCycleGpu()
+    {
+        GetData(out var exponents, out var divisors, out var cycles);
+        ulong checksum = 0UL;
+
+        for (int i = 0; i < SampleCount; i++)
+        {
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]);
         }
 
         return checksum;

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);
+		}
+	}
 
-    /// <summary>
-    /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
-    /// large set.
-    /// </summary>
-    // [Benchmark(Baseline = true)]
-    public ulong MontgomeryWithoutCycle()
-    {
-        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
-        ulong checksum = 0UL;
+	/// <summary>
+	/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
+	/// large set.
+	/// </summary>
+	// [Benchmark(Baseline = true)]
+	public ulong MontgomeryWithoutCycleCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
 
-        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;
+	}
 
-    /// <summary>
-    /// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
-    /// on the large one (0.08×).
-    /// </summary>
-    [Benchmark]
-    public ulong MontgomeryWithPrecomputedCycle()
-    {
-        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
-        ulong checksum = 0UL;
+	/// <summary>
+	/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
+	/// large set.
+	/// </summary>
+	// [Benchmark(Baseline = true)]
+	public ulong MontgomeryWithoutCycleGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
 
-        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;
+	}
 
-    /// <summary>
-    /// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
-    /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
-    /// </summary>
-    // [Benchmark]
-    public ulong MontgomeryWithGpuCycleComputation()
-    {
-        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
-        ulong checksum = 0UL;
+	/// <summary>
+	/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
+	/// on the large one (0.08×).
+	/// </summary>
+	[Benchmark]
+	public ulong MontgomeryWithPrecomputedCycleCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
+		ulong checksum = 0UL;
 
-        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;
+	}
+
+	/// <summary>
+	/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
+	/// on the large one (0.08×).
+	/// </summary>
+	[Benchmark]
+	public ulong MontgomeryWithPrecomputedCycleGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	/// <summary>
+	/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
+	/// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
+	/// </summary>
+	// [Benchmark]
+	public ulong MontgomeryWithGpuCycleComputationCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	/// <summary>
+	/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
+	/// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
+	/// </summary>
+	// [Benchmark]
+	public ulong MontgomeryWithGpuCycleComputationGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	[Benchmark]
+	public ulong MontgomeryWithHeuristicCycleComputationCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	[Benchmark]
+	public ulong MontgomeryWithHeuristicCycleComputationGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	[Benchmark]
+	public UInt128 HeuristicWideCycleComputation()
+	{
+		UInt128 checksum = UInt128.Zero;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
+			_wideCycles[i] = cycle;
+			checksum ^= cycle;
+		}
+
+		return checksum;
+	}
+
+	private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
+	{
+		switch (Scale)
+		{
+			case InputScale.Small:
+				exponents = _smallExponents;
+				divisors = _smallDivisors;
+				cycles = _smallCycles;
+				break;
+			case InputScale.Large:
+				exponents = _largeExponents;
+				divisors = _largeDivisors;
+				cycles = _largeCycles;
+				break;
+			case InputScale.VeryLarge:
+				exponents = _veryLargeExponents;
+				divisors = _veryLargeDivisors;
+				cycles = _veryLargeCycles;
+				break;
+			default:
+				throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null);
+		}
+	}
+
+	private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
+
+	private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue);
+
+	private ulong NextVeryLargeExponent()
+	{
+		Span<byte> buffer = stackalloc byte[8];
+		_random.NextBytes(buffer);
+		ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
+		return value;
+	}
+
+	private ulong NextSmallOddModulus()
+	{
+		ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
+		return value | 1UL;
+	}
+
+	private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
+	{
+		int bitLength = _random.Next(48, 64);
+		ulong modulus = (1UL << bitLength) - 1UL;
+		return (modulus, (ulong)bitLength);
+	}
+
+	private ulong NextVeryLargeOddModulus()
+	{
+		Span<byte> buffer = stackalloc byte[8];
+		_random.NextBytes(buffer);
+		ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
+		return value;
+	}
+
+	private UInt128 NextWideOddModulus()
+	{
+		Span<byte> buffer = stackalloc byte[16];
+		_random.NextBytes(buffer);
+		ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
+		Span<byte> highSlice = buffer.Slice(8, 8);
+		ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
+		UInt128 value = ((UInt128)high << 64) | low;
+		return value | UInt128.One;
+	}
+
+	private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
+	{
+		if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+		{
+			return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL);
+		}
+
+		ulong nPrime = ComputeMontgomeryNPrime(modulus);
+		ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
+		ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
+		ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
+
+		return new MontgomeryDivisorData(
+			modulus,
+			nPrime,
+			montgomeryOne,
+			montgomeryTwo,
+			montgomeryTwoSquared);
+	}
 
 
-    [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<byte> buffer = stackalloc byte[8];
-        _random.NextBytes(buffer);
-        ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
-        return value;
-    }
-
-    private ulong NextSmallOddModulus()
-    {
-        ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
-        return value | 1UL;
-    }
-
-    private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
-    {
-        int bitLength = _random.Next(48, 64);
-        ulong modulus = (1UL << bitLength) - 1UL;
-        return (modulus, (ulong)bitLength);
-    }
-
-    private ulong NextVeryLargeOddModulus()
-    {
-        Span<byte> buffer = stackalloc byte[8];
-        _random.NextBytes(buffer);
-        ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
-        return value;
-    }
-
-    private UInt128 NextWideOddModulus()
-    {
-        Span<byte> buffer = stackalloc byte[16];
-        _random.NextBytes(buffer);
-        ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
-        Span<byte> highSlice = buffer.Slice(8, 8);
-        ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
-        UInt128 value = ((UInt128)high << 64) | low;
-        return value | UInt128.One;
-    }
-
-    private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
-    {
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL);
-        }
-
-        ulong nPrime = ComputeMontgomeryNPrime(modulus);
-        ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
-        ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
-        ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
-
-        return new MontgomeryDivisorData(
-            modulus,
-            nPrime,
-            montgomeryOne,
-            montgomeryTwo,
-            montgomeryTwoSquared);
-    }
-
-
-    private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
-    {
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
-        }
-
-        // When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
-
-        // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(modulus);
-        // // bool isPrime = PrimeTester.IsPrimeInternal(modulus, CancellationToken.None);
-        // if (!isPrime)
-        // {
-        //     return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
-        // }
+		// 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);
+	}
 }
 

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)
         {

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();
     }

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<Accelerator> 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)

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<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
-        var twiddlesInvMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesInvMont")!.GetValue(cache)!;
-        var r2 = (ulong)cacheType.GetProperty("MontR2Mod64")!.GetValue(cache)!;
-        var mod = (ulong)cacheType.GetProperty("ModulusLow")!.GetValue(cache)!;
-
-        int count = 7; // length 8 -> 7 twiddles
-        var forward = new GpuUInt128[count];
-        var inverse = new GpuUInt128[count];
-        twiddlesMont.CopyToCPU(forward);
-        twiddlesInvMont.CopyToCPU(inverse);
-
-        for (int i = 0; i < count; i++)
+        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<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
+            var twiddlesInvMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesInvMont")!.GetValue(cache)!;
+            var r2 = (ulong)cacheType.GetProperty("MontR2Mod64")!.GetValue(cache)!;
+            var mod = (ulong)cacheType.GetProperty("ModulusLow")!.GetValue(cache)!;
+
+            int count = 7; // length 8 -> 7 twiddles
+            var forward = new GpuUInt128[count];
+            var inverse = new GpuUInt128[count];
+            twiddlesMont.CopyToCPU(forward);
+            twiddlesInvMont.CopyToCPU(inverse);
+
+            for (int i = 0; i < count; i++)
+            {
+                ulong prod = (ulong)((((BigInteger)forward[i].Low) * inverse[i].Low) % mod);
+                prod.Should().Be(r2);
+            }
+        }
+        finally
+        {
+            gpu.Dispose();
         }
     }
 
@@ -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<GpuUInt128>(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.

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<GpuUInt128>(gpu.Length);
-        buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
-        NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
-        buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
-
-        for (int i = 0; i < cpu.Length; i++)
+        GpuContextPool.GpuContextLease context = GpuContextPool.Rent();
+        try
         {
-            gpu[i].Should().Be(cpu[i]);
+            var accelerator = context.Accelerator;
+            using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
+            buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
+            NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
+            buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
+
+            for (int i = 0; i < cpu.Length; i++)
+            {
+                gpu[i].Should().Be(cpu[i]);
+            }
+        }
+        finally
+        {
+            context.Dispose();
         }
     }
 

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
         {

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]

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]

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -14,8 +14,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
     private ulong _lastStatusDivisor;
     private bool _isConfigured;
     private int _batchSize = 1_024;
-    private const int PollardBrentMaxRestarts = 4;
-    private const int PollardBrentIterationBudget = 262_144;
 
 
     public int BatchSize
@@ -176,31 +174,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             return false;
         }
 
-        ulong pollardIterationCount = 0UL;
-        if (prime > 3UL)
-        {
-            var pollardBrent = new MersennePollardBrent(prime, _batchSize);
-            if (pollardBrent.IsSupported)
-            {
-                bool pollardFound = pollardBrent.TryFindFactor(
-                    allowedMax,
-                    PollardBrentMaxRestarts,
-                    PollardBrentIterationBudget,
-                    out ulong pollardFactor,
-                    out pollardIterationCount);
-
-                if (pollardFound)
-                {
-                    processedCount = pollardIterationCount;
-                    lastProcessed = pollardFactor;
-                    processedAll = true;
-                    return true;
-                }
-            }
-        }
-
-        processedCount = pollardIterationCount;
-
         UInt128 step = (UInt128)prime << 1;
         if (step == UInt128.Zero)
         {
@@ -265,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;
@@ -296,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;
     }
 

PerfectNumbers.Core/Cpu/MersennePollardBrent.cs

@@ -1,926 +0,0 @@
-using System;
-using System.Buffers;
-using System.Numerics;
-using System.Security.Cryptography;
-using PerfectNumbers.Core;
-using UInt128 = System.UInt128;
-
-namespace PerfectNumbers.Core.Cpu;
-
-internal sealed class MersennePollardBrent
-{
-    private enum ArithmeticMode
-    {
-        Unsupported,
-        UInt64,
-        UInt128,
-        BigInteger,
-    }
-
-    private readonly ArithmeticMode _mode;
-    private readonly int _bitLength;
-    private readonly int _byteLength;
-    private readonly int _bitRemainder;
-    private readonly byte _maskByte;
-    private readonly BigInteger _modulus;
-    private readonly BigInteger _mask;
-    private readonly ulong _modulus64;
-    private readonly ulong _mask64;
-    private readonly UInt128 _modulus128;
-    private readonly UInt128 _mask128;
-    private readonly int _batchSize;
-    private readonly bool _isSupported;
-
-    public MersennePollardBrent(ulong exponent, int batchSize)
-    {
-        if (exponent < 3UL || exponent > int.MaxValue)
-        {
-            _mode = ArithmeticMode.Unsupported;
-            _isSupported = false;
-            _bitLength = 0;
-            _byteLength = 0;
-            _bitRemainder = 0;
-            _maskByte = 0;
-            _modulus = BigInteger.Zero;
-            _mask = BigInteger.Zero;
-            _modulus64 = 0UL;
-            _mask64 = 0UL;
-            _modulus128 = UInt128.Zero;
-            _mask128 = UInt128.Zero;
-            _batchSize = 0;
-            return;
-        }
-
-        _isSupported = true;
-        _bitLength = (int)exponent;
-        _byteLength = (_bitLength + 7) / 8;
-        _bitRemainder = _bitLength & 7;
-        _maskByte = _bitRemainder == 0 ? (byte)0xFF : (byte)((1 << _bitRemainder) - 1);
-        _batchSize = Math.Clamp(batchSize, 8, 1024);
-
-        if (_bitLength <= 64)
-        {
-            _mode = ArithmeticMode.UInt64;
-            if (_bitLength == 64)
-            {
-                _modulus64 = ulong.MaxValue;
-                _mask64 = ulong.MaxValue;
-            }
-            else
-            {
-                _modulus64 = (1UL << _bitLength) - 1UL;
-                _mask64 = _modulus64;
-            }
-
-            _modulus128 = UInt128.Zero;
-            _mask128 = UInt128.Zero;
-            _modulus = BigInteger.Zero;
-            _mask = BigInteger.Zero;
-            return;
-        }
-
-        if (_bitLength <= 128)
-        {
-            _mode = ArithmeticMode.UInt128;
-            if (_bitLength == 128)
-            {
-                _modulus128 = UInt128.MaxValue;
-                _mask128 = UInt128.MaxValue;
-            }
-            else
-            {
-                _modulus128 = (UInt128.One << _bitLength) - UInt128.One;
-                _mask128 = _modulus128;
-            }
-
-            _modulus64 = 0UL;
-            _mask64 = 0UL;
-            _modulus = BigInteger.Zero;
-            _mask = BigInteger.Zero;
-            return;
-        }
-
-        _mode = ArithmeticMode.BigInteger;
-        _modulus = (BigInteger.One << _bitLength) - BigInteger.One;
-        _mask = _modulus;
-        _modulus64 = 0UL;
-        _mask64 = 0UL;
-        _modulus128 = UInt128.Zero;
-        _mask128 = UInt128.Zero;
-    }
-
-    public bool IsSupported => _isSupported;
-
-    public bool TryFindFactor(
-        ulong allowedMax,
-        int maxRestarts,
-        int maxIterationsPerRestart,
-        out ulong factor,
-        out ulong iterationCount)
-    {
-        factor = 0UL;
-        iterationCount = 0UL;
-
-        if (!_isSupported)
-        {
-            return false;
-        }
-
-        if (allowedMax <= 1UL)
-        {
-            return false;
-        }
-
-        return _mode switch
-        {
-            ArithmeticMode.UInt64 => TryFindFactorUInt64(allowedMax, maxRestarts, maxIterationsPerRestart, out factor, out iterationCount),
-            ArithmeticMode.UInt128 => TryFindFactorUInt128(allowedMax, maxRestarts, maxIterationsPerRestart, out factor, out iterationCount),
-            ArithmeticMode.BigInteger => TryFindFactorBigInteger(allowedMax, maxRestarts, maxIterationsPerRestart, out factor, out iterationCount),
-            _ => false,
-        };
-    }
-
-    private bool TryFindFactorUInt64(
-        ulong allowedMax,
-        int maxRestarts,
-        int maxIterationsPerRestart,
-        out ulong factor,
-        out ulong iterationCount)
-    {
-        factor = 0UL;
-        iterationCount = 0UL;
-
-        int restarts = Math.Max(1, maxRestarts);
-        int iterationBudget = Math.Max(_batchSize, maxIterationsPerRestart);
-
-        for (int restart = 0; restart < restarts; restart++)
-        {
-            ulong x = NextRandomSeed64();
-            ulong y = x;
-            ulong c = NextRandomConstant64();
-            ulong q = 1UL;
-            ulong g = 1UL;
-            ulong ys = 0UL;
-            int r = 1;
-            ulong localIterations = 0UL;
-
-            while (g == 1UL && localIterations < (ulong)iterationBudget)
-            {
-                x = y;
-                for (int i = 0; i < r && localIterations < (ulong)iterationBudget; i++)
-                {
-                    y = Iterate64(y, c, ref localIterations);
-                }
-
-                if (localIterations >= (ulong)iterationBudget)
-                {
-                    break;
-                }
-
-                int k = 0;
-                while (k < r && g == 1UL && localIterations < (ulong)iterationBudget)
-                {
-                    ys = y;
-                    int t = Math.Min(_batchSize, r - k);
-                    for (int i = 0; i < t; i++)
-                    {
-                        y = Iterate64(y, c, ref localIterations);
-                        ulong diff = AbsSub64(y, x);
-                        if (diff == 0UL)
-                        {
-                            continue;
-                        }
-
-                        q = MulMod64(q, diff);
-                    }
-
-                    g = BinaryGcd64(q, _modulus64);
-                    k += t;
-                }
-
-                r <<= 1;
-            }
-
-            if (g == _modulus64)
-            {
-                do
-                {
-                    ys = Iterate64(ys, c, ref localIterations);
-                    ulong diff = AbsSub64(ys, x);
-                    if (diff == 0UL)
-                    {
-                        continue;
-                    }
-
-                    g = BinaryGcd64(diff, _modulus64);
-                }
-                while (g == 1UL);
-            }
-
-            iterationCount += localIterations;
-
-            if (g == 1UL || g == _modulus64)
-            {
-                continue;
-            }
-
-            if (g > 1UL && g <= allowedMax)
-            {
-                factor = g;
-                return true;
-            }
-        }
-
-        return false;
-    }
-
-    private bool TryFindFactorUInt128(
-        ulong allowedMax,
-        int maxRestarts,
-        int maxIterationsPerRestart,
-        out ulong factor,
-        out ulong iterationCount)
-    {
-        factor = 0UL;
-        iterationCount = 0UL;
-
-        int restarts = Math.Max(1, maxRestarts);
-        int iterationBudget = Math.Max(_batchSize, maxIterationsPerRestart);
-        UInt128 allowedMaxWide = (UInt128)allowedMax;
-
-        for (int restart = 0; restart < restarts; restart++)
-        {
-            UInt128 x = NextRandomSeed128();
-            UInt128 y = x;
-            UInt128 c = NextRandomConstant128();
-            UInt128 q = UInt128.One;
-            UInt128 g = UInt128.One;
-            UInt128 ys = UInt128.Zero;
-            int r = 1;
-            ulong localIterations = 0UL;
-
-            while (g == UInt128.One && localIterations < (ulong)iterationBudget)
-            {
-                x = y;
-                for (int i = 0; i < r && localIterations < (ulong)iterationBudget; i++)
-                {
-                    y = Iterate128(y, c, ref localIterations);
-                }
-
-                if (localIterations >= (ulong)iterationBudget)
-                {
-                    break;
-                }
-
-                int k = 0;
-                while (k < r && g == UInt128.One && localIterations < (ulong)iterationBudget)
-                {
-                    ys = y;
-                    int t = Math.Min(_batchSize, r - k);
-                    for (int i = 0; i < t; i++)
-                    {
-                        y = Iterate128(y, c, ref localIterations);
-                        UInt128 diff = AbsSub128(y, x);
-                        if (diff == UInt128.Zero)
-                        {
-                            continue;
-                        }
-
-                        q = MulMod128(q, diff);
-                    }
-
-                    g = q.BinaryGcd(_modulus128);
-                    k += t;
-                }
-
-                r <<= 1;
-            }
-
-            if (g == _modulus128)
-            {
-                do
-                {
-                    ys = Iterate128(ys, c, ref localIterations);
-                    UInt128 diff = AbsSub128(ys, x);
-                    if (diff == UInt128.Zero)
-                    {
-                        continue;
-                    }
-
-                    g = diff.BinaryGcd(_modulus128);
-                }
-                while (g == UInt128.One);
-            }
-
-            iterationCount += localIterations;
-
-            if (g == UInt128.One || g == _modulus128)
-            {
-                continue;
-            }
-
-            if (g > UInt128.One && g <= allowedMaxWide)
-            {
-                factor = (ulong)g;
-                return true;
-            }
-        }
-
-        return false;
-    }
-
-    private bool TryFindFactorBigInteger(
-        ulong allowedMax,
-        int maxRestarts,
-        int maxIterationsPerRestart,
-        out ulong factor,
-        out ulong iterationCount)
-    {
-        factor = 0UL;
-        iterationCount = 0UL;
-
-        int restarts = Math.Max(1, maxRestarts);
-        int iterationBudget = Math.Max(_batchSize, maxIterationsPerRestart);
-
-        for (int restart = 0; restart < restarts; restart++)
-        {
-            BigInteger x = NextRandomSeedBigInteger();
-            BigInteger y = x;
-            BigInteger c = NextRandomConstantBigInteger();
-            BigInteger q = BigInteger.One;
-            BigInteger g = BigInteger.One;
-            BigInteger ys = BigInteger.Zero;
-            int r = 1;
-            ulong localIterations = 0UL;
-
-            while (g.IsOne && localIterations < (ulong)iterationBudget)
-            {
-                x = y;
-                for (int i = 0; i < r && localIterations < (ulong)iterationBudget; i++)
-                {
-                    y = IterateBigInteger(y, c, ref localIterations);
-                }
-
-                if (localIterations >= (ulong)iterationBudget)
-                {
-                    break;
-                }
-
-                int k = 0;
-                while (k < r && g.IsOne && localIterations < (ulong)iterationBudget)
-                {
-                    ys = y;
-                    int t = Math.Min(_batchSize, r - k);
-                    for (int i = 0; i < t; i++)
-                    {
-                        y = IterateBigInteger(y, c, ref localIterations);
-                        BigInteger diff = AbsSubBigInteger(y, x);
-                        if (diff.IsZero)
-                        {
-                            continue;
-                        }
-
-                        q = MulModBigInteger(q, diff);
-                    }
-
-                    g = BigInteger.GreatestCommonDivisor(q, _modulus);
-                    k += t;
-                }
-
-                r <<= 1;
-            }
-
-            if (g == _modulus)
-            {
-                do
-                {
-                    ys = IterateBigInteger(ys, c, ref localIterations);
-                    BigInteger diff = AbsSubBigInteger(ys, x);
-                    if (diff.IsZero)
-                    {
-                        continue;
-                    }
-
-                    g = BigInteger.GreatestCommonDivisor(diff, _modulus);
-                }
-                while (g.IsOne);
-            }
-
-            iterationCount += localIterations;
-
-            if (g.IsOne || g == _modulus)
-            {
-                continue;
-            }
-
-            if (g <= ulong.MaxValue)
-            {
-                ulong candidate = (ulong)g;
-                if (candidate > 1UL && candidate <= allowedMax)
-                {
-                    factor = candidate;
-                    return true;
-                }
-            }
-        }
-
-        return false;
-    }
-
-    private ulong Iterate64(ulong value, ulong c, ref ulong iterationCounter)
-    {
-        iterationCounter++;
-        value = SquareMod64(value);
-        value = AddMod64(value, c);
-        return value;
-    }
-
-    private UInt128 Iterate128(UInt128 value, UInt128 c, ref ulong iterationCounter)
-    {
-        iterationCounter++;
-        value = SquareMod128(value);
-        value = AddMod128(value, c);
-        return value;
-    }
-
-    private BigInteger IterateBigInteger(BigInteger value, BigInteger c, ref ulong iterationCounter)
-    {
-        iterationCounter++;
-        value = SquareModBigInteger(value);
-        value = AddModBigInteger(value, c);
-        return value;
-    }
-
-    private ulong AddMod64(ulong a, ulong b)
-    {
-        UInt128 sum = (UInt128)a + b;
-        return Reduce64(sum);
-    }
-
-    private UInt128 AddMod128(UInt128 a, UInt128 b)
-    {
-        UInt128 sum = a + b;
-        bool carryOut = sum < a;
-        if (!carryOut && sum < _modulus128)
-        {
-            return sum;
-        }
-
-        ulong p3 = 0UL;
-        ulong p2 = carryOut ? 1UL : 0UL;
-        ulong p1 = (ulong)(sum >> 64);
-        ulong p0 = (ulong)sum;
-        return ReduceProductBitwise(p3, p2, p1, p0, _modulus128);
-    }
-
-    private BigInteger AddModBigInteger(BigInteger a, BigInteger b)
-    {
-        return ReduceBigInteger(a + b);
-    }
-
-    private ulong SquareMod64(ulong value)
-    {
-        UInt128 square = (UInt128)value * value;
-        return Reduce64(square);
-    }
-
-    private UInt128 SquareMod128(UInt128 value)
-    {
-        return MulMod128(value, value);
-    }
-
-    private BigInteger SquareModBigInteger(BigInteger value)
-    {
-        return ReduceBigInteger(value * value);
-    }
-
-    private ulong MulMod64(ulong a, ulong b)
-    {
-        UInt128 product = (UInt128)a * b;
-        return Reduce64(product);
-    }
-
-    private UInt128 MulMod128(UInt128 a, UInt128 b)
-    {
-        MultiplyFull(a, b, out ulong p3, out ulong p2, out ulong p1, out ulong p0);
-        return ReduceProductBitwise(p3, p2, p1, p0, _modulus128);
-    }
-
-    private BigInteger MulModBigInteger(BigInteger a, BigInteger b)
-    {
-        if (b == _modulus)
-        {
-            return a;
-        }
-
-        if (a == _modulus)
-        {
-            a = BigInteger.Zero;
-        }
-
-        return ReduceBigInteger(a * b);
-    }
-
-    private ulong Reduce64(UInt128 value)
-    {
-        if (_bitLength == 64)
-        {
-            UInt128 combined = (UInt128)(ulong)value + (value >> 64);
-            UInt128 modulus = _modulus64;
-            if (combined >= modulus)
-            {
-                combined -= modulus;
-                if (combined >= modulus)
-                {
-                    combined -= modulus;
-                }
-            }
-
-            return (ulong)combined;
-        }
-
-        UInt128 mask = _mask64;
-        UInt128 modulusValue = _modulus64;
-        UInt128 working = value;
-        while (true)
-        {
-            UInt128 high = working >> _bitLength;
-            if (high == UInt128.Zero)
-            {
-                break;
-            }
-
-            working = (working & mask) + high;
-        }
-
-        if (working >= modulusValue)
-        {
-            working -= modulusValue;
-            if (working >= modulusValue)
-            {
-                working %= modulusValue;
-                if (working >= modulusValue)
-                {
-                    working -= modulusValue;
-                }
-            }
-        }
-
-        return (ulong)working;
-    }
-
-    private BigInteger ReduceBigInteger(BigInteger value)
-    {
-        if (value.Sign < 0)
-        {
-            value %= _modulus;
-            if (value.Sign < 0)
-            {
-                value += _modulus;
-            }
-
-            return value;
-        }
-
-        while (true)
-        {
-            BigInteger high = value >> _bitLength;
-            if (high.IsZero)
-            {
-                break;
-            }
-
-            value &= _mask;
-            value += high;
-        }
-
-        if (value >= _modulus)
-        {
-            value -= _modulus;
-            if (value >= _modulus)
-            {
-                value %= _modulus;
-                if (value.Sign < 0)
-                {
-                    value += _modulus;
-                }
-            }
-        }
-
-        return value;
-    }
-
-    private ulong NextRandomSeed64()
-    {
-        ulong value;
-        do
-        {
-            value = NextRandomBelowModulus64();
-        }
-        while (value <= 1UL);
-
-        return value;
-    }
-
-    private UInt128 NextRandomSeed128()
-    {
-        UInt128 value;
-        do
-        {
-            value = NextRandomBelowModulus128();
-        }
-        while (value <= UInt128.One);
-
-        return value;
-    }
-
-    private BigInteger NextRandomSeedBigInteger()
-    {
-        BigInteger value;
-        do
-        {
-            value = NextRandomBelowModulusBigInteger();
-        }
-        while (value <= BigInteger.One);
-
-        return value;
-    }
-
-    private ulong NextRandomConstant64()
-    {
-        ulong value;
-        do
-        {
-            value = NextRandomBelowModulus64();
-        }
-        while (value == 0UL);
-
-        return value;
-    }
-
-    private UInt128 NextRandomConstant128()
-    {
-        UInt128 value;
-        do
-        {
-            value = NextRandomBelowModulus128();
-        }
-        while (value == UInt128.Zero);
-
-        return value;
-    }
-
-    private BigInteger NextRandomConstantBigInteger()
-    {
-        BigInteger value;
-        do
-        {
-            value = NextRandomBelowModulusBigInteger();
-        }
-        while (value.IsZero);
-
-        return value;
-    }
-
-    private ulong NextRandomBelowModulus64()
-    {
-        if (_byteLength == 0)
-        {
-            return 0UL;
-        }
-
-        Span<byte> buffer = stackalloc byte[_byteLength];
-        ulong candidate;
-        do
-        {
-            RandomNumberGenerator.Fill(buffer);
-            if (_bitRemainder != 0)
-            {
-                buffer[_byteLength - 1] &= _maskByte;
-            }
-
-            candidate = (ulong)ReadUInt128LittleEndian(buffer);
-        }
-        while (candidate >= _modulus64);
-
-        return candidate;
-    }
-
-    private UInt128 NextRandomBelowModulus128()
-    {
-        if (_byteLength == 0)
-        {
-            return UInt128.Zero;
-        }
-
-        Span<byte> buffer = stackalloc byte[_byteLength];
-        UInt128 candidate;
-        do
-        {
-            RandomNumberGenerator.Fill(buffer);
-            if (_bitRemainder != 0)
-            {
-                buffer[_byteLength - 1] &= _maskByte;
-            }
-
-            candidate = ReadUInt128LittleEndian(buffer);
-        }
-        while (candidate >= _modulus128);
-
-        return candidate;
-    }
-
-    private BigInteger NextRandomBelowModulusBigInteger()
-    {
-        if (_byteLength == 0)
-        {
-            return BigInteger.Zero;
-        }
-
-        byte[] buffer = ArrayPool<byte>.Shared.Rent(_byteLength);
-        try
-        {
-            BigInteger value;
-            do
-            {
-                RandomNumberGenerator.Fill(buffer.AsSpan(0, _byteLength));
-                if (_bitRemainder != 0)
-                {
-                    buffer[_byteLength - 1] &= _maskByte;
-                }
-
-                value = new BigInteger(buffer, isUnsigned: true, isBigEndian: false);
-            }
-            while (value >= _modulus);
-
-            return value;
-        }
-        finally
-        {
-            ArrayPool<byte>.Shared.Return(buffer, clearArray: true);
-        }
-    }
-
-    private static ulong AbsSub64(ulong a, ulong b)
-    {
-        return a >= b ? a - b : b - a;
-    }
-
-    private static UInt128 AbsSub128(UInt128 a, UInt128 b)
-    {
-        return a >= b ? a - b : b - a;
-    }
-
-    private static BigInteger AbsSubBigInteger(BigInteger a, BigInteger b)
-    {
-        BigInteger diff = a - b;
-        return diff.Sign < 0 ? -diff : diff;
-    }
-
-    private static ulong BinaryGcd64(ulong u, ulong v)
-    {
-        if (u == 0UL)
-        {
-            return v;
-        }
-
-        if (v == 0UL)
-        {
-            return u;
-        }
-
-        int shift = BitOperations.TrailingZeroCount(u | v);
-        u >>= BitOperations.TrailingZeroCount(u);
-        do
-        {
-            v >>= BitOperations.TrailingZeroCount(v);
-            if (u > v)
-            {
-                (u, v) = (v, u);
-            }
-
-            v -= u;
-        }
-        while (v != 0UL);
-
-        return u << shift;
-    }
-
-    private static void MultiplyFull(UInt128 left, UInt128 right, out ulong p3, out ulong p2, out ulong p1, out ulong p0)
-    {
-        ulong leftLow = (ulong)left;
-        ulong leftHigh = (ulong)(left >> 64);
-        ulong rightLow = (ulong)right;
-        ulong rightHigh = (ulong)(right >> 64);
-
-        (ulong High, ulong Low) partial0 = Multiply64(leftLow, rightLow);
-        (ulong High, ulong Low) partial1 = Multiply64(leftLow, rightHigh);
-        (ulong High, ulong Low) partial2 = Multiply64(leftHigh, rightLow);
-        (ulong High, ulong Low) partial3 = Multiply64(leftHigh, rightHigh);
-
-        p0 = partial0.Low;
-        ulong carry = partial0.High;
-        ulong sum = carry + partial1.Low;
-        ulong carryMid = sum < partial1.Low ? 1UL : 0UL;
-        sum += partial2.Low;
-        if (sum < partial2.Low)
-        {
-            carryMid++;
-        }
-
-        p1 = sum;
-        sum = partial1.High + partial2.High;
-        ulong carryHigh = sum < partial2.High ? 1UL : 0UL;
-        sum += partial3.Low;
-        if (sum < partial3.Low)
-        {
-            carryHigh++;
-        }
-
-        sum += carryMid;
-        if (sum < carryMid)
-        {
-            carryHigh++;
-        }
-
-        p2 = sum;
-        p3 = partial3.High + carryHigh;
-    }
-
-    private static (ulong High, ulong Low) Multiply64(ulong left, ulong right)
-    {
-        ulong a0 = (uint)left;
-        ulong a1 = left >> 32;
-        ulong b0 = (uint)right;
-        ulong b1 = right >> 32;
-
-        ulong lo = a0 * b0;
-        ulong mid1 = a1 * b0;
-        b0 = a0 * b1;
-        b1 *= a1;
-
-        a0 = (lo >> 32) + (uint)mid1 + (uint)b0;
-        a1 = (lo & 0xFFFFFFFFUL) | (a0 << 32);
-        b1 += (mid1 >> 32) + (b0 >> 32) + (a0 >> 32);
-
-        return (b1, a1);
-    }
-
-    private static UInt128 ReduceProductBitwise(ulong p3, ulong p2, ulong p1, ulong p0, UInt128 modulus)
-    {
-        UInt128 remainder = UInt128.Zero;
-        for (int bit = 63; bit >= 0; bit--)
-        {
-            remainder <<= 1;
-            remainder |= (UInt128)((p3 >> bit) & 1UL);
-            if (remainder >= modulus)
-            {
-                remainder -= modulus;
-            }
-        }
-
-        for (int bit = 63; bit >= 0; bit--)
-        {
-            remainder <<= 1;
-            remainder |= (UInt128)((p2 >> bit) & 1UL);
-            if (remainder >= modulus)
-            {
-                remainder -= modulus;
-            }
-        }
-
-        for (int bit = 63; bit >= 0; bit--)
-        {
-            remainder <<= 1;
-            remainder |= (UInt128)((p1 >> bit) & 1UL);
-            if (remainder >= modulus)
-            {
-                remainder -= modulus;
-            }
-        }
-
-        for (int bit = 63; bit >= 0; bit--)
-        {
-            remainder <<= 1;
-            remainder |= (UInt128)((p0 >> bit) & 1UL);
-            if (remainder >= modulus)
-            {
-                remainder -= modulus;
-            }
-        }
-
-        return remainder;
-    }
-
-    private static UInt128 ReadUInt128LittleEndian(ReadOnlySpan<byte> buffer)
-    {
-        UInt128 value = UInt128.Zero;
-        for (int i = buffer.Length - 1; i >= 0; i--)
-        {
-            value <<= 8;
-            value |= buffer[i];
-        }
-
-        return value;
-    }
-}

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;
+    }
+}

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;

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;
+    }
+}

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;

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;
 

PerfectNumbers.Core/Gpu/GpuUInt128.cs

@@ -267,7 +267,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     public void Add(in ReadOnlyGpuUInt128 other) => AddInternal(other.High, other.Low);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Add(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<GpuUInt128>, IEquatable<GpuUInt128>
     public void AddMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void AddMod(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<GpuUInt128>, IEquatable<GpuUInt128>
     }
 
     [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<GpuUInt128>, IEquatable<GpuUInt128>
     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<GpuUInt128>, IEquatable<GpuUInt128>
     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<GpuUInt128>, IEquatable<GpuUInt128>
     }
 
     [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<GpuUInt128>, IEquatable<GpuUInt128>
     public void SubMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void SubMod(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<GpuUInt128>, IEquatable<GpuUInt128>
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Mul(in GpuUInt128 other)
+    public void Mul(GpuUInt128 other)
     {
         ReadOnlyGpuUInt128 readOnlyThis = AsReadOnly();
         ReadOnlyGpuUInt128 readOnlyOther = other.AsReadOnly();

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

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<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
@@ -1812,7 +1811,7 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong exponent = exponents[index];
-        remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
+        remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)

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;
             }
         }

PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs

@@ -112,9 +112,12 @@ public static class MersenneNumberDivisorByDivisorTester
         ulong maxPrime = 0UL;
         int primeWriteIndex = 0;
 
-        using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
-        bool hasPrime = primeEnumerator.MoveNext();
-        ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
+		// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
+		// We don't need to additionally check it.
+
+        // using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
+        // bool hasPrime = primeEnumerator.MoveNext();
+        // ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
 
         for (int index = startIndex; index < candidateCount; index++)
         {
@@ -127,25 +130,29 @@ public static class MersenneNumberDivisorByDivisorTester
                 continue;
             }
 
-            while (hasPrime && currentPrime < candidate)
-            {
-                hasPrime = primeEnumerator.MoveNext();
-                if (!hasPrime)
-                {
-                    break;
-                }
+			// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
+			// We don't need to additionally check it.
+			
+			// while (hasPrime && currentPrime < candidate)
+			// {
+			//     hasPrime = primeEnumerator.MoveNext();
+			//     if (!hasPrime)
+			//     {
+			//         break;
+			//     }
 
-                currentPrime = primeEnumerator.Current;
-            }
+			//     currentPrime = primeEnumerator.Current;
+			// }
 
-            if (!hasPrime || currentPrime != candidate)
-            {
-                markComposite();
-                printResult(candidate, false, false, false);
-                continue;
-            }
+			// if (!hasPrime || currentPrime != candidate)
+			// if (!Number.IsPrime(candidate))
+			// {
+			//     markComposite();
+			//     printResult(candidate, false, false, false);
+			//     continue;
+			// }
 
-            primesSpan[primeWriteIndex++] = candidate;
+			primesSpan[primeWriteIndex++] = candidate;
 
             if (candidate > maxPrime)
             {

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;
 

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<FactorEntry>.Shared.Rent(1);
+                local[0] = new FactorEntry(prime, 1);
+                return Rent(local, 1UL, true, 1);
+            }
+
+            FactorEntry[] extended = ArrayPool<FactorEntry>.Shared.Rent(Count + 1);
+            Array.Copy(source, 0, extended, 0, Count);
+            extended[Count] = new FactorEntry(prime, 1);
+            Span<FactorEntry> 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<FactorEntry>.Shared.Return(factors, clearArray: false);
+                }
+
+                _next = s_poolHead;
+                s_poolHead = this;
+            }
+        }
+    }
+}

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

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)

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)]

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<ulong>.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<ulong>.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)
     {