Merge pull request #130 from mlemanczyk/codex/optimize-number-types-for-performance

Add native-width Pollard-Brent variants

EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs

@@ -50,7 +50,7 @@ public class GpuPow2ModBenchmarks
         return Pow2ModSingleBit(Input.Exponent, Input.Modulus);
     }
 
-    private static GpuUInt128 Pow2ModSingleBit(ulong exponent, GpuUInt128 modulus)
+    private static GpuUInt128 Pow2ModSingleBit(ulong exponent, in GpuUInt128 modulus)
     {
         if (modulus.IsZero || modulus == GpuUInt128.One)
         {

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs

@@ -52,31 +52,14 @@ public class Pow2MontgomeryModBenchmarks
     /// Baseline right-to-left Montgomery ladder; measured 36.36 μs on the small sample set and 117.32 μs on the large one.
     /// </summary>
     [Benchmark(Baseline = true)]
-    public ulong BaselineCpu()
+    public ulong Baseline()
     {
         GetData(out var exponents, out var divisors, out _);
         ulong checksum = 0UL;
 
         for (int i = 0; i < SampleCount; i++)
         {
-            checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false);
+            checksum ^= exponents[i].Pow2MontgomeryModWindowed(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;
@@ -86,31 +69,14 @@ public class Pow2MontgomeryModBenchmarks
     /// Uses a known divisor cycle to fold the exponent first, dropping runtimes to 26.04 μs on small inputs (0.72× baseline) and 9.53 μs on large inputs (0.08× baseline).
     /// </summary>
     [Benchmark]
-    public ulong BaselineWithKnownCycleCpu()
+    public ulong BaselineWithKnownCycle()
     {
         GetData(out var exponents, out var divisors, out var cycles);
         ulong checksum = 0UL;
 
         for (int i = 0; i < SampleCount; i++)
         {
-            checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(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

@@ -89,32 +89,14 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
     /// large set.
     /// </summary>
     // [Benchmark(Baseline = true)]
-	public ulong MontgomeryWithoutCycleCpu()
+    public ulong MontgomeryWithoutCycle()
     {
         GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
         ulong checksum = 0UL;
 
         for (int i = 0; i < SampleCount; i++)
         {
-			checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false);
+            checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false);
-		}
-
-		return checksum;
-	}
-
-	/// <summary>
-	/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
-	/// large set.
-	/// </summary>
-	// [Benchmark(Baseline = true)]
-	public ulong MontgomeryWithoutCycleGpu()
-	{
-		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
-		ulong checksum = 0UL;
-
-		for (int i = 0; i < SampleCount; i++)
-		{
-			checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false);
         }
 
         return checksum;
@@ -125,32 +107,14 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
     /// on the large one (0.08×).
     /// </summary>
     [Benchmark]
-	public ulong MontgomeryWithPrecomputedCycleCpu()
+    public ulong MontgomeryWithPrecomputedCycle()
     {
         GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
         ulong checksum = 0UL;
 
         for (int i = 0; i < SampleCount; i++)
         {
-			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
-		}
-
-		return checksum;
-	}
-
-	/// <summary>
-	/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
-	/// on the large one (0.08×).
-	/// </summary>
-	[Benchmark]
-	public ulong MontgomeryWithPrecomputedCycleGpu()
-	{
-		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
-		ulong checksum = 0UL;
-
-		for (int i = 0; i < SampleCount; i++)
-		{
-			checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]);
         }
 
         return checksum;
@@ -161,7 +125,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
     /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
     /// </summary>
     // [Benchmark]
-	public ulong MontgomeryWithGpuCycleComputationCpu()
+    public ulong MontgomeryWithGpuCycleComputation()
     {
         GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
         ulong checksum = 0UL;
@@ -169,33 +133,15 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
         for (int i = 0; i < SampleCount; i++)
         {
             ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
-			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(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()
+    public ulong MontgomeryWithHeuristicCycleComputation()
     {
         GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
         ulong checksum = 0UL;
@@ -203,22 +149,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
         for (int i = 0; i < SampleCount; i++)
         {
             ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
-			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(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;
@@ -348,17 +279,16 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
 #if DEBUG
         Console.WriteLine("Trying heuristic. Prime order calculation");
 #endif
-		ulong order = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
             modulus,
             _previousPrimeOrder,
             divisorData,
-			PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-			PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-		if (order != 0UL)
+        if (orderResult.Order != 0UL)
         {
-			_previousPrimeOrder = order;
+            _previousPrimeOrder = orderResult.Order;
-			return order;
+            return orderResult.Order;
         }
 
         _previousPrimeOrder = null;
@@ -377,15 +307,15 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
 #if DEBUG
         Console.WriteLine("Trying heuristic. Wide prime order calculation");
 #endif
-		UInt128 order = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResultWide orderResult = PrimeOrderCalculator.Calculate(
             modulus,
             _previousWidePrimeOrder,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
 
-		if (order != UInt128.Zero)
+        if (orderResult.Order != UInt128.Zero)
         {
-			_previousWidePrimeOrder = order;
+            _previousWidePrimeOrder = orderResult.Order;
-			return order;
+            return orderResult.Order;
         }
 
         _previousWidePrimeOrder = null;

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs

@@ -176,7 +176,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
         return (UInt128)result;
     }
 
-    private static bool ShouldUseSingleBit(GpuUInt128 exponent)
+    private static bool ShouldUseSingleBit(in GpuUInt128 exponent)
     {
         if (exponent.High == 0UL && exponent.Low <= Pow2WindowFallbackThreshold)
         {
@@ -257,7 +257,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
         return result;
     }
 
-    private static bool IsBitSet(GpuUInt128 value, int bitIndex)
+    private static bool IsBitSet(in GpuUInt128 value, int bitIndex)
     {
         if (bitIndex >= 64)
         {
@@ -267,7 +267,7 @@ public class Pow2MontgomeryModUInt128Benchmarks
         return ((value.Low >> bitIndex) & 1UL) != 0UL;
     }
 
-    private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount)
+    private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount)
     {
         if (windowStart != 0)
         {

PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs

@@ -359,9 +359,10 @@ public class PrimeOrderGpuHeuristicsTests
             PrimeOrderGpuHeuristics.OverflowRegistry[prime] = 0;
 
             MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(prime);
-            ulong result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault, PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
+            PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(prime, null, divisorData, PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
 
-            result.Should().Be(3UL);
+            result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+            result.Order.Should().Be(3UL);
         }
         finally
         {

PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs

@@ -45,14 +45,14 @@ public class PrimeOrderCalculatorTests
     {
         ulong expected = ComputeOrderByDoubling(prime);
 
-        ulong result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Should().Be(expected);
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        result.Order.Should().Be(expected);
     }
 
     [Theory]
@@ -62,21 +62,21 @@ public class PrimeOrderCalculatorTests
     {
         ulong previousOrder = ComputeOrderByDoubling(previousPrime);
 
-        ulong heuristic = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        ulong strict = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        heuristic.Should().Be(strict);
+        heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        heuristic.Order.Should().Be(strict.Order);
     }
 
     [Theory]
@@ -85,14 +85,14 @@ public class PrimeOrderCalculatorTests
     [InlineData(3UL, 2UL)]
     public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder)
     {
-        ulong result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Should().Be(expectedOrder);
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        result.Order.Should().Be(expectedOrder);
     }
 
     [Fact]
@@ -105,14 +105,14 @@ public class PrimeOrderCalculatorTests
             maxPowChecks: 8,
             mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
 
-        ulong result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
             prime: 13UL,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(13UL),
-            config,
+            config);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Should().Be(ComputeOrderByDoubling(13UL));
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
+        result.Order.Should().Be(ComputeOrderByDoubling(13UL));
     }
 
     [Fact]
@@ -126,14 +126,14 @@ public class PrimeOrderCalculatorTests
             mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
 
         MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
-        ulong heuristic = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
             prime: 239UL,
             previousOrder: null,
             divisorData,
-            heuristicConfig,
+            heuristicConfig);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        heuristic.Should().Be(ComputeOrderByDoubling(239UL));
+        heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
+        heuristic.Order.Should().Be(ComputeOrderByDoubling(239UL));
 
         var strictConfig = new PrimeOrderCalculator.PrimeOrderSearchConfig(
             smallFactorLimit: 2,
@@ -141,14 +141,14 @@ public class PrimeOrderCalculatorTests
             maxPowChecks: 1,
             mode: PrimeOrderCalculator.PrimeOrderMode.Strict);
 
-        ulong strict = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
             prime: 239UL,
             previousOrder: null,
             divisorData,
-            strictConfig,
+            strictConfig);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        strict.Should().Be(heuristic);
+        strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        strict.Order.Should().Be(heuristic.Order);
     }
 
     [Fact]
@@ -157,12 +157,13 @@ public class PrimeOrderCalculatorTests
     {
         UInt128 prime = UInt128.Parse("18446744073709641691");
 
-        UInt128 result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResultWide result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
 
-        result.Should().Be((UInt128)1229782938247309446UL);
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        result.Order.Should().Be((UInt128)1229782938247309446UL);
     }
 
     [Fact]

PerfectNumbers.Core/Cpu/MersennePollardBrent.cs

@@ -0,0 +1,926 @@
+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

@@ -1,49 +0,0 @@
-using System;
-using System.Threading;
-
-namespace PerfectNumbers.Core;
-
-internal static class DeterministicRandom
-{
-    [ThreadStatic]
-    private static ulong s_threadState;
-
-    private static long s_seedCounter;
-
-    public static ulong NextUInt64()
-    {
-        ulong state = s_threadState;
-        if (state == 0UL)
-        {
-            state = InitializeThreadState();
-        }
-
-        state ^= state >> 12;
-        state ^= state << 25;
-        state ^= state >> 27;
-        s_threadState = state;
-        return state * 2685821657736338717UL;
-    }
-
-    public static UInt128 NextUInt128()
-    {
-        UInt128 high = (UInt128)NextUInt64();
-        UInt128 low = (UInt128)NextUInt64();
-        return (high << 64) | low;
-    }
-
-    private static ulong InitializeThreadState()
-    {
-        ulong seed = (ulong)Interlocked.Increment(ref s_seedCounter);
-        seed += 0x9E3779B97F4A7C15UL;
-        seed = (seed ^ (seed >> 30)) * 0xBF58476D1CE4E5B9UL;
-        seed = (seed ^ (seed >> 27)) * 0x94D049BB133111EBUL;
-        seed ^= seed >> 31;
-        if (seed == 0UL)
-        {
-            seed = 0x9E3779B97F4A7C15UL;
-        }
-
-        return seed;
-    }
-}

PerfectNumbers.Core/ExponentRemainderStepper.cs

@@ -1,20 +1,31 @@
 namespace PerfectNumbers.Core;
 
-internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor)
+internal struct ExponentRemainderStepper
 {
-    private readonly MontgomeryDivisorData _divisor = divisor;
+    private readonly MontgomeryDivisorData _divisor;
-    private readonly ulong _modulus = divisor.Modulus;
+    private readonly ulong _modulus;
-    private readonly ulong _nPrime = divisor.NPrime;
+    private readonly ulong _nPrime;
-    private readonly ulong _montgomeryOne = divisor.MontgomeryOne;
+    private readonly ulong _montgomeryOne;
-    private ulong _previousExponent = 0UL;
+    private ulong _previousExponent;
-    private ulong _currentMontgomery = divisor.MontgomeryOne;
+    private ulong _currentMontgomery;
-    private bool _hasState = false;
+    private bool _hasState;
 
-	public readonly bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL;
+    public ExponentRemainderStepper(in MontgomeryDivisorData divisor)
+    {
+        _divisor = divisor;
+        _modulus = divisor.Modulus;
+        _nPrime = divisor.NPrime;
+        _montgomeryOne = divisor.MontgomeryOne;
+        _previousExponent = 0UL;
+        _currentMontgomery = divisor.MontgomeryOne;
+        _hasState = false;
+    }
 
-    public readonly bool HasState => _hasState;
+    public bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL;
 
-    public readonly ulong PreviousExponent => _previousExponent;
+    public bool HasState => _hasState;
+
+    public ulong PreviousExponent => _previousExponent;
 
     public void Reset()
     {
@@ -36,7 +47,7 @@ internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor)
             // TODO: Route these state resets through the ProcessEightBitWindows helper once the scalar
             // Pow2MontgomeryMod implementation adopts it so delta stepping inherits the benchmarked
             // 2× gains for large exponents.
-            _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
             _previousExponent = exponent;
             _hasState = true;
             return ReduceCurrent();
@@ -48,7 +59,7 @@ internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor)
         // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
         // TODO: Replace this per-delta powmod with the upcoming windowed ladder so incremental
         // updates stop paying the single-bit cost highlighted in GpuPow2ModBenchmarks.
-        ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+        ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
         _previousExponent = exponent;
         return ReduceCurrent();
@@ -66,7 +77,7 @@ internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor)
         {
             // TODO: Switch this reload to the shared windowed pow2 helper once available so CPU
             // residue checks align with the optimized ProcessEightBitWindows timings.
-            _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
             _previousExponent = exponent;
             _hasState = true;
             return _currentMontgomery == _montgomeryOne;
@@ -77,7 +88,7 @@ internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor)
         // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor.
         // TODO: Use the windowed delta pow2 helper here as well to avoid the single-bit ladder that
         // currently lags behind the benchmarked implementation for huge divisor cycles.
-        ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+        ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
         _previousExponent = exponent;
         return _currentMontgomery == _montgomeryOne;
@@ -112,7 +123,7 @@ internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor)
         {
             // TODO: Replace this fallback path with the upcoming ProcessEightBitWindows helper so
             // fresh Montgomery states also benefit from the faster pow2 ladder measured on CPUs.
-            _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
             _previousExponent = exponent;
             _hasState = true;
             isUnity = _currentMontgomery == _montgomeryOne;

PerfectNumbers.Core/FactorEntry.cs

@@ -1,8 +0,0 @@
-namespace PerfectNumbers.Core
-{
-    public readonly struct FactorEntry(ulong value, int exponent)
-    {
-        public readonly ulong Value = value;
-        public readonly int Exponent = exponent;
-    }
-}

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(GpuUInt128 other) => AddInternal(other.High, other.Low);
+    public void Add(in GpuUInt128 other) => AddInternal(other.High, other.Low);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private void AddInternal(ulong otherHigh, ulong otherLow)
@@ -294,17 +294,17 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     public void AddMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void AddMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
+    public void AddMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => AddModInternal(value.High, value.Low, in modulus);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void AddMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus)
+    public void AddMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus)
     {
         ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
         AddModInternal(value.High, value.Low, in readOnlyMod);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void AddMod(GpuUInt128 value, GpuUInt128 modulus)
+    public void AddMod(in GpuUInt128 value, in GpuUInt128 modulus)
     {
         ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
         AddModInternal(value.High, value.Low, in readOnlyMod);
@@ -390,7 +390,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void AddMod(GpuUInt128 value, ulong modulus) => AddMod(value.Low, modulus);
+    public void AddMod(in GpuUInt128 value, ulong modulus) => AddMod(value.Low, modulus);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public void AddMod(ulong value, ulong modulus)
@@ -414,7 +414,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     public void Sub(in ReadOnlyGpuUInt128 other) => SubInternal(other.High, other.Low);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Sub(GpuUInt128 other) => SubInternal(other.High, other.Low);
+    public void Sub(in GpuUInt128 other) => SubInternal(other.High, other.Low);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private void SubInternal(ulong otherHigh, ulong otherLow)
@@ -461,7 +461,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     public void Mul64(in ReadOnlyGpuUInt128 other) => Mul64Internal(other.High, other.Low);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Mul64(GpuUInt128 other) => Mul64Internal(other.High, other.Low);
+    public void Mul64(in GpuUInt128 other) => Mul64Internal(other.High, other.Low);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private void Mul64Internal(ulong otherHigh, ulong otherLow)
@@ -588,7 +588,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static int TrailingZeroCount(GpuUInt128 value)
+    public static int TrailingZeroCount(in GpuUInt128 value)
     {
         if (value.Low != 0UL)
         {
@@ -646,17 +646,17 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     public void SubMod(in ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void SubMod(GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
+    public void SubMod(in GpuUInt128 value, in ReadOnlyGpuUInt128 modulus) => SubModInternal(value.High, value.Low, in modulus);
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void SubMod(in ReadOnlyGpuUInt128 value, GpuUInt128 modulus)
+    public void SubMod(in ReadOnlyGpuUInt128 value, in GpuUInt128 modulus)
     {
         ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
         SubModInternal(value.High, value.Low, in readOnlyMod);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void SubMod(GpuUInt128 value, GpuUInt128 modulus)
+    public void SubMod(in GpuUInt128 value, in GpuUInt128 modulus)
     {
         ReadOnlyGpuUInt128 readOnlyMod = modulus.AsReadOnly();
         SubModInternal(value.High, value.Low, in readOnlyMod);
@@ -727,7 +727,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Mul(GpuUInt128 other)
+    public void Mul(in GpuUInt128 other)
     {
         ReadOnlyGpuUInt128 readOnlyThis = AsReadOnly();
         ReadOnlyGpuUInt128 readOnlyOther = other.AsReadOnly();

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -551,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulong exponent = exponents[index];
-        hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = exponent.Pow2MontgomeryModWindowed(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.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
+        results[index] = exponent.Pow2MontgomeryModWindowed(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 ulong order)
+        out PrimeOrderCalculator.PrimeOrderResult result)
     {
-        order = 0UL;
+        result = default;
 
         var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
         var execution = lease.EnterExecutionScope();
@@ -497,17 +497,18 @@ internal static partial class PrimeOrderGpuHeuristics
             throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity.");
         }
 
+        ulong order = 0UL;
         resultBuffer.View.CopyToCPU(ref order, 1);
 
-        if (kernelStatus == PrimeOrderKernelStatus.FactoringFailure)
+        PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found
-        {
+            ? PrimeOrderCalculator.PrimeOrderStatus.Found
-            order = 0UL;
+            : PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved;
-        }
 
+        result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order);
         DisposeResources();
         lease.Dispose();
 
-        return order != 0UL;
+        return true;
 
         void DisposeResources()
         {
@@ -811,7 +812,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             ulong reduced = phi / factor;
-            if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL)
+            if (reduced.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL)
             {
                 return false;
             }
@@ -1595,7 +1596,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
     {
-        return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL;
     }
 
     private static ulong CalculateByDoublingKernel(ulong prime)
@@ -1795,7 +1796,7 @@ internal static partial class PrimeOrderGpuHeuristics
             return 0UL;
         }
 
-        return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
+        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
@@ -1811,7 +1812,7 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong exponent = exponents[index];
-        remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false);
+        remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)

PerfectNumbers.Core/HeuristicFailureLog.cs

@@ -1,98 +0,0 @@
-using System.Collections.Concurrent;
-using System.Globalization;
-using System.Text;
-
-namespace PerfectNumbers.Core
-{
-    internal static class HeuristicFailureLog
-    {
-        private const string LogFileName = "prime-order-heuristic-fallbacks.log";
-        private const int FlushThreshold = 4096;
-        private const int MaxBufferedEntries = 16;
-        private static readonly ThreadLocal<StringBuilder?> s_threadBuilder = new(() => StringBuilderPool.Rent(), trackAllValues: true);
-        private static readonly ThreadLocal<int> s_threadEntryCount = new(() => 0, trackAllValues: true);
-        private static readonly object s_fileLock = new();
-        private static readonly string s_logPath = Path.Combine(AppContext.BaseDirectory, LogFileName);
-        private static readonly UTF8Encoding s_utf8NoBom = new(false);
-
-        static HeuristicFailureLog()
-        {
-            AppDomain.CurrentDomain.ProcessExit += static (_, _) => FlushAll(force: true);
-        }
-
-        public static void Record(ulong prime, ulong? candidateOrder, HeuristicFailureReason reason)
-        {
-            string primeText = prime.ToString(CultureInfo.InvariantCulture);
-            string? candidateText = candidateOrder.HasValue ? candidateOrder.Value.ToString(CultureInfo.InvariantCulture) : null;
-            RecordInternal(primeText, candidateText, reason);
-        }
-
-        public static void Record(UInt128 prime, UInt128? candidateOrder, HeuristicFailureReason reason)
-        {
-            string primeText = prime.ToString(CultureInfo.InvariantCulture);
-            string? candidateText = candidateOrder.HasValue ? candidateOrder.Value.ToString(CultureInfo.InvariantCulture) : null;
-            RecordInternal(primeText, candidateText, reason);
-        }
-
-        private static void RecordInternal(string primeText, string? candidateText, HeuristicFailureReason reason)
-        {
-            StringBuilder? existing = s_threadBuilder.Value;
-            StringBuilder builder = existing ?? StringBuilderPool.Rent();
-            int count = s_threadEntryCount.Value + 1;
-            bool forceFlush = count >= MaxBufferedEntries;
-            builder.Append(DateTime.UtcNow.ToString("o", CultureInfo.InvariantCulture));
-            builder.Append(" | prime=");
-            builder.Append(primeText);
-            builder.Append(" | reason=");
-            builder.Append(reason);
-            if (candidateText is not null)
-            {
-                builder.Append(" | candidateOrder=");
-                builder.Append(candidateText);
-            }
-
-            builder.AppendLine();
-            FlushBuilderIfNeeded(builder, forceFlush);
-            if (builder.Length == 0)
-            {
-                count = 0;
-            }
-
-            s_threadEntryCount.Value = count;
-            s_threadBuilder.Value = builder;
-        }
-
-        private static void FlushAll(bool force)
-        {
-            foreach (StringBuilder? builder in s_threadBuilder.Values)
-            {
-                if (builder is not null)
-                {
-                    FlushBuilderIfNeeded(builder, force);
-                }
-            }
-        }
-
-        private static void FlushBuilderIfNeeded(StringBuilder builder, bool force)
-        {
-            if (builder.Length == 0)
-            {
-                return;
-            }
-
-            if (!force && builder.Length < FlushThreshold)
-            {
-                return;
-            }
-
-            string text = builder.ToString();
-            builder.Clear();
-            lock (s_fileLock)
-            {
-                using FileStream stream = new FileStream(s_logPath, FileMode.Append, FileAccess.Write, FileShare.ReadWrite);
-                using StreamWriter writer = new StreamWriter(stream, s_utf8NoBom);
-                writer.Write(text);
-            }
-        }
-    }
-}

PerfectNumbers.Core/HeuristicFailureReason.cs

@@ -1,11 +0,0 @@
-namespace PerfectNumbers.Core
-{
-    internal enum HeuristicFailureReason
-    {
-        PhiPartialFactorizationFailed,
-        StrictModeRequested,
-        HeuristicPipelineUnresolved,
-        HeuristicPipelineAssumedFullOrder,
-    }
-
-}

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -1,8 +1,11 @@
 using System.Buffers;
+using System.Buffers.Binary;
 using System.Collections.Generic;
 using System.Numerics;
 using System.Runtime.CompilerServices;
+using System.Security.Cryptography;
 using System.Text;
+using System.Threading;
 using PerfectNumbers.Core.Gpu;
 using ILGPU;
 using ILGPU.Runtime;
@@ -51,15 +54,15 @@ public class MersenneDivisorCycles
         return snapshot;
     }
 
-    public static bool CycleEqualsExponentCpu(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
+    public static bool CycleEqualsExponent(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
     {
         Dictionary<ulong, FactorCacheEntry>? cache = null;
-        return CycleEqualsExponentCpu(divisor, exponent, divisorData, ref cache);
+        return CycleEqualsExponent(divisor, exponent, divisorData, ref cache);
     }
 
-    public static bool CycleEqualsExponentCpu(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary<ulong, FactorCacheEntry>? factorCache)
+    public static bool CycleEqualsExponent(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary<ulong, FactorCacheEntry>? factorCache)
     {
-        if (CycleEqualsExponentForMersenneCandidateCpu(divisor, divisorData, exponent))
+        if (CycleEqualsExponentForMersenneCandidate(divisor, divisorData, exponent))
         {
             return true;
         }
@@ -94,15 +97,15 @@ public class MersenneDivisorCycles
             factorCache = cache;
         }
 
-        if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL)
+        if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL)
         {
             return cycleLength == exponent;
         }
 
-        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
     }
 
-    public static bool CycleEqualsExponentForMersenneCandidateCpu(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
+    public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
     {
         if (!IsValidMersenneDivisorCandidate(divisor, exponent))
         {
@@ -123,7 +126,7 @@ public class MersenneDivisorCycles
             }
         }
 
-        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
     }
 
     private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
@@ -231,14 +234,14 @@ public class MersenneDivisorCycles
             return (UInt128)cycle;
         }
 
-        UInt128 wideOrder = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
             PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-        if (wideOrder != UInt128.Zero)
+        if (wideResult.Order != UInt128.Zero)
         {
-            return wideOrder;
+            return wideResult.Order;
         }
 
         // Otherwise, find order of 2 mod divisor
@@ -261,7 +264,7 @@ public class MersenneDivisorCycles
         return order;
     }
 
-    public static bool TryCalculateCycleLengthForExponentCpu(
+    public static bool TryCalculateCycleLengthForExponent(
         ulong divisor,        
         ulong exponent,
         in MontgomeryDivisorData divisorData,
@@ -282,15 +285,15 @@ public class MersenneDivisorCycles
             return true;
         }
 
-        ulong order = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
             divisorData,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
             PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-        if (order != 0UL)
+        if (orderResult.Order != 0UL)
         {
-            cycleLength = order;
+            cycleLength = orderResult.Order;
             return true;
         }
 
@@ -324,7 +327,7 @@ public class MersenneDivisorCycles
             return false;
         }
 
-        cycleLength = ReduceOrderCpu(divisorData, phi, factorCounts);
+        cycleLength = ReduceOrder(divisorData, phi, factorCounts);
         return true;
     }
 
@@ -468,7 +471,7 @@ public class MersenneDivisorCycles
         }
     }
 
-    private static ulong ReduceOrderCpu(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
+    private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
     {
         if (factorCounts.Count == 0)
         {
@@ -502,7 +505,7 @@ public class MersenneDivisorCycles
                     }
 
                     ulong candidate = order / prime;
-                    if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL)
+                    if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
                     {
                         order = candidate;
                         continue;
@@ -527,31 +530,20 @@ public class MersenneDivisorCycles
             return 2UL;
         }
 
-        PollardRhoMontgomeryReducer reducer = PollardRhoMontgomeryReducer.Create(n);
-        ulong diff = 0UL;
-        ulong d = 0UL;
-
         while (true)
         {
-            ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+            ulong c = (NextRandomUInt64() % (n - 1UL)) + 1UL;
-            ulong xStandard = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
+            ulong x = (NextRandomUInt64() % (n - 2UL)) + 2UL;
-            ulong yStandard = xStandard;
+            ulong y = x;
-            ulong cMontgomery = reducer.ConvertToMontgomery(c);
+            ulong d = 1UL;
-            ulong xMontgomery = reducer.ConvertToMontgomery(xStandard);
-            ulong yMontgomery = xMontgomery;
-            d = 1UL;
 
             while (d == 1UL)
             {
-                xMontgomery = reducer.Advance(xMontgomery, cMontgomery);
+                x = AdvancePolynomial(x, c, n);
-                xStandard = reducer.AdvanceStandard(xStandard, c);
+                y = AdvancePolynomial(y, c, n);
+                y = AdvancePolynomial(y, c, n);
 
-                yMontgomery = reducer.Advance(yMontgomery, cMontgomery);
+                ulong diff = x > y ? x - y : y - x;
-                yStandard = reducer.AdvanceStandard(yStandard, c);
-                yMontgomery = reducer.Advance(yMontgomery, cMontgomery);
-                yStandard = reducer.AdvanceStandard(yStandard, c);
-
-                diff = xStandard > yStandard ? xStandard - yStandard : yStandard - xStandard;
                 d = BinaryGcd(diff, n);
             }
 
@@ -562,6 +554,13 @@ public class MersenneDivisorCycles
         }
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus)
+    {
+        UInt128 value = (UInt128)x * x + c;
+        return (ulong)(value % modulus);
+    }
+
     private static ulong BinaryGcd(ulong a, ulong b)
     {
         if (a == 0UL)
@@ -593,6 +592,13 @@ public class MersenneDivisorCycles
         }
     }
 
+    private static ulong NextRandomUInt64()
+    {
+        Span<byte> buffer = stackalloc byte[8];
+        RandomNumberGenerator.Fill(buffer);
+        return BinaryPrimitives.ReadUInt64LittleEndian(buffer);
+    }
+
     public readonly struct FactorCacheEntry
     {
         public FactorCacheEntry(ulong[] primes, byte[] exponents, int count)
@@ -934,15 +940,15 @@ public class MersenneDivisorCycles
 
         if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
         {
-            ulong order = PrimeOrderCalculator.Calculate(
+            PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
                     divisor,
                     previousOrder: null,
                     divisorData,
                     PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
                     PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-            if (order != 0UL)
+            if (orderResult.Order != 0UL)
             {
-                cycleLength = order;
+                cycleLength = orderResult.Order;
                 return true;
             }
         }

PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs

@@ -26,6 +26,8 @@ public static class MersenneNumberDivisorByDivisorTester
             return;
         }
 
+        _ = threadCount;
+
         candidates.Sort();
         Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
         int candidateCount = candidateSpan.Length;
@@ -36,9 +38,6 @@ public static class MersenneNumberDivisorByDivisorTester
         if (previousResults is not null && previousResults.Count > 0)
         {
             int recordedCount = previousResults.Count;
-            bool preferSortedMerge = recordedCount <= 4096 || (recordedCount * 4L) <= candidateCount;
-            if (preferSortedMerge)
-            {
             ulong[] recordedCandidates = ArrayPool<ulong>.Shared.Rent(recordedCount);
             try
             {
@@ -76,7 +75,7 @@ public static class MersenneNumberDivisorByDivisorTester
                 if (writeIndex < candidateCount)
                 {
                     candidates.RemoveRange(writeIndex, candidateCount - writeIndex);
-                        candidateSpan = CollectionsMarshal.AsSpan(candidates); // candidateSpan now tracks the reduced list.
+                    candidateSpan = CollectionsMarshal.AsSpan(candidates);
                     candidateCount = candidateSpan.Length;
                 }
             }
@@ -85,30 +84,6 @@ public static class MersenneNumberDivisorByDivisorTester
                 ArrayPool<ulong>.Shared.Return(recordedCandidates, clearArray: true);
             }
         }
-            else
-            {
-                int writeIndex = 0;
-                for (int readIndex = 0; readIndex < candidateCount; readIndex++)
-                {
-                    ulong candidate = candidateSpan[readIndex];
-
-                    if (previousResults.ContainsKey(candidate))
-                    {
-                        skippedByPreviousResults++;
-                        continue;
-                    }
-
-                    candidateSpan[writeIndex++] = candidate;
-                }
-
-                if (writeIndex < candidateCount)
-                {
-                    candidates.RemoveRange(writeIndex, candidateCount - writeIndex);
-                    candidateSpan = CollectionsMarshal.AsSpan(candidates); // candidateSpan now tracks the reduced list.
-                    candidateCount = candidateSpan.Length;
-                }
-            }
-        }
 
         if (skippedByPreviousResults > 0)
         {
@@ -288,32 +263,7 @@ public static class MersenneNumberDivisorByDivisorTester
                 MaxDegreeOfParallelism = workerCount
             };
 
-            primeCount = filteredPrimes.Count; // Reusing primeCount to track the filtered count.
+            Parallel.ForEach(filteredPrimes, options, ProcessPrime);
-            int primesPerWorker = (primeCount + workerCount - 1) / workerCount;
-
-            Parallel.For(
-                    0,
-                    workerCount,
-                    options,
-                    workerIndex =>
-                    {
-                        int start = workerIndex * primesPerWorker;
-                        if (start >= primeCount)
-                        {
-                            return;
-                        }
-
-                        int end = start + primesPerWorker;
-                        if (end > primeCount)
-                        {
-                            end = primeCount;
-                        }
-
-                        for (int i = start; i < end; i++)
-                        {
-                            ProcessPrime(filteredPrimes[i]);
-                        }
-                    });
         }
     }
 }

PerfectNumbers.Core/MersennePrimeFactorTester.cs

@@ -371,6 +371,10 @@ public static class MersennePrimeFactorTester
             return 2UL;
         }
 
+        var rng = new Random(1234567); // TODO: Replace this RNG with the deterministic span-based
+                                       // sequence measured faster in PollardRhoBenchmarks so we avoid
+                                       // per-call Random allocations.
+
         ulong c, d, x, y, diff;
         while (true)
         {
@@ -379,8 +383,8 @@ public static class MersennePrimeFactorTester
                 return 0UL;
             }
 
-            c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+            c = (ulong)rng.Next(2, int.MaxValue);
-            x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
+            x = (ulong)rng.Next(2, int.MaxValue);
             y = x;
             d = 1UL;
 
@@ -519,6 +523,9 @@ public static class MersennePrimeFactorTester
             return 2;
         }
 
+        var rng = new Random(1234567); // TODO: Replace with the deterministic UInt128 sequence from
+                                       // PollardRho128Benchmarks so we eliminate Random allocations.
+
         UInt128 c, d, x, y, diff;
         while (true)
         {
@@ -527,8 +534,8 @@ public static class MersennePrimeFactorTester
                 return 0;
             }
 
-            c = (DeterministicRandom.NextUInt128() % (n - 1)) + 1;
+            c = (UInt128)rng.NextInt64(2, long.MaxValue);
-            x = (DeterministicRandom.NextUInt128() % (n - 2)) + 2;
+            x = (UInt128)rng.NextInt64(2, long.MaxValue);
             y = x;
             d = 1;
 

PerfectNumbers.Core/PartialFactorResult.cs

@@ -1,76 +0,0 @@
-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();
-			}
-			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() { 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,7 +8,6 @@ public static class PerfectNumberConstants
     public const int MaxQForDivisorCycles = 4_000_000;
     public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000;
     public const int MaxOddPowersCount = 128;
-    public const int PooledArrayThreshold = 64;
     // TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal
     // divisor-cycle datasets we generate offline.
     // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune

PerfectNumbers.Core/PollardRhoMontgomeryReducer.cs

@@ -1,119 +0,0 @@
-using System;
-using System.Runtime.CompilerServices;
-
-namespace PerfectNumbers.Core;
-
-internal readonly struct PollardRhoMontgomeryReducer
-{
-    private readonly MontgomeryDivisorData _divisorData;
-    private readonly ulong _reciprocal;
-
-    [ThreadStatic]
-    private static bool s_hasCachedReciprocal;
-
-    [ThreadStatic]
-    private static ulong s_cachedReciprocal;
-
-    public PollardRhoMontgomeryReducer(in MontgomeryDivisorData divisorData, ulong reciprocal)
-    {
-        _divisorData = divisorData;
-        _reciprocal = reciprocal == 0UL || reciprocal == ulong.MaxValue
-            ? ReciprocalMath.ComputeReciprocalEstimate(divisorData.Modulus)
-            : reciprocal;
-    }
-
-    public ulong Modulus => _divisorData.Modulus;
-    private ulong NPrime => _divisorData.NPrime;
-    private ulong MontgomeryTwoSquared => _divisorData.MontgomeryTwoSquared;
-    private ulong MontgomeryOne => _divisorData.MontgomeryOne;
-    private ulong MontgomeryTwo => _divisorData.MontgomeryTwo;
-    public ulong Reciprocal => _reciprocal;
-
-    public static PollardRhoMontgomeryReducer Create(ulong modulus)
-    {
-        MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(modulus);
-        ulong reciprocal;
-        if (s_hasCachedReciprocal)
-        {
-            reciprocal = ReciprocalMath.RefineReciprocalEstimate(s_cachedReciprocal, modulus);
-            if (reciprocal == 0UL || reciprocal == ulong.MaxValue)
-            {
-                reciprocal = ReciprocalMath.ComputeReciprocalEstimate(modulus);
-            }
-        }
-        else
-        {
-            reciprocal = ReciprocalMath.ComputeReciprocalEstimate(modulus);
-        }
-
-        s_cachedReciprocal = reciprocal;
-        s_hasCachedReciprocal = true;
-        return new PollardRhoMontgomeryReducer(divisorData, reciprocal);
-    }
-
-    public static PollardRhoMontgomeryReducer Create(in MontgomeryDivisorData divisorData, ulong previousReciprocal)
-    {
-        ulong reciprocal = ReciprocalMath.RefineReciprocalEstimate(previousReciprocal, divisorData.Modulus);
-        if (reciprocal == 0UL || reciprocal == ulong.MaxValue)
-        {
-            reciprocal = ReciprocalMath.ComputeReciprocalEstimate(divisorData.Modulus);
-        }
-
-        s_cachedReciprocal = reciprocal;
-        s_hasCachedReciprocal = true;
-        return new PollardRhoMontgomeryReducer(divisorData, reciprocal);
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public ulong ConvertToMontgomery(ulong value)
-    {
-        if (value == 0UL)
-        {
-            return 0UL;
-        }
-
-        if (value == 1UL)
-        {
-            return MontgomeryOne;
-        }
-
-        if (value == 2UL)
-        {
-            return MontgomeryTwo;
-        }
-
-        return value.MontgomeryMultiply(MontgomeryTwoSquared, Modulus, NPrime);
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public ulong ReduceFromMontgomery(ulong montgomeryValue)
-    {
-        return montgomeryValue.MontgomeryMultiply(1UL, Modulus, NPrime);
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public ulong Advance(ulong montgomeryValue, ulong cMontgomery)
-    {
-        ulong squared = montgomeryValue.MontgomeryMultiply(montgomeryValue, Modulus, NPrime);
-        return MontgomeryAdd(squared, cMontgomery, Modulus);
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public ulong AdvanceStandard(ulong standardValue, ulong c)
-    {
-        UInt128 polynomial = (UInt128)standardValue * standardValue + c;
-        return ReciprocalMath.Reduce(polynomial, Modulus, _reciprocal);
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong MontgomeryAdd(ulong a, ulong b, ulong modulus)
-    {
-        UInt128 total = (UInt128)a + b;
-        if (total >= modulus)
-        {
-            total -= modulus;
-        }
-
-        return (ulong)total;
-    }
-}

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -9,18 +9,30 @@ namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
-    public static UInt128 Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
+    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 UInt128 Calculate(
+    public static PrimeOrderResultWide Calculate(
         UInt128 prime,
         UInt128? previousOrder,
         PrimeOrderSearchConfig config,
         PrimeOrderHeuristicDevice device)
     {
-        var scope = UsePow2Mode(device);
+        using var scope = UsePow2Mode(device);
         MontgomeryDivisorData divisorData;
-        UInt128 result;
         if (prime <= ulong.MaxValue)
         {
             ulong? previous = null;
@@ -39,49 +51,41 @@ internal static partial class PrimeOrderCalculator
 
             ulong prime64 = (ulong)prime;
             divisorData = MontgomeryDivisorData.FromModulus(prime64);
-            ulong order64 = Calculate(prime64, previous, divisorData, config, device);
+            PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device);
-            result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
+            return new PrimeOrderResultWide(result.Status, result.Order);
         }
         else
         {
             divisorData = default;
-            result = CalculateWideInternal(prime, previousOrder, divisorData, config);
         }
 
-        scope.Dispose();
+        return CalculateWideInternal(prime, previousOrder, divisorData, config);
-        return result;
     }
 
-    private static UInt128 CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+    private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
     {
         if (prime <= UInt128.One)
         {
-            return UInt128.One;
+            return new PrimeOrderResultWide(PrimeOrderStatus.Found, UInt128.One);
         }
 
         if (prime == (UInt128)3UL)
         {
-            return (UInt128)2UL;
+            return new PrimeOrderResultWide(PrimeOrderStatus.Found, (UInt128)2UL);
         }
 
         UInt128 phi = prime - UInt128.One;
-        var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
+        using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
         PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
-        UInt128 result;
         if (phiFactors.Factors is null)
         {
-            result = FinishStrictlyWide(prime, divisorData);
+            return FinishStrictlyWide(prime, divisorData, config.Mode);
-        }
-        else
-        {
-            result = RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
         }
 
-        debugScope.Dispose();
+        return RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
-        return result;
     }
 
-    private static UInt128 RunHeuristicPipelineWide(
+    private static PrimeOrderResultWide RunHeuristicPipelineWide(
         UInt128 prime,
         UInt128? previousOrder,
         in MontgomeryDivisorData divisorData,
@@ -91,7 +95,7 @@ internal static partial class PrimeOrderCalculator
     {
         if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
         {
-            return phi;
+            return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi);
         }
 
         UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData);
@@ -99,25 +103,27 @@ internal static partial class PrimeOrderCalculator
 
         if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config))
         {
-            return candidateOrder;
+            return new PrimeOrderResultWide(PrimeOrderStatus.Found, candidateOrder);
         }
 
         if (config.Mode == PrimeOrderMode.Strict)
         {
-            return FinishStrictlyWide(prime, divisorData);
+            return FinishStrictlyWide(prime, divisorData, PrimeOrderMode.Strict);
         }
 
         if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
         {
-            return order;
+            return new PrimeOrderResultWide(PrimeOrderStatus.Found, order);
         }
 
-        return candidateOrder;
+        return FinishStrictlyWide(prime, divisorData, config.Mode);
     }
 
-    private static UInt128 FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData)
+    private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData, PrimeOrderMode mode)
     {
-        return CalculateByFactorizationWide(prime, divisorData);
+        UInt128 strictOrder = CalculateByFactorizationWide(prime, divisorData);
+        PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved;
+        return new PrimeOrderResultWide(status, strictOrder);
     }
 
     private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -1,30 +1,65 @@
+using System;
 using System.Buffers;
 using System.Buffers.Binary;
+using System.Collections.Generic;
 using System.Diagnostics;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+using System.Security.Cryptography;
 using PerfectNumbers.Core.Gpu;
+using System.Threading;
 
 namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
+    internal enum PrimeOrderStatus
+    {
+        Found,
+        HeuristicUnresolved,
+    }
+
+    internal readonly struct PrimeOrderResult
+    {
+        public PrimeOrderResult(PrimeOrderStatus status, ulong order)
+        {
+            Status = status;
+            Order = order;
+        }
+
+        public PrimeOrderStatus Status { get; }
+
+        public ulong Order { get; }
+    }
+
     internal enum PrimeOrderMode
     {
         Heuristic,
         Strict,
     }
 
-    internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode)
+    internal readonly struct PrimeOrderSearchConfig
     {
-		public static PrimeOrderSearchConfig HeuristicDefault => new(smallFactorLimit: 100_000, pollardRhoMilliseconds: 128, maxPowChecks: 24, PrimeOrderMode.Heuristic);
+        public PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode)
-        public static PrimeOrderSearchConfig StrictDefault => new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict);
+        {
+            SmallFactorLimit = smallFactorLimit;
+            PollardRhoMilliseconds = pollardRhoMilliseconds;
+            MaxPowChecks = maxPowChecks;
+            Mode = mode;
+        }
 
-		public readonly uint SmallFactorLimit = smallFactorLimit;
+        public uint SmallFactorLimit { get; }
-		public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds;
+
-		public readonly int MaxPowChecks = maxPowChecks;
+        public int PollardRhoMilliseconds { get; }
-		public readonly PrimeOrderMode Mode = mode;		
+
+        public int MaxPowChecks { get; }
+
+        public PrimeOrderMode Mode { get; }
+
+        public static PrimeOrderSearchConfig HeuristicDefault => new(100_000, 128, 24, PrimeOrderMode.Heuristic);
+
+        public static PrimeOrderSearchConfig StrictDefault => new(1_000_000, 0, 0, PrimeOrderMode.Strict);
     }
 
     internal enum PrimeOrderHeuristicDevice
@@ -48,18 +83,8 @@ internal static partial class PrimeOrderCalculator
     [ThreadStatic]
     private static bool s_debugLoggingEnabled;
 
-    [ThreadStatic]
-    private static List<ulong>? s_partialFactorPending;
-
-    [ThreadStatic]
-    private static Stack<ulong>? s_partialFactorCompositeStack;
-
-    [ThreadStatic]
-    private static Dictionary<ulong, int>? s_partialFactorCounts;
-
     private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu;
-
+    private readonly struct Pow2ModeScope : IDisposable
-	private readonly struct Pow2ModeScope
     {
         private readonly bool _previousInitialized;
         private readonly bool _previousAllow;
@@ -83,7 +108,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private readonly struct DebugLoggingScope
+    private readonly struct DebugLoggingScope : IDisposable
     {
         private readonly bool _previous;
 
@@ -121,95 +146,93 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    public static ulong Calculate(
+    public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
-        ulong prime,
+        => Calculate(prime, previousOrder, divisorData, config, PrimeOrderHeuristicDevice.Gpu);
-        ulong? previousOrder,
-        in MontgomeryDivisorData divisorData,
-        in PrimeOrderSearchConfig config)
-        => Calculate(prime, previousOrder, divisorData, config, PrimeOrderHeuristicDevice.Cpu);
 
-    public static ulong Calculate(
+    public static PrimeOrderResult Calculate(
         ulong prime,
         ulong? previousOrder,
         in MontgomeryDivisorData divisorData,
-        in PrimeOrderSearchConfig config,
+        PrimeOrderSearchConfig config,
         PrimeOrderHeuristicDevice device)
     {
-        var scope = UsePow2Mode(device);
+        using var scope = UsePow2Mode(device);
-        ulong order = CalculateInternal(prime, previousOrder, divisorData, config);
+        return CalculateInternal(prime, previousOrder, divisorData, config);
-        scope.Dispose();
-        return order;
     }
 
-    private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
+    private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
     {
         if (prime <= 3UL)
         {
-            return prime == 3UL ? 2UL : 1UL;
+            return new PrimeOrderResult(PrimeOrderStatus.Found, prime == 3UL ? 2UL : 1UL);
         }
 
         ulong phi = prime - 1UL;
 
-        if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder))
+        if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out PrimeOrderResult gpuResult))
         {
-            return gpuOrder;
+            return gpuResult;
         }
 
+        using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
+        DebugLog("Partial factoring φ(p)");
+
         PartialFactorResult phiFactors = PartialFactor(phi, config);
 
-		ulong result;
         if (phiFactors.Factors is null)
         {
-			result = FinishStrictlyCpu(prime, divisorData);
+            DebugLog("No factors found");
-			phiFactors.Dispose();
+            return FinishStrictly(prime, config.Mode);
-			return result;
         }
 
-		result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors);
+        return RunHeuristicPipeline(prime, previousOrder, config, divisorData, phi, phiFactors);
-		phiFactors.Dispose();
-		return result;
     }
 
-    private static ulong RunHeuristicPipelineCpu(
+    private static PrimeOrderResult RunHeuristicPipeline(
         ulong prime,
         ulong? previousOrder,
-        in PrimeOrderSearchConfig config,
+        PrimeOrderSearchConfig config,
         in MontgomeryDivisorData divisorData,
         ulong phi,
         PartialFactorResult phiFactors)
     {
-        if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData))
+        DebugLog("Trying special max check");
+        if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData))
         {
-            return phi;
+            return new PrimeOrderResult(PrimeOrderStatus.Found, phi);
         }
 
-        ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData);
+        DebugLog("Initializing starting order");
-        candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData);
+        ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData);
+        candidateOrder = ExponentLowering(candidateOrder, prime, phiFactors, divisorData);
 
-        if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config))
+        DebugLog("Trying to confirm order");
+        if (TryConfirmOrder(prime, candidateOrder, divisorData, config))
         {
-            return candidateOrder;
+            return new PrimeOrderResult(PrimeOrderStatus.Found, candidateOrder);
         }
 
         if (config.Mode == PrimeOrderMode.Strict)
         {
-            return FinishStrictlyCpu(prime, divisorData);
+            return FinishStrictly(prime, PrimeOrderMode.Strict);
         }
 
-        if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order))
+        if (TryHeuristicFinish(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order))
         {
-            return order;
+            return new PrimeOrderResult(PrimeOrderStatus.Found, order);
         }
 
-        return candidateOrder;
+        DebugLog("Heuristic unresolved, finishing strictly");
+        return FinishStrictly(prime, config.Mode);
     }
 
-    private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData)
+    private static PrimeOrderResult FinishStrictly(ulong prime, PrimeOrderMode mode)
     {
-        return CalculateByFactorizationCpu(prime, divisorData);
+        ulong strictOrder = CalculateByDoubling(prime);
+        return new PrimeOrderResult(mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved, strictOrder);
     }
 
-    private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
+    private static bool TrySpecialMax(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
     {
         ReadOnlySpan<FactorEntry> factorSpan = factors.Factors;
         int length = factors.Count;
@@ -217,7 +240,7 @@ internal static partial class PrimeOrderCalculator
         {
             ulong factor = factorSpan[i].Value;
             ulong reduced = phi / factor;
-            if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+            if (Pow2EqualsOne(reduced, prime, divisorData))
             {
                 return false;
             }
@@ -226,13 +249,13 @@ internal static partial class PrimeOrderCalculator
         return true;
     }
 
-    private static ulong InitializeStartingOrderCpu(ulong prime, ulong phi, in MontgomeryDivisorData divisorData)
+    private static ulong InitializeStartingOrder(ulong prime, ulong phi, in MontgomeryDivisorData divisorData)
     {
         ulong order = phi;
         if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL)
         {
             ulong half = phi >> 1;
-            if (Pow2EqualsOneCpu(half, prime, divisorData))
+            if (Pow2EqualsOne(half, prime, divisorData))
             {
                 order = half;
             }
@@ -241,7 +264,7 @@ internal static partial class PrimeOrderCalculator
         return order;
     }
 
-    private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
+    private static ulong ExponentLowering(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
     {
         FactorEntry[]? tempArray = null;
         try
@@ -274,7 +297,7 @@ internal static partial class PrimeOrderCalculator
                     }
 
                     ulong reduced = order / primeFactor;
-                    if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+                    if (Pow2EqualsOne(reduced, prime, divisorData))
                     {
                         order = reduced;
                         continue;
@@ -295,7 +318,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
+    private static bool TryConfirmOrder(ulong prime, ulong order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
     {
         if (order == 0UL)
         {
@@ -303,18 +326,16 @@ internal static partial class PrimeOrderCalculator
         }
 
         // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`.
-        // DebugLog("Verifying a^order ≡ 1 (mod p)");
+        DebugLog("Verifying a^order ≡ 1 (mod p)");
-        if (!Pow2EqualsOneCpu(order, prime, divisorData))
+        if (!Pow2EqualsOne(order, prime, divisorData))
         {
             return false;
         }
 
-        // DebugLog("Partial factoring order");
+        DebugLog("Partial factoring order");
 
         // TODO: Do we do partial factoring of order multiple times?
         PartialFactorResult factorization = PartialFactor(order, config);
-        try
-        {
         if (factorization.Factors is null)
         {
             return false;
@@ -324,10 +345,12 @@ internal static partial class PrimeOrderCalculator
         {
             if (factorization.Cofactor <= 1UL)
             {
+            DebugLog("Cofactor <= 1. No factors found");
                 return false;
             }
 
-                // DebugLog("Cofactor > 1. Testing primality of cofactor");
+            // TODO: Use Open.Numerics.Primality for this final check once it's available.
+            DebugLog("Cofactor > 1. Testing primality of cofactor");
             bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
             // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
             if (!isPrime)
@@ -335,14 +358,12 @@ internal static partial class PrimeOrderCalculator
                 return false;
             }
 
-                // DebugLog("Adding cofactor as prime factor");
+            DebugLog("Adding cofactor as prime factor");
-                PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
+            factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
-                factorization.Dispose();
-                factorization = extended;
         }
 
         ReadOnlySpan<FactorEntry> span = factorization.Factors;
-            // DebugLog("Verifying prime-power reductions");
+        DebugLog("Verifying prime-power reductions");
         int length = factorization.Count;
         for (int i = 0; i < length; i++)
         {
@@ -356,7 +377,7 @@ internal static partial class PrimeOrderCalculator
                 }
 
                 reduced /= primeFactor;
-                    if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+                if (Pow2EqualsOne(reduced, prime, divisorData))
                 {
                     return false;
                 }
@@ -365,18 +386,13 @@ internal static partial class PrimeOrderCalculator
 
         return true;
     }
-        finally
-        {
-            factorization.Dispose();
-        }
-    }
 
-    private static bool TryHeuristicFinishCpu(
+    private static bool TryHeuristicFinish(
         ulong prime,
         ulong order,
         ulong? previousOrder,
         in MontgomeryDivisorData divisorData,
-        in PrimeOrderSearchConfig config,
+        PrimeOrderSearchConfig config,
         PartialFactorResult phiFactors,
         out ulong result)
     {
@@ -387,10 +403,8 @@ internal static partial class PrimeOrderCalculator
         }
 
         // TODO: Do we do partial factoring of order multiple times?
-        // DebugLog("Trying heuristic. Partial factoring order");
+        DebugLog("Trying heuristic. Partial factoring order");
         PartialFactorResult orderFactors = PartialFactor(order, config);
-        try
-        {
         if (orderFactors.Factors is null)
         {
             return false;
@@ -410,22 +424,20 @@ internal static partial class PrimeOrderCalculator
                 return false;
             }
 
-                PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
+            orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
-                orderFactors.Dispose();
-                orderFactors = extended;
         }
 
         int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4;
         List<ulong> candidates = new(capacity);
         FactorEntry[] factorArray = orderFactors.Factors!;
-            // DebugLog("Building candidates list");
+        DebugLog("Building candidates list");
         BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
         if (candidates.Count == 0)
         {
             return false;
         }
 
-            // DebugLog("Sorting candidates");
+        DebugLog("Sorting candidates");
         SortCandidates(prime, previousOrder, candidates);
 
         int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
@@ -434,7 +446,7 @@ internal static partial class PrimeOrderCalculator
         bool allowGpuBatch = true;
         Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
 
-            // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
+        DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
         int index = 0;
         const int MaxGpuBatchSize = 1024;
         const int StackGpuBatchSize = 1024;
@@ -504,14 +516,14 @@ internal static partial class PrimeOrderCalculator
                 }
                 else
                 {
-                        equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData);
+                    equalsOne = Pow2EqualsOne(candidate, prime, divisorData);
                 }
                 if (!equalsOne)
                 {
                     continue;
                 }
 
-                    if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget))
+                if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget))
                 {
                     continue;
                 }
@@ -533,14 +545,9 @@ internal static partial class PrimeOrderCalculator
             index += batchSize;
         }
 
-            // DebugLog("No candidate confirmed");
+        DebugLog("No candidate confirmed");
         return false;
     }
-        finally
-        {
-            orderFactors.Dispose();
-        }
-    }
 
     private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates)
     {
@@ -709,11 +716,9 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget)
+    private static bool TryConfirmCandidate(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, ref int powUsed, int powBudget)
     {
         PartialFactorResult factorization = PartialFactor(candidate, config);
-        try
-        {
         if (factorization.Factors is null)
         {
             return false;
@@ -733,9 +738,7 @@ internal static partial class PrimeOrderCalculator
                 return false;
             }
 
-                PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
+            factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
-                factorization.Dispose();
-                factorization = extended;
         }
 
         ReadOnlySpan<FactorEntry> span = factorization.Factors;
@@ -758,7 +761,7 @@ internal static partial class PrimeOrderCalculator
                 }
 
                 powUsed++;
-                    if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+                if (Pow2EqualsOne(reduced, prime, divisorData))
                 {
                     return false;
                 }
@@ -767,13 +770,8 @@ internal static partial class PrimeOrderCalculator
 
         return true;
     }
-        finally
-        {
-            factorization.Dispose();
-        }
-    }
 
-    private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
+    private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
     {
         if (IsGpuPow2Allowed)
         {
@@ -784,59 +782,10 @@ internal static partial class PrimeOrderCalculator
             }
         }
 
-        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
     }
 
-    private static List<ulong> AcquirePendingCompositeList(int capacityHint)
+    private static PartialFactorResult PartialFactor(ulong value, PrimeOrderSearchConfig config)
-    {
-        List<ulong>? list = s_partialFactorPending;
-        if (list is null)
-        {
-            list = new List<ulong>(capacityHint);
-            s_partialFactorPending = list;
-            return list;
-        }
-
-        list.Clear();
-        if (list.Capacity < capacityHint)
-        {
-            list.Capacity = capacityHint;
-        }
-
-        return list;
-    }
-
-    private static Stack<ulong> AcquireCompositeStack(int capacityHint)
-    {
-        Stack<ulong>? stack = s_partialFactorCompositeStack;
-        if (stack is null)
-        {
-            stack = new Stack<ulong>(capacityHint);
-            s_partialFactorCompositeStack = stack;
-            return stack;
-        }
-
-        stack.Clear();
-        stack.EnsureCapacity(capacityHint);
-        return stack;
-    }
-
-    private static Dictionary<ulong, int> AcquireFactorCountDictionary(int capacityHint)
-    {
-        Dictionary<ulong, int>? dictionary = s_partialFactorCounts;
-        if (dictionary is null)
-        {
-            dictionary = new Dictionary<ulong, int>(capacityHint);
-            s_partialFactorCounts = dictionary;
-            return dictionary;
-        }
-
-        dictionary.Clear();
-        dictionary.EnsureCapacity(capacityHint);
-        return dictionary;
-    }
-
-    private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config)
     {
         if (value <= 1UL)
         {
@@ -854,19 +803,18 @@ internal static partial class PrimeOrderCalculator
         bool useDictionary = false;
 
         uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
-        ulong remaining = value;
+        ulong remaining;
 
         bool gpuFactored = false;
         if (IsGpuHeuristicDevice)
         {
             gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _);
         }
-
+        else
-        List<ulong> pending = AcquirePendingCompositeList(2);
-        Stack<ulong>? compositeStack = null;
-
-        try
         {
+            remaining = value;
+        }
+
         if (!gpuFactored)
         {
             primeSlots.Clear();
@@ -874,11 +822,12 @@ internal static partial class PrimeOrderCalculator
             if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining))
             {
                 useDictionary = true;
-                    counts = AcquireFactorCountDictionary(Math.Max(factorCount, 8));
+                counts = new Dictionary<ulong, int>(capacity: 8);
                 remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
             }
         }
 
+        List<ulong> pending = new();
         if (remaining > 1UL)
         {
             pending.Add(remaining);
@@ -886,15 +835,16 @@ internal static partial class PrimeOrderCalculator
 
         if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
         {
+            DebugLog("Processing pending composites with Pollard's Rho");
             Stopwatch stopwatch = Stopwatch.StartNew();
             long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks;
-                compositeStack = AcquireCompositeStack(Math.Max(pending.Count * 2, 4));
+            Stack<ulong> stack = new();
-                compositeStack.Push(remaining);
+            stack.Push(remaining);
             pending.Clear();
 
-                while (compositeStack.Count > 0)
+            while (stack.Count > 0)
             {
-                    ulong composite = compositeStack.Pop();
+                ulong composite = stack.Pop();
                 if (composite <= 1UL)
                 {
                     continue;
@@ -920,8 +870,8 @@ internal static partial class PrimeOrderCalculator
                 }
 
                 ulong quotient = composite / factor;
-                    compositeStack.Push(factor);
+                stack.Push(factor);
-                    compositeStack.Push(quotient);
+                stack.Push(quotient);
             }
         }
 
@@ -945,16 +895,18 @@ internal static partial class PrimeOrderCalculator
         {
             if ((counts is null || counts.Count == 0) && cofactor == value)
             {
-                    return PartialFactorResult.Rent(null, value, false, 0);
+                DebugLog("cofactor is the same as value, no factors found");
+                return new PartialFactorResult(null, value, false, 0);
             }
 
             if (counts is null || counts.Count == 0)
             {
-                    return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
+                return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
             }
 
             FactorEntry[] factors = ArrayPool<FactorEntry>.Shared.Rent(counts.Count);
             int index = 0;
+            DebugLog(() => $"Collecting {counts.Count} prime factors");
             foreach (KeyValuePair<ulong, int> entry in counts)
             {
                 factors[index] = new FactorEntry(entry.Key, entry.Value);
@@ -964,8 +916,12 @@ internal static partial class PrimeOrderCalculator
             Span<FactorEntry> span = factors.AsSpan(0, index);
             span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
 
+            FactorEntry[] resultArray = new FactorEntry[index];
+            span.CopyTo(resultArray);
+            ArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
+
             bool fullyFactored = cofactor == 1UL;
-                return PartialFactorResult.Rent(factors, cofactor, fullyFactored, index);
+            return new PartialFactorResult(resultArray, cofactor, fullyFactored, index);
         }
 
         int actualCount = 0;
@@ -981,14 +937,16 @@ internal static partial class PrimeOrderCalculator
         {
             if (cofactor == value)
             {
-                    return PartialFactorResult.Rent(null, value, false, 0);
+                DebugLog("cofactor is the same as value, no factors found");
+                return new PartialFactorResult(null, value, false, 0);
             }
 
-                return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
+            return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
         }
 
-            FactorEntry[] array = ArrayPool<FactorEntry>.Shared.Rent(actualCount);
+        FactorEntry[] array = new FactorEntry[actualCount];
         int arrayIndex = 0;
+        DebugLog(() => $"Collecting {actualCount} prime factors");
         for (int i = 0; i < factorCount; i++)
         {
             ulong primeValue = primeSlots[i];
@@ -1002,17 +960,9 @@ internal static partial class PrimeOrderCalculator
             arrayIndex++;
         }
 
-            Span<FactorEntry> arraySpan = array.AsSpan(0, arrayIndex);
+        Array.Sort(array, static (a, b) => a.Value.CompareTo(b.Value));
-            arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value));
         bool fullyFactoredArray = cofactor == 1UL;
-            return PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex);
+        return new PartialFactorResult(array, cofactor, fullyFactoredArray, array.Length);
-        }
-        finally
-        {
-            pending.Clear();
-            counts?.Clear();
-            compositeStack?.Clear();
-        }
     }
 
     private const int GpuSmallPrimeFactorSlots = 64;
@@ -1164,9 +1114,7 @@ internal static partial class PrimeOrderCalculator
             return true;
         }
 
-        PollardRhoMontgomeryReducer reducer = PollardRhoMontgomeryReducer.Create(n);
+        Span<byte> buffer = stackalloc byte[8];
-        ulong diff = 0UL;
-        ulong d = 0UL;
 
         while (true)
         {
@@ -1175,13 +1123,12 @@ internal static partial class PrimeOrderCalculator
                 return false;
             }
 
-            ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+            RandomNumberGenerator.Fill(buffer);
-            ulong xStandard = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
+            ulong c = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 1UL)) + 1UL;
-            ulong yStandard = xStandard;
+            RandomNumberGenerator.Fill(buffer);
-            ulong cMontgomery = reducer.ConvertToMontgomery(c);
+            ulong x = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 2UL)) + 2UL;
-            ulong xMontgomery = reducer.ConvertToMontgomery(xStandard);
+            ulong y = x;
-            ulong yMontgomery = xMontgomery;
+            ulong d = 1UL;
-            d = 1UL;
 
             while (d == 1UL)
             {
@@ -1190,15 +1137,10 @@ internal static partial class PrimeOrderCalculator
                     return false;
                 }
 
-                xMontgomery = reducer.Advance(xMontgomery, cMontgomery);
+                x = AdvancePolynomial(x, c, n);
-                xStandard = reducer.AdvanceStandard(xStandard, c);
+                y = AdvancePolynomial(y, c, n);
-
+                y = AdvancePolynomial(y, c, n);
-                yMontgomery = reducer.Advance(yMontgomery, cMontgomery);
+                ulong diff = x > y ? x - y : y - x;
-                yStandard = reducer.AdvanceStandard(yStandard, c);
-                yMontgomery = reducer.Advance(yMontgomery, cMontgomery);
-                yStandard = reducer.AdvanceStandard(yStandard, c);
-
-                diff = xStandard > yStandard ? xStandard - yStandard : yStandard - xStandard;
                 d = BinaryGcd(diff, n);
             }
 
@@ -1210,6 +1152,13 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus)
+    {
+        UInt128 value = (UInt128)x * x + c;
+        return (ulong)(value % modulus);
+    }
+
     private static ulong BinaryGcd(ulong a, ulong b)
     {
         if (a == 0UL)
@@ -1305,7 +1254,7 @@ internal static partial class PrimeOrderCalculator
             return;
         }
 
-        counts ??= AcquireFactorCountDictionary(Math.Max(count, 8));
+        counts ??= new Dictionary<ulong, int>(Math.Max(count, 8));
         CopyFactorsToDictionary(primes, exponents, count, counts);
         count = 0;
         useDictionary = true;
@@ -1331,110 +1280,59 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData)
+    private static ulong CalculateByDoubling(ulong prime)
     {
-        ulong phi = prime - 1UL;
+        ulong order = 1UL;
-        Dictionary<ulong, int> counts = new(capacity: 8);
+        ulong pow = 2UL;
-        FactorCompletely(phi, counts);
+
-        if (counts.Count == 0)
+        while (pow != 1UL)
         {
-            return phi;
+            pow <<= 1;
+            if (pow >= prime)
+            {
+                pow -= prime;
             }
 
-        List<KeyValuePair<ulong, int>> entries = new(counts);
+            order++;
-        entries.Sort(static (a, b) => a.Key.CompareTo(b.Key));
-
-        ulong order = phi;
-        int entryCount = entries.Count;
-        for (int i = 0; i < entryCount; i++)
-        {
-            ulong primeFactor = entries[i].Key;
-            int exponent = entries[i].Value;
-            for (int iteration = 0; iteration < exponent; iteration++)
-            {
-                if ((order % primeFactor) != 0UL)
-                {
-                    break;
-                }
-
-                ulong candidate = order / primeFactor;
-                if (Pow2EqualsOneCpu(candidate, prime, divisorData))
-                {
-                    order = candidate;
-                    continue;
-                }
-
-                break;
-            }
         }
 
         return order;
     }
 
-    private static void FactorCompletely(ulong value, Dictionary<ulong, int> counts)
-    {
-        if (value <= 1UL)
-        {
-            return;
-        }
-
-        if (Open.Numeric.Primes.Prime.Numbers.IsPrime(value))
-        {
-            AddFactor(counts, value, 1);
-            return;
-        }
-
-        ulong factor = PollardRhoStrict(value);
-        FactorCompletely(factor, counts);
-        FactorCompletely(value / factor, counts);
-    }
-
-    private static ulong PollardRhoStrict(ulong n)
-    {
-        if ((n & 1UL) == 0UL)
-        {
-            return 2UL;
-        }
-
-        PollardRhoMontgomeryReducer reducer = PollardRhoMontgomeryReducer.Create(n);
-        ulong diff = 0UL;
-        ulong d = 0UL;
-
-        while (true)
-        {
-            ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
-            ulong xStandard = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
-            ulong yStandard = xStandard;
-            ulong cMontgomery = reducer.ConvertToMontgomery(c);
-            ulong xMontgomery = reducer.ConvertToMontgomery(xStandard);
-            ulong yMontgomery = xMontgomery;
-            d = 1UL;
-
-            while (d == 1UL)
-            {
-                xMontgomery = reducer.Advance(xMontgomery, cMontgomery);
-                xStandard = reducer.AdvanceStandard(xStandard, c);
-
-                yMontgomery = reducer.Advance(yMontgomery, cMontgomery);
-                yStandard = reducer.AdvanceStandard(yStandard, c);
-                yMontgomery = reducer.Advance(yMontgomery, cMontgomery);
-                yStandard = reducer.AdvanceStandard(yStandard, c);
-
-                diff = xStandard > yStandard ? xStandard - yStandard : yStandard - xStandard;
-                d = BinaryGcd(diff, n);
-            }
-
-            if (d != n)
-            {
-                return d;
-            }
-        }
-    }
-
     private readonly struct CandidateKey(int primary, long secondary, long tertiary)
     {
         public readonly int Primary = primary;
         public readonly long Secondary = secondary;
         public readonly long Tertiary = tertiary;
     }
+
+    private readonly struct FactorEntry(ulong value, int exponent)
+    {
+        public readonly ulong Value = value;
+        public readonly int Exponent = exponent;
+    }
+
+    private readonly struct PartialFactorResult(FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count)
+    {
+        public readonly FactorEntry[]? Factors = factors;
+        public readonly ulong Cofactor = cofactor;
+        public readonly bool FullyFactored = fullyFactored;
+        public readonly int Count = count;
+        public readonly static PartialFactorResult Empty = new(null, 1UL, true, 0);
+
+        public PartialFactorResult WithAdditionalPrime(ulong prime)
+        {
+            if (Factors is null)
+            {
+                FactorEntry[] local = [new FactorEntry(prime, 1)];
+                return new PartialFactorResult(local, 1UL, true, 1);
+            }
+
+            FactorEntry[] extended = new FactorEntry[Count + 1];
+            Array.Copy(Factors, 0, extended, 0, Count);
+            extended[Count] = new FactorEntry(prime, 1);
+            Array.Sort(extended, static (a, b) => a.Value.CompareTo(b.Value));
+            return new PartialFactorResult(extended, 1UL, true, Count + 1);
+        }
+    }
 }

PerfectNumbers.Core/ReciprocalMath.cs

@@ -1,125 +0,0 @@
-using System.Numerics;
-using System.Runtime.CompilerServices;
-
-namespace PerfectNumbers.Core;
-
-internal static class ReciprocalMath
-{
-    private const int ReciprocalShift = 64;
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong ComputeReciprocalEstimate(ulong divisor)
-    {
-        if (divisor <= 1UL)
-        {
-            return ulong.MaxValue;
-        }
-
-        UInt128 scaledOne = UInt128.One << ReciprocalShift;
-        UInt128 estimate = scaledOne / divisor;
-        if ((scaledOne % divisor) != UInt128.Zero)
-        {
-            estimate += UInt128.One;
-        }
-
-        if (estimate == UInt128.Zero)
-        {
-            return 1UL;
-        }
-
-        if (estimate > ulong.MaxValue)
-        {
-            return ulong.MaxValue;
-        }
-
-        return (ulong)estimate;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong RefineReciprocalEstimate(ulong reciprocal, ulong divisor)
-    {
-        if (divisor <= 1UL)
-        {
-            return ulong.MaxValue;
-        }
-
-        if (reciprocal == 0UL || reciprocal == ulong.MaxValue)
-        {
-            return ComputeReciprocalEstimate(divisor);
-        }
-
-        UInt128 reciprocal128 = reciprocal;
-        UInt128 product = (UInt128)divisor * reciprocal128;
-        UInt128 scaled = product >> ReciprocalShift;
-        UInt128 twoScaled = UInt128.One << (ReciprocalShift + 1);
-        if (scaled >= twoScaled)
-        {
-            return ComputeReciprocalEstimate(divisor);
-        }
-
-        UInt128 correction = twoScaled - scaled;
-        UInt128 refined = (reciprocal128 * correction) >> ReciprocalShift;
-        if (refined == UInt128.Zero)
-        {
-            return 1UL;
-        }
-
-        return refined > ulong.MaxValue ? ulong.MaxValue : (ulong)refined;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Reduce(UInt128 numerator, ulong divisor, ulong reciprocal)
-    {
-        if (divisor <= 1UL)
-        {
-            return 0UL;
-        }
-
-        if (numerator == UInt128.Zero)
-        {
-            return 0UL;
-        }
-
-        if (reciprocal == 0UL || reciprocal == ulong.MaxValue)
-        {
-            return (ulong)(numerator % divisor);
-        }
-
-        ulong high = (ulong)(numerator >> ReciprocalShift);
-        ulong low = (ulong)numerator;
-
-        if (high >= divisor)
-        {
-            return (ulong)(numerator % divisor);
-        }
-
-        UInt128 highProduct = (UInt128)high * reciprocal;
-        UInt128 lowProduct = (UInt128)low * reciprocal;
-        UInt128 estimate = highProduct + (lowProduct >> ReciprocalShift);
-
-        if (estimate == UInt128.Zero)
-        {
-            return (ulong)(numerator % divisor);
-        }
-
-        UInt128 product = estimate * divisor;
-        while (product > numerator)
-        {
-            if (estimate == UInt128.Zero)
-            {
-                return (ulong)(numerator % divisor);
-            }
-
-            estimate--;
-            product -= divisor;
-        }
-
-        UInt128 remainder = numerator - product;
-        while (remainder >= divisor)
-        {
-            remainder -= divisor;
-        }
-
-        return (ulong)remainder;
-    }
-}

PerfectNumbers.Core/UInt128Extensions.cs

@@ -218,7 +218,7 @@ public static class UInt128Extensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static bool IsBitSet(GpuUInt128 value, int bitIndex)
+    private static bool IsBitSet(in GpuUInt128 value, int bitIndex)
     {
         if (bitIndex >= 64)
         {
@@ -229,16 +229,18 @@ public static class UInt128Extensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong ExtractWindowValue(GpuUInt128 exponent, int windowStart, int windowBitCount)
+    private static ulong ExtractWindowValue(in GpuUInt128 exponent, int windowStart, int windowBitCount)
     {
-		ulong mask = (1UL << windowBitCount) - 1UL;
         if (windowStart != 0)
         {
-            exponent.ShiftRight(windowStart);
+            GpuUInt128 shifted = exponent;
-            return exponent.Low & mask;
+            shifted.ShiftRight(windowStart);
+            ulong mask = (1UL << windowBitCount) - 1UL;
+            return shifted.Low & mask;
         }
 
-        return exponent.Low & mask;
+        ulong directMask = (1UL << windowBitCount) - 1UL;
+        return exponent.Low & directMask;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]

PerfectNumbers.Core/ULongExtensions.cs

@@ -1,5 +1,4 @@
 using System;
-using System.Buffers;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using ILGPU.Algorithms;
@@ -327,7 +326,7 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+    public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
     {
         ulong modulus = divisor.Modulus;
         if (exponent == 0UL)
@@ -346,78 +345,7 @@ public static class ULongExtensions
 
         ulong result = divisor.MontgomeryOne;
         ulong nPrime = divisor.NPrime;
-        ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount);
+        ulong[] oddPowers = InitializeMontgomeryOddPowers(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)
@@ -533,30 +461,7 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
+    private static ulong[] InitializeMontgomeryOddPowers(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;
@@ -576,31 +481,17 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
+    public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
     {
         ulong rotationCount = exponent % cycleLength;
-        return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false);
+        return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+    public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
     {
-		return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false);
+        return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false);
     }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
-    {
-        ulong rotationCount = exponent % cycleLength;
-        return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false);
-    }
-
-	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
-	{
-		return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false);
-	}
-	
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
     {

docs/bydivisor_cpu_path.md

@@ -1,196 +0,0 @@
-# CPU `--mersenne=bydivisor` pipeline notes
-
-This document summarizes how the `EvenPerfectBitScanner` executable executes the
-`--mersenne=bydivisor` mode on the CPU, and how the analytic observations we
-derived for the quotients
-
-\[
-Q_k = \frac{N}{d_k}, \qquad d_k = 2pk + 1, \qquad N = 2^p - 1,
-\]
-
-can help us identify concrete optimization opportunities.
-
-## High-level control flow
-
-* `Program.Main` wires `_byDivisorTester` to a CPU implementation whenever the
-  CLI switches resolve to `--mersenne=bydivisor` and
-  `--mersenne-device=cpu`. The CLI also forces a block size of 1 so that the
-  pipeline pushes exponents to the tester one at a time.
-  (`Program.cs`, lines 459-520.)
-* When the mode starts, `RunByDivisorMode` delegates to
-  `MersenneNumberDivisorByDivisorTester.Run`, passing the candidate list from
-  `--filter-p`, the CPU tester instance and the optional replay of previous
-  outcomes. (`Program.cs`, lines 809-826.)
-* `MersenneNumberDivisorByDivisorTester.Run` sorts and deduplicates the
-  candidate exponents, drops the ones already confirmed in prior runs, validates
-  primality with the `Open.Numeric.Primes` enumerator and configures the tester
-  with the largest remaining exponent. (`MersenneNumberDivisorByDivisorTester.cs`,
-  lines 14-162.)
-* For each confirmed prime exponent, the CPU tester computes the largest divisor
-  it needs to probe through `ComputeAllowedMaxDivisor`. If the range is non-empty
-  it enters the `CheckDivisors` loop. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
-  lines 83-110 and 161-266.)
-
-## Inner loop behaviour (`CheckDivisors`)
-
-The hot path is the `while (divisor <= limit)` loop in
-`MersenneNumberDivisorByDivisorCpuTester.CheckDivisors`.
-
-1. The loop enumerates divisor candidates of the form `d_k = 2pk + 1` by
-   incrementing `divisor` by `step = 2p` every iteration. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
-   lines 177-256.)
-2. Small modulus filters track `d_k (mod 3·5·7·8·10·11)` to discard impossible
-   candidates without touching the divisor-cycle cache. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
-   lines 195-263.)
-3. For every admissible divisor the tester fetches its Montgomery data and the
-   cached cycle length. When the cache does not contain an entry the
-   implementation recomputes the cycle by factoring the divisor. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
-   lines 223-254.)
-4. A cycle equal to the tested exponent proves divisibility, so the loop aborts
-   as soon as one hit is found. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
-   lines 241-247.)
-
-`DivisorScanSession.CheckDivisor` follows the same `d_k` progression when a
-batch API is used. It keeps a running Montgomery residue via
-`ExponentRemainderStepper`, but the stepper currently recomputes the delta by
-calling a full `Pow2MontgomeryModWindowed` whenever the exponent jumps.
-(`MersenneNumberDivisorByDivisorCpuTester.cs`, lines 358-405 and
-`ExponentRemainderStepper.cs`, lines 26-95.)
-
-## Where the analytic model fits
-
-The formulas we established for the quotient sequence
-
-\[
-\frac{Q_{k+1}}{Q_k} = \frac{2pk + 1}{2p(k+1) + 1}, \qquad
-Q_k - Q_{k+1} = \frac{N·2p}{(2pk + 1)(2p(k+1) + 1)}
-\]
-
-tell us that the successive quotients and remainders shrink smoothly and that the
-per-step update depends only on the previous `(Q_k, r_k)` state and the next
-`d_{k+1}`. This lines up with several TODO comments that highlight missing
-incremental updates:
-
-* `ExponentRemainderStepper` already caches the previous exponent and a running
-  Montgomery residue, but it falls back to a full modular exponentiation to
-  compute the delta between consecutive divisors. (`ExponentRemainderStepper.cs`,
-  lines 49-72.)
-* `CheckDivisors` rebuilds the Montgomery data for every hit candidate even
-  though the sequence of `d_k` values is strictly monotone.
-
-The quotient identity gives us a way to compute `Δ_k = Q_k - Q_{k+1}` by reading
-`Q_k` and the inexpensive division `(2pQ_k - r_k) / d_{k+1}`. Once we know `Δ_k`
-we can update
-
-\[
-Q_{k+1} = Q_k - Δ_k, \qquad r_{k+1} = r_k - 2p·Q_k + Δ_k · d_{k+1},
-\]
-
-without ever evaluating `2^p mod d_{k+1}` from scratch. Incorporating this into
-`ExponentRemainderStepper` would let the CPU path reuse a single Montgomery state
-across the entire scan and swap the current delta exponentiation (which costs a
-full modular ladder) for one multiply-high with a precomputed reciprocal of
-`d_{k+1}`. The loop already maintains `d_k` and `step`, so the extra bookkeeping
-fits naturally.
-
-## Optimization ideas grounded in the code
-
-1. **Delta-driven Montgomery updates.** Extend `DivisorScanSession` so it tracks
-   `(Q_k, r_k)` alongside the Montgomery residue and uses the analytic delta
-   formula to derive the next residue. This replaces the per-step powmod in
-   `ExponentRemainderStepper.ComputeNext` with a handful of 64-bit operations.
-   (`MersenneNumberDivisorByDivisorCpuTester.cs`, lines 358-405;
-   `ExponentRemainderStepper.cs`, lines 49-72.)
-2. **Bit-range skipping.** Because `Q_{k+1} / Q_k` is exactly the rational
-   `(2pk + 1)/(2p(k+1) + 1)`, we can predict when `⌊log₂ Q_k⌋` loses a bit and
-   jump the loop counter to the first `k` that crosses a target quotient bound
-   instead of checking each divisor. This is compatible with the existing
-   `allowedMax` guard because it only modifies how we arrive at the stopping
-   divisor. (`MersenneNumberDivisorByDivisorCpuTester.cs`, lines 161-266.)
-3. **Reuse Montgomery metadata.** The loop already calls `MontgomeryDivisorDataCache.Get`
-   for each divisor, but the cache lookup is repeated even though the next
-   divisor is a small increment away. Persisting the previous `MontgomeryDivisorData`
-   together with the `(Q_k, r_k)` tuple lets us skip redundant cache reads when
-   the divisor stays within the hot range. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
-   lines 223-256.)
-4. **Reciprocal-based delta evaluation.** `DivisorScanSession.TryComputeAnalyticHit` already keeps the `(Q_k, r_k)` tuple
-   per exponent, but the current update loop advances one step at a time and divides by `nextDivisor` through `Int128`
-   arithmetic. (`MersenneNumberDivisorByDivisorCpuTester.cs`, lines 531-594.) We can remove the per-iteration division by
-   storing a 64-bit reciprocal of `d_k` alongside the state and updating it with a single Newton step when the divisor grows
-   by `2p`. Once the reciprocal is known, `Δ_k` follows from a multiply-high of `2p·Q_k - r_k`. This turns the loop into a
-   straight-line sequence of multiplies and subtractions and makes the state update viable even when we unroll multiple
-   increments at once.
-5. **Carry the Montgomery residue forward.** `PrimeDivisorState` already tracks `MontgomeryResidue`, yet the analytic path
-   currently recomputes it after the quotient update by calling `ToMontgomeryResidue`. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
-   lines 602-618.) We can instead reuse the residue by multiplying it with the precomputed `R^2 mod d_k` factor that the
-   cache returns for every divisor. Because the divisor grows deterministically, we only need to fetch the updated Montgomery
-   constants when `d_k` crosses a cache boundary, letting the fast path avoid modular reductions entirely.
-
-These adjustments directly leverage the quotient progression we analysed and
-align with the existing TODO list, giving us a clear path to reduce the number of
-modular exponentiations and cache lookups in the CPU by-divisor path.
-
-## Implementation status
-
-* **Delta-driven Montgomery updates:** Implemented. The CPU session now keeps per-prime quotients, remainders, and last divisors so it can advance the analytic state across consecutive candidates, the reciprocal-corrected loop avoids per-step `UInt128` division, and zero cycle remainders short-circuit to direct Montgomery hits without invoking powmods.
-* **Bit-range skipping:** Implemented. The CPU divisor sweep now predicts when the quotient loses a bit across both 64-bit and wider strides by deriving the next jump with `UInt128` arithmetic,
-  and it only falls back to single-step increments when a projected jump would overshoot the configured limit or exceed `ulong.MaxValue`.
-* **Reuse Montgomery metadata:** Implemented. Each prime state now persists the divisor's `n'`, `R`, and `R^2` factors so resynchronization seeds the exponent stepper directly, and divisor sessions rebuild the Montgomery tuple from the cached state before touching the shared cache, meaning cache lookups only occur when no stored snapshot matches the new divisor.
-* **Reciprocal-based delta evaluation:** Implemented. The analytic updater now retries with a direct 128-bit division when the reciprocal corrections exhaust their tolerance, refreshes the cached reciprocal from that exact quotient, and keeps the progression on the analytic path without rebuilding the state.
-* **Carry the Montgomery residue forward:** Implemented. Each analytic update now normalizes the refreshed remainder and multiplies it by the cached `R^2 mod d` factor immediately, so the exponent stepper always has a Montgomery residue ready without recomputing modular reductions.
-
-
-## Prime-order factoring bottleneck and planned fixes
-
-The Pollard rho fallback in `PrimeOrderCalculator` still performs every quadratic
-polynomial advance by issuing a full-width division. The tight loop in
-`AdvancePolynomial` squares the iterate and takes a `% modulus`, which triggers
-hardware division for each step and dominates the runtime once we test large
-Mersenne divisors.【F:PerfectNumbers.Core/PrimeOrderCalculator.cs†L1180-L1211】
-The by-divisor CPU path consequently spends a measurable portion of its wall
-clock time inside the factorization helper whenever it needs to refresh a cycle
-or validate that a divisor candidate is prime.
-
-To bring this code in line with the Montgomery-first approach we use elsewhere
-we should:
-
-* **Introduce a Montgomery-backed polynomial step.** `PrimeDivisorState` already
-  carries the full set of Montgomery constants (`n'`, `R`, `R^2`) for each
-  divisor we touch in the by-divisor scan, and `DivisorScanSession` rebuilds a
-  `MontgomeryDivisorData` snapshot whenever it seeds the analytic stepper or the
-  exponent synchronizer.【F:PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs†L645-L779】【F:PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs†L1246-L1398】
-  We can surface the same data to the Pollard rho helper and rewrite
-  `AdvancePolynomial` so it accepts a precomputed reducer `(modulus, n', R^2)`.
-  The update then becomes `x = MontgomerySquareAdd(x, cMontgomery)`—a multiply
-  followed by one Montgomery reduction instead of a general-purpose remainder.
-  This eliminates the `% modulus` call entirely and lets the CPU factorizer reuse
-  the same fast path as the analytic delta updates.
-* **Reuse reciprocal estimates between neighbouring divisors.** The analytic path
-  already stores a reciprocal for each divisor and refines it with Newton steps
-  to divide by `d_k` without hitting the 128-bit fallback.【F:PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs†L876-L939】
-  Pollard rho can share the very same multiplier. We can stash the `FastDiv64`
-  multiplier inside the Pollard rho state and update it when `q` grows by `2p`
-  using a single Newton refinement, mirroring `DivideRoundUpUsingReciprocal`.
-  Once the reciprocal is in place, the `(x * x + c) mod q` step reduces to a
-  pair of `UInt128` products and one multiply-high to estimate the quotient—a
-  drop-in use of `ULongExtensions.FastDiv64`.【F:PerfectNumbers.Core/ULongExtensions.cs†L69-L101】
-* **Warm-start the Montgomery context from the monotone `q` stream.** Each scan
-  step increments the divisor by the fixed stride `2p`, so the cached Montgomery
-  tuple only changes slowly. We can therefore hand the next divisor’s
-  `MontgomeryDivisorData` directly to the Pollard rho routine instead of forcing
-  it to recompute the reducer. This is already how `DivisorScanSession` avoids
-  redundant cache lookups; exporting the cached tuple keeps the factorizer in
-  lock-step with the rest of the pipeline.【F:PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs†L645-L779】
-
-With these changes the divisor-refresh path would no longer spend most of its
-cycles on integer divisions, and the monotone growth of `q = 2kp + 1` for
-`p ≥ 138,000,000` would guarantee that the reciprocal refinements stay stable
-across thousands of consecutive candidates.
-
-### Pollard rho factoring implementation status
-
-* **Montgomery-backed polynomial step:** Completed. `PollardRhoMontgomeryReducer` now squares and adds in Montgomery space,
-  and the Pollard rho loops in `PrimeOrderCalculator` and `MersenneDivisorCycles` reuse the shared reducer instead of issuing `%` operations.
-* **Reciprocal reuse:** Completed. The reducer tracks a cached FastDiv64 reciprocal, refines it for each monotone modulus, and exposes `AdvanceStandard` so the standard-domain iterates reuse the cached multiplier without `%` or `UInt128` division.
-* **Warm-started Montgomery context:** Completed. The reducer accepts precomputed `MontgomeryDivisorData` snapshots, letting callers seed Pollard rho directly from the cached divisor state when the by-divisor pipeline advances to the next `q`.
-

todo_list.txt

@@ -123,6 +123,7 @@ PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207:
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214:                // TODO: Remove this no-cycle branch and require on-demand cycle computation (without
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240:            // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262:                : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the
+PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283:        // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356:                // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396:            // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here
 PerfectNumbers.Core\IMersenneNumberDivisorByDivisorTester.cs:7:    bool UseDivisorCycles { get; set; } // TODO: Remove the setter once divisor cycle usage becomes mandatory so all implementations always leverage the cached cycles.
@@ -266,6 +267,7 @@ todo_list.txt:122:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTeste
 todo_list.txt:123:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214:                // TODO: Remove this no-cycle branch and require on-demand cycle computation (without
 todo_list.txt:124:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240:            // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the
 todo_list.txt:125:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262:                : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the
+todo_list.txt:126:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283:        // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
 todo_list.txt:127:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356:                // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips
 todo_list.txt:128:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396:            // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here
 todo_list.txt:129:PerfectNumbers.Core\IMersenneNumberDivisorByDivisorTester.cs:7:    bool UseDivisorCycles { get; set; } // TODO: Remove the setter once divisor cycle usage becomes mandatory so all implementations always leverage the cached cycles.