Merge pull request #144 from mlemanczyk/codex/optimize-division-in-primeordercalculator

Cache Pollard rho reciprocals for faster Montgomery steps

EvenPerfectBitScanner.Benchmarks/EvenPerfectBitScanner.Benchmarks.csproj

@@ -15,7 +15,4 @@
     <PackageReference Include="BenchmarkDotNet" Version="0.13.12" />
   </ItemGroup>
 
-  <ItemGroup>
-    <Compile Remove=".temp.cs" />
-  </ItemGroup>
 </Project>

EvenPerfectBitScanner.Benchmarks/MutableUInt128AddBenchmarks.cs

@@ -1,75 +0,0 @@
-using System.Runtime.CompilerServices;
-using BenchmarkDotNet.Attributes;
-using BenchmarkDotNet.Jobs;
-using PerfectNumbers.Core;
-
-namespace EvenPerfectBitScanner.Benchmarks;
-
-[ShortRunJob(RuntimeMoniker.Net80)]
-[MemoryDiagnoser]
-public class MutableUInt128AddBenchmarks
-{
-	private static readonly BenchmarkInput[] Inputs =
-	[
-		new(Create(0UL, 0x0000_0000_0000_0007UL), Create(0, 0x3UL), "TinyOperands"),
-		new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), Create(0xFFFF_FFFF, 0x0000_0001UL), "MixedOperands"),
-		new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), Create(0xFFFF_FFFF,  0xFFFF_FFFBUL), "NearMaxOperands"),
-    ];
-
-    [ParamsSource(nameof(GetInputs))]
-    public BenchmarkInput Input { get; set; }
-
-    public static IEnumerable<BenchmarkInput> GetInputs() => Inputs;
-
-	private static MutableUInt128 Create(ulong high, ulong low) => new(high, low);
-
-	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	private static void StandardInternal(ref MutableUInt128 source, MutableUInt128 value)
-	{
-        ulong temp = source.Low + value.Low;
-        source.Low = temp;
-        temp = temp < source.Low ? 1UL : 0UL;
-		temp += source.High;
-		temp += value.High;
-        source.High = temp;
-	}
-
-	/// <summary>
-	/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.370 ns with mixed operands,
-	/// 3.361 ns when both limbs and inputs are near their maxima, and 3.368 ns with tiny operands.
-	/// </summary>
-	/// <remarks>
-	/// Observed means: MixedOperands 3.370 ns, NearMaxOperands 3.361 ns, TinyOperands 3.368 ns.
-	/// </remarks>
-	[Benchmark(Baseline = true)]
-	public MutableUInt128 Standard()
-	{
-		MutableUInt128 value = Input.InitialValue;
-		StandardInternal(ref value, Input.Addend);
-		return value;
-	}
-
-    /// <summary>
-    /// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 3.175 ns on mixed operands,
-    /// 3.177 ns on near-max inputs, and 3.178 ns on tiny operands.
-    /// </summary>
-    /// <remarks>
-    /// Observed means: MixedOperands 3.175 ns, NearMaxOperands 3.177 ns, TinyOperands 3.178 ns.
-	/// Status: Implemented
-    /// </remarks>
-    [Benchmark]
-    public MutableUInt128 Optimized()
-    {
-        MutableUInt128 value = Input.InitialValue;
-        value.Add(Input.Addend);
-        return value;
-    }
-
-    public readonly record struct BenchmarkInput(MutableUInt128 InitialValue, MutableUInt128 Addend, string Name)
-    {
-        public override string ToString()
-        {
-            return Name;
-        }
-    }
-}

EvenPerfectBitScanner.Benchmarks/MutableUInt128Mul64Benchmarks.cs

@@ -1,100 +0,0 @@
-using System.Runtime.CompilerServices;
-using BenchmarkDotNet.Attributes;
-using BenchmarkDotNet.Jobs;
-using PerfectNumbers.Core;
-
-namespace EvenPerfectBitScanner.Benchmarks;
-
-[ShortRunJob(RuntimeMoniker.Net80)]
-[MemoryDiagnoser]
-public class MutableUInt128Mul64Benchmarks
-{
-	private static readonly BenchmarkInput[] Inputs =
-	[
-		new(Create(0UL, 0x0000_0000_0000_0007UL), 0x3UL, "TinyOperands"),
-		new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), 0xFFFF_FFFF_0000_0001UL, "MixedOperands"),
-		new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), 0xFFFF_FFFF_FFFF_FFFBUL, "NearMaxOperands"),
-    ];
-
-    [ParamsSource(nameof(GetInputs))]
-    public BenchmarkInput Input { get; set; }
-
-    public static IEnumerable<BenchmarkInput> GetInputs() => Inputs;
-
-	private static MutableUInt128 Create(ulong high, ulong low) => new(high, low);
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong Mul64(ulong left, ulong right, out ulong low)
-    {
-        const ulong Mask32 = 0xFFFF_FFFFUL;
-
-        ulong leftLow = left & Mask32;
-        ulong leftHigh = left >> 32;
-        ulong rightLow = right & Mask32;
-        ulong rightHigh = right >> 32;
-
-        ulong lowLow = leftLow * rightLow;
-        ulong lowHigh = leftLow * rightHigh;
-        ulong highLow = leftHigh * rightLow;
-        ulong highHigh = leftHigh * rightHigh;
-
-        ulong carry = (lowLow >> 32) + (lowHigh & Mask32) + (highLow & Mask32);
-        ulong high = highHigh + (lowHigh >> 32) + (highLow >> 32) + (carry >> 32);
-
-        low = ((carry & Mask32) << 32) | (lowLow & Mask32);
-        return high;
-    }
-
-	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	private static void StandardInternal(ref MutableUInt128 source, ulong value)
-	{
-        ulong currentHigh = source.High;
-        ulong currentLow = source.Low;
-
-        ulong carry = Mul64(currentLow, value, out ulong newLow);
-        Mul64(currentHigh, value, out ulong shiftedHigh);
-        ulong newHigh = shiftedHigh + carry;
-
-        source.Low = newLow;
-        source.High = newHigh;
-	}
-
-	/// <summary>
-	/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.113 ns with mixed operands,
-	/// 3.128 ns when both limbs and inputs are near their maxima, and 3.131 ns with tiny operands.
-	/// </summary>
-	/// <remarks>
-	/// Observed means: MixedOperands 3.113 ns, NearMaxOperands 3.128 ns, TinyOperands 3.131 ns.
-	/// </remarks>
-	[Benchmark(Baseline = true)]
-	public MutableUInt128 Standard()
-	{
-		MutableUInt128 value = Input.InitialValue;
-		StandardInternal(ref value, Input.Addend);
-		return value;
-	}
-
-    /// <summary>
-    /// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 2.969 ns on mixed operands,
-    /// 2.912 ns on near-max inputs, and 2.935 ns on tiny operands.
-    /// </summary>
-    /// <remarks>
-    /// Observed means: MixedOperands 2.969 ns, NearMaxOperands 2.912 ns, TinyOperands 2.935 ns.
-	/// Status: Implemented
-    /// </remarks>
-    [Benchmark]
-    public MutableUInt128 Optimized()
-    {
-        MutableUInt128 value = Input.InitialValue;
-        value.Mul(Input.Addend);
-        return value;
-    }
-
-    public readonly record struct BenchmarkInput(MutableUInt128 InitialValue, ulong Addend, string Name)
-    {
-        public override string ToString()
-        {
-            return Name;
-        }
-    }
-}

EvenPerfectBitScanner.Benchmarks/MutableUInt128MultiplyAddBenchmarks.cs

@@ -1,76 +0,0 @@
-using System.Collections.Generic;
-using System.Numerics;
-using System.Runtime.CompilerServices;
-using BenchmarkDotNet.Attributes;
-using BenchmarkDotNet.Jobs;
-using PerfectNumbers.Core;
-
-namespace EvenPerfectBitScanner.Benchmarks;
-
-[ShortRunJob(RuntimeMoniker.Net80)]
-[MemoryDiagnoser]
-public class MutableUInt128MultiplyAddBenchmarks
-{
-    private static readonly MultiplyAddInput[] Inputs =
-    [
-        new(Create(0UL, 0x0000_0000_0000_0007UL), 0x3UL, 0x5UL, "TinyOperands"),
-        new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), 0xFFFF_FFFF_0000_0001UL, 0x1357_9BDF_2468_ACEEUL, "MixedOperands"),
-        new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), 0xFFFF_FFFF_FFFF_FFFBUL, 0xFFFF_FFFF_FFFF_FFF7UL, "NearMaxOperands"),
-    ];
-
-    [ParamsSource(nameof(GetInputs))]
-    public MultiplyAddInput Input { get; set; }
-
-    public static IEnumerable<MultiplyAddInput> GetInputs() => Inputs;
-
-    private static UInt128 Create(ulong high, ulong low)
-    {
-        return ((UInt128)high << 64) | low;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-	private static void MultiplyAdd(ref MutableUInt128 source, ulong multiplier, ulong addend)
-	{
-		source.Mul(multiplier);
-		source.Add(addend);
-	}
-
-	/// <summary>
-	/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.370 ns with mixed operands,
-	/// 3.361 ns when both limbs and inputs are near their maxima, and 3.368 ns with tiny operands.
-	/// </summary>
-	/// <remarks>
-	/// Observed means: MixedOperands 3.370 ns, NearMaxOperands 3.361 ns, TinyOperands 3.368 ns.
-	/// </remarks>
-	[Benchmark(Baseline = true)]
-	public MutableUInt128 MultiplyAddChained()
-	{
-		MutableUInt128 value = new(Input.InitialValue);
-		MultiplyAdd(ref value, Input.Multiplier, Input.Addend);
-		return value;
-	}
-
-    /// <summary>
-    /// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 3.175 ns on mixed operands,
-    /// 3.177 ns on near-max inputs, and 3.178 ns on tiny operands.
-    /// </summary>
-    /// <remarks>
-    /// Observed means: MixedOperands 3.175 ns, NearMaxOperands 3.177 ns, TinyOperands 3.178 ns.
-	/// Status: Implemented
-    /// </remarks>
-    [Benchmark]
-    public MutableUInt128 MultiplyAddFullyInlined()
-    {
-        MutableUInt128 value = new(Input.InitialValue);
-        value.MulAdd(Input.Multiplier, Input.Addend);
-        return value;
-    }
-
-    public readonly record struct MultiplyAddInput(UInt128 InitialValue, ulong Multiplier, ulong Addend, string Name)
-    {
-        public override string ToString()
-        {
-            return Name;
-        }
-    }
-}

EvenPerfectBitScanner.Benchmarks/MutableUInt128SubMutableBenchmarks.cs

@@ -1,75 +0,0 @@
-using System.Runtime.CompilerServices;
-using BenchmarkDotNet.Attributes;
-using BenchmarkDotNet.Jobs;
-using PerfectNumbers.Core;
-
-namespace EvenPerfectBitScanner.Benchmarks;
-
-[ShortRunJob(RuntimeMoniker.Net80)]
-[MemoryDiagnoser]
-public class MutableUInt128SubMutableBenchmarks
-{
-	private static readonly BenchmarkInput[] Inputs =
-	[
-		new(Create(0UL, 0x0000_0000_0000_0007UL), Create(0, 0x3UL), "TinyOperands"),
-		new(Create(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL), Create(0xFFFF_FFFF, 0x0000_0001UL), "MixedOperands"),
-		new(Create(0xFFFF_FFFF_FFFF_FFFFUL, 0xFFFF_FFFF_FFFF_FFFFUL), Create(0xFFFF_FFFF,  0xFFFF_FFFBUL), "NearMaxOperands"),
-    ];
-
-    [ParamsSource(nameof(GetInputs))]
-    public BenchmarkInput Input { get; set; }
-
-    public static IEnumerable<BenchmarkInput> GetInputs() => Inputs;
-
-	private static MutableUInt128 Create(ulong high, ulong low) => new(high, low);
-
-	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	private static void StandardInternal(ref MutableUInt128 source, MutableUInt128 value)
-	{
-        ulong low = source.Low;
-        ulong result = low - value.Low;
-        ulong borrow = low < value.Low ? 1UL : 0UL;
-
-        source.Low = result;
-        source.High = source.High - value.High - borrow;
-	}
-
-	/// <summary>
-	/// Baseline MultiplyAdd implementation that chains Multiply and Add; measured 3.370 ns with mixed operands,
-	/// 3.361 ns when both limbs and inputs are near their maxima, and 3.368 ns with tiny operands.
-	/// </summary>
-	/// <remarks>
-	/// Observed means: MixedOperands 3.370 ns, NearMaxOperands 3.361 ns, TinyOperands 3.368 ns.
-	/// </remarks>
-	[Benchmark(Baseline = true)]
-	public MutableUInt128 Standard()
-	{
-		MutableUInt128 value = Input.InitialValue;
-		StandardInternal(ref value, Input.Addend);
-		return value;
-	}
-
-    /// <summary>
-    /// Fully inlined MultiplyAdd implementation that kept all arithmetic inside the method; landed at 3.175 ns on mixed operands,
-    /// 3.177 ns on near-max inputs, and 3.178 ns on tiny operands.
-    /// </summary>
-    /// <remarks>
-    /// Observed means: MixedOperands 3.175 ns, NearMaxOperands 3.177 ns, TinyOperands 3.178 ns.
-	/// Status: Implemented
-    /// </remarks>
-    [Benchmark]
-    public MutableUInt128 Optimized()
-    {
-        MutableUInt128 value = Input.InitialValue;
-        value.Sub(Input.Addend);
-        return value;
-    }
-
-    public readonly record struct BenchmarkInput(MutableUInt128 InitialValue, MutableUInt128 Addend, string Name)
-    {
-        public override string ToString()
-        {
-            return Name;
-        }
-    }
-}

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs

@@ -377,11 +377,10 @@ public class Pow2MontgomeryModCycleComputationBenchmarks
 #if DEBUG
 		Console.WriteLine("Trying heuristic. Wide prime order calculation");
 #endif
-        UInt128 order = PrimeOrderCalculator.Calculate(
-                modulus,
-                _previousWidePrimeOrder,
-                PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
-                PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
+		UInt128 order = PrimeOrderCalculator.Calculate(
+			modulus,
+			_previousWidePrimeOrder,
+			PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
 
 		if (order != UInt128.Zero)
 		{

EvenPerfectBitScanner.Benchmarks/Program.cs

@@ -8,22 +8,16 @@ public static class Program
 {
     public static void Main(string[] args)
     {
-		if (args.Length > 0)
-		{
-			IConfig config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary | ConfigOptions.DisableOptimizationsValidator);
-			BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config);
-			return;
-		}
+        // var config = DefaultConfig.Instance.WithOptions(ConfigOptions.JoinSummary | ConfigOptions.DisableOptimizationsValidator);
+        // BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run(args, config);
 
-		BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly).Run();
+        Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new()
+        {
+            Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large
+        };
 
-        // Pow2MontgomeryModCycleComputationBenchmarks benchmarks = new()
-        // {
-        //     Scale = Pow2MontgomeryModCycleComputationBenchmarks.InputScale.Large
-        // };
-
-        // GpuPrimeWorkLimiter.SetLimit(1);
-        // benchmarks.Setup();
+        GpuPrimeWorkLimiter.SetLimit(1);
+        benchmarks.Setup();
 
         // Console.WriteLine("Press ENTER to run");
         // Console.ReadLine();

EvenPerfectBitScanner.ResultsParser/EvenPerfectBitScanner.ResultsParser.csproj

@@ -11,7 +11,4 @@
     <PackageReference Include="Open.Numeric.Primes" Version="4.0.4" />
   </ItemGroup>
 
-  <ItemGroup>
-    <Compile Remove=".temp.cs" />
-  </ItemGroup>
 </Project>

EvenPerfectBitScanner.Tests/EvenPerfectBitScanner.Tests.csproj

@@ -22,7 +22,4 @@
   <ItemGroup>
     <ProjectReference Include="../EvenPerfectBitScanner/EvenPerfectBitScanner.csproj" />
   </ItemGroup>
-  <ItemGroup>
-    <Compile Remove=".temp.cs" />
-  </ItemGroup>
 </Project>

EvenPerfectBitScanner/EvenPerfectBitScanner.csproj

@@ -16,7 +16,4 @@
     <PackageReference Include="Open.Numeric" Version="3.0.1" />
     <PackageReference Include="Open.Numeric.Primes" Version="4.0.4" />
   </ItemGroup>
-  <ItemGroup>
-    <Compile Remove=".temp.cs" />
-  </ItemGroup>
 </Project>

EvenPerfectBitScanner/Program.cs

@@ -44,7 +44,6 @@ internal static class Program
     private static string? _resultsPrefix;
     private static readonly Optimized PrimeIterator = new();
     private static ulong _byDivisorStartPrime;
-    private static TimeSpan? _primeTestLimit;
 
     [ThreadStatic]
     private static bool _lastCompositeP;
@@ -75,8 +74,6 @@ internal static class Program
         bool useResidue = false;    // M_p test via residue divisors (replaces LL/incremental)
         bool useDivisor = false;     // M_p divisibility by specific divisor
         bool useByDivisor = false;   // Iterative divisor scan across primes		UInt128 divisor = UInt128.Zero;
-        ByDivisorDeltasDevice byDivisorDeltasDevice = ByDivisorDeltasDevice.Cpu;
-        ByDivisorMontgomeryDevice byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Cpu;
                                      // Device routing
         bool useGpuCycles = true;
         bool mersenneOnGpu = true;   // controls Lucas/incremental/pow2mod device
@@ -184,51 +181,6 @@ internal static class Program
             {
                 mersenneOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase);
             }
-            else if (arg.StartsWith("--bydivisor-deltas-device=", StringComparison.OrdinalIgnoreCase))
-            {
-                ReadOnlySpan<char> deviceValue = arg.AsSpan(arg.IndexOf('=') + 1);
-                if (deviceValue.Equals("gpu", StringComparison.OrdinalIgnoreCase))
-                {
-                    byDivisorDeltasDevice = ByDivisorDeltasDevice.Gpu;
-                }
-                else if (deviceValue.Equals("cpu", StringComparison.OrdinalIgnoreCase))
-                {
-                    byDivisorDeltasDevice = ByDivisorDeltasDevice.Cpu;
-                }
-                else
-                {
-                    Console.WriteLine("Invalid --bydivisor-deltas-device value. Expected cpu or gpu.");
-                    return;
-                }
-            }
-            else if (arg.StartsWith("--bydivisor-montgomery-device=", StringComparison.OrdinalIgnoreCase))
-            {
-                ReadOnlySpan<char> deviceValue = arg.AsSpan(arg.IndexOf('=') + 1);
-                if (deviceValue.Equals("gpu", StringComparison.OrdinalIgnoreCase))
-                {
-                    byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Gpu;
-                }
-                else if (deviceValue.Equals("cpu", StringComparison.OrdinalIgnoreCase))
-                {
-                    byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Cpu;
-                }
-                else
-                {
-                    Console.WriteLine("Invalid --bydivisor-montgomery-device value. Expected cpu or gpu.");
-                    return;
-                }
-            }
-            else if (arg.StartsWith("--prime-test-limit=", StringComparison.OrdinalIgnoreCase))
-            {
-                ReadOnlySpan<char> value = arg.AsSpan(arg.IndexOf('=') + 1);
-                if (!TryParsePrimeTestLimit(value, out TimeSpan parsedLimitValue))
-                {
-                    Console.WriteLine("Failed to parse --prime-test-limit value.");
-                    return;
-                }
-
-                _primeTestLimit = parsedLimitValue;
-            }
             else if (arg.Equals("--workaround-mod", StringComparison.OrdinalIgnoreCase))
             {
                 useModuloWorkaround = true;
@@ -425,8 +377,6 @@ internal static class Program
             }
         }
 
-        GpuKernelPool.CpuCycleUsesGpu = orderOnGpu;
-
         if (useByDivisor)
         {
             _byDivisorStartPrime = startPrimeProvided ? currentP : 0UL;
@@ -567,11 +517,6 @@ internal static class Program
                     ? new MersenneNumberDivisorByDivisorGpuTester()
                     : new MersenneNumberDivisorByDivisorCpuTester();
             _byDivisorTester.BatchSize = scanBatchSize;
-            if (_byDivisorTester is MersenneNumberDivisorByDivisorCpuTester cpuByDivisorTester)
-            {
-                cpuByDivisorTester.DeltasDevice = byDivisorDeltasDevice;
-                cpuByDivisorTester.MontgomeryDevice = byDivisorMontgomeryDevice;
-            }
         }
 
         // Load RLE blacklist (optional)
@@ -607,18 +552,18 @@ internal static class Program
                         useDivisor,
                         useByDivisor,
                         mersenneOnGpu,
-                        useOrder,
-                        useModuloWorkaround,
-                        _useGcdFilter,
-                        NttGpuMath.GpuTransformBackend,
-                        gpuPrimeThreads,
-                        sliceSize,
-                        scanBatchSize,
-                        _orderWarmupLimitOverride ?? 5_000_000UL,
-                        NttGpuMath.ReductionMode,
-                        mersenneOnGpu ? "gpu" : "cpu",
-                        (GpuContextPool.ForceCpu ? "cpu" : "gpu"),
-                        orderOnGpu ? "gpu" : "cpu");
+        useOrder,
+        useModuloWorkaround,
+        _useGcdFilter,
+        NttGpuMath.GpuTransformBackend,
+        gpuPrimeThreads,
+        sliceSize,
+        scanBatchSize,
+        _orderWarmupLimitOverride ?? 5_000_000UL,
+        NttGpuMath.ReductionMode,
+        mersenneOnGpu ? "gpu" : "cpu",
+        (GpuContextPool.ForceCpu ? "cpu" : "gpu"),
+        orderOnGpu ? "gpu" : "cpu");
 
         if (!string.IsNullOrEmpty(_resultsPrefix))
         {
@@ -720,7 +665,7 @@ internal static class Program
 
         if (_useByDivisorMode)
         {
-            RunByDivisorMode(byDivisorCandidates, threadCount, _primeTestLimit);
+            RunByDivisorMode(byDivisorCandidates, threadCount);
             FlushBuffer();
             StringBuilderPool.Return(_outputBuilder!);
             return;
@@ -861,7 +806,7 @@ internal static class Program
         return candidates;
     }
 
-    private static void RunByDivisorMode(List<ulong> candidates, int threadCount, TimeSpan? primeTestLimit)
+    private static void RunByDivisorMode(List<ulong> candidates, int threadCount)
     {
         if (_byDivisorTester is null)
         {
@@ -878,8 +823,7 @@ internal static class Program
                 static () => _lastCompositeP = true,
                 static () => _lastCompositeP = false,
                 PrintResult,
-                threadCount,
-                primeTestLimit);
+                threadCount);
     }
 
     private static unsafe void LoadResultsFile(string resultsFileName, Action<ulong, bool, bool> lineProcessorAction)
@@ -946,51 +890,6 @@ internal static class Program
         }
     }
 
-    private static bool TryParsePrimeTestLimit(ReadOnlySpan<char> value, out TimeSpan result)
-    {
-        result = TimeSpan.Zero;
-        if (value.IsEmpty)
-        {
-            return false;
-        }
-
-        ReadOnlySpan<char> magnitudeSpan = value;
-        double magnitude;
-
-        if (value.Length >= 2 && value.EndsWith("ms", StringComparison.OrdinalIgnoreCase))
-        {
-            magnitudeSpan = value[..^2];
-            if (!double.TryParse(magnitudeSpan, NumberStyles.Float, CultureInfo.InvariantCulture, out magnitude))
-            {
-                return false;
-            }
-
-            result = TimeSpan.FromMilliseconds(magnitude);
-            return true;
-        }
-
-        char suffix = char.ToLowerInvariant(value[^1]);
-        if (suffix == 's' || suffix == 'm')
-        {
-            magnitudeSpan = value[..^1];
-            if (!double.TryParse(magnitudeSpan, NumberStyles.Float, CultureInfo.InvariantCulture, out magnitude))
-            {
-                return false;
-            }
-
-            result = suffix == 's' ? TimeSpan.FromSeconds(magnitude) : TimeSpan.FromMinutes(magnitude);
-            return true;
-        }
-
-        if (!double.TryParse(value, NumberStyles.Float, CultureInfo.InvariantCulture, out magnitude))
-        {
-            return false;
-        }
-
-        result = TimeSpan.FromMilliseconds(magnitude);
-        return true;
-    }
-
     private static void PrintHelp()
     {
         Console.WriteLine("Usage: EvenPerfectBitScanner [options]");
@@ -1024,9 +923,6 @@ internal static class Program
         Console.WriteLine("  --divisor-cycles-batch=<value> batch size for cycles generation (default 512)");
         Console.WriteLine("  --divisor-cycles-continue  continue divisor cycles generation");
         Console.WriteLine("  --divisor-cycles-limit=<value> cycle search iterations when --mersenne=divisor");
-        Console.WriteLine("  --prime-test-limit=<value>   Time limit for --mersenne=bydivisor prime checks (e.g. 5s, 500ms)");
-        Console.WriteLine("  --bydivisor-deltas-device=cpu|gpu  Device for --mersenne=bydivisor residue batches (default: cpu)");
-        Console.WriteLine("  --bydivisor-montgomery-device=cpu|gpu  Device for Montgomery data when using --mersenne=bydivisor (default: cpu)");
         Console.WriteLine("  --use-order            test primality via q order");
         Console.WriteLine("  --workaround-mod       avoid '%' operator on the GPU");
         // mod-automaton removed

PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs

@@ -15,15 +15,9 @@ public class GpuContextPoolTests
     {
         GpuContextPool.DisposeAll();
 
-        GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
-        {
-            lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
-        }
-        finally
-        {
-            lease.Dispose();
-        }
+        using var lease = GpuContextPool.RentPreferred(preferCpu: true);
+
+        lease.Accelerator.AcceleratorType.Should().Be(AcceleratorType.CPU);
 
         GpuContextPool.DisposeAll();
     }
@@ -36,27 +30,17 @@ public class GpuContextPoolTests
 
         Context firstContext;
         Accelerator firstAccelerator;
-        GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
+        using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
         {
             firstContext = lease.Context;
             firstAccelerator = lease.Accelerator;
         }
-        finally
-        {
-            lease.Dispose();
-        }
 
-        lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
+        using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
         {
             lease.Context.Should().BeSameAs(firstContext);
             lease.Accelerator.Should().BeSameAs(firstAccelerator);
         }
-        finally
-        {
-            lease.Dispose();
-        }
 
         GpuContextPool.DisposeAll();
     }
@@ -68,27 +52,17 @@ public class GpuContextPoolTests
         GpuContextPool.DisposeAll();
 
         Context firstContext;
-        GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
+        using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
         {
             firstContext = lease.Context;
         }
-        finally
-        {
-            lease.Dispose();
-        }
 
         GpuContextPool.DisposeAll();
 
-        lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
+        using (var lease = GpuContextPool.RentPreferred(preferCpu: true))
         {
             lease.Context.Should().NotBeSameAs(firstContext);
         }
-        finally
-        {
-            lease.Dispose();
-        }
 
         GpuContextPool.DisposeAll();
     }

PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs

@@ -199,17 +199,10 @@ public class GpuKernelCompilationTests
     [MemberData(nameof(KernelLoaders))]
     public void Kernel_compiles_without_internal_compiler_errors(string name, Action<Accelerator> compile)
     {
-        GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
-        {
-            var accelerator = lease.Accelerator;
-            Action action = () => compile(accelerator);
-            action.Should().NotThrow($"kernel {name} should compile");
-        }
-        finally
-        {
-            lease.Dispose();
-        }
+        using var lease = GpuContextPool.RentPreferred(preferCpu: true);
+        var accelerator = lease.Accelerator;
+        Action action = () => compile(accelerator);
+        action.Should().NotThrow($"kernel {name} should compile");
     }
 
     private static void CompileDivisorCycleCacheKernel(Accelerator accelerator)

PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs

@@ -150,36 +150,29 @@ public class NttGpuMathTests
         var modulus = new GpuUInt128(0UL, 17UL);
         var primitiveRoot = new GpuUInt128(0UL, 3UL);
 
-        GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
-        try
+        using var gpu = GpuContextPool.Rent();
+        var accelerator = gpu.Accelerator;
+
+        var method = typeof(NttGpuMath).GetMethod("GetSquareCache", BindingFlags.NonPublic | BindingFlags.Static);
+        method.Should().NotBeNull();
+        var cache = method!.Invoke(null, new object[] { accelerator, 8, modulus, primitiveRoot })!;
+        var cacheType = cache.GetType();
+
+        var twiddlesMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
+        var twiddlesInvMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesInvMont")!.GetValue(cache)!;
+        var r2 = (ulong)cacheType.GetProperty("MontR2Mod64")!.GetValue(cache)!;
+        var mod = (ulong)cacheType.GetProperty("ModulusLow")!.GetValue(cache)!;
+
+        int count = 7; // length 8 -> 7 twiddles
+        var forward = new GpuUInt128[count];
+        var inverse = new GpuUInt128[count];
+        twiddlesMont.CopyToCPU(forward);
+        twiddlesInvMont.CopyToCPU(inverse);
+
+        for (int i = 0; i < count; i++)
         {
-            var accelerator = gpu.Accelerator;
-
-            var method = typeof(NttGpuMath).GetMethod("GetSquareCache", BindingFlags.NonPublic | BindingFlags.Static);
-            method.Should().NotBeNull();
-            var cache = method!.Invoke(null, new object[] { accelerator, 8, modulus, primitiveRoot })!;
-            var cacheType = cache.GetType();
-
-            var twiddlesMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesMont")!.GetValue(cache)!;
-            var twiddlesInvMont = (MemoryBuffer1D<GpuUInt128, Stride1D.Dense>)cacheType.GetProperty("TwiddlesInvMont")!.GetValue(cache)!;
-            var r2 = (ulong)cacheType.GetProperty("MontR2Mod64")!.GetValue(cache)!;
-            var mod = (ulong)cacheType.GetProperty("ModulusLow")!.GetValue(cache)!;
-
-            int count = 7; // length 8 -> 7 twiddles
-            var forward = new GpuUInt128[count];
-            var inverse = new GpuUInt128[count];
-            twiddlesMont.CopyToCPU(forward);
-            twiddlesInvMont.CopyToCPU(inverse);
-
-            for (int i = 0; i < count; i++)
-            {
-                ulong prod = (ulong)((((BigInteger)forward[i].Low) * inverse[i].Low) % mod);
-                prod.Should().Be(r2);
-            }
-        }
-        finally
-        {
-            gpu.Dispose();
+            ulong prod = (ulong)((((BigInteger)forward[i].Low) * inverse[i].Low) % mod);
+            prod.Should().Be(r2);
         }
     }
 
@@ -345,23 +338,12 @@ public class NttGpuMathTests
         values[1] = new GpuUInt128(0UL, 2UL);
         var modulus = new GpuUInt128(0UL, 17UL);
         var primitiveRoot = new GpuUInt128(0UL, 3UL);
-        GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
-        try
+        using (var gpu = GpuContextPool.Rent())
         {
             var buffer = gpu.Accelerator.Allocate1D<GpuUInt128>(values.Length);
-            try
-            {
-                buffer.View.CopyFromCPU(ref values[0], values.Length);
-                NttGpuMath.SquareDevice(gpu.Accelerator, buffer.View, modulus, primitiveRoot);
-            }
-            finally
-            {
-                buffer.Dispose();
-            }
-        }
-        finally
-        {
-            gpu.Dispose();
+            buffer.View.CopyFromCPU(ref values[0], values.Length);
+            NttGpuMath.SquareDevice(gpu.Accelerator, buffer.View, modulus, primitiveRoot);
+            buffer.Dispose();
         }
 
         // Verify that the square cache now holds entries for the accelerator.

PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs

@@ -33,23 +33,16 @@ public class NttMontgomeryTests
         NttGpuMath.Square(cpu, modulus, primitiveRoot);
 
         // Device execution
-        GpuContextPool.GpuContextLease context = GpuContextPool.Rent();
-        try
-        {
-            var accelerator = context.Accelerator;
-            using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
-            buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
-            NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
-            buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
+        using var context = GpuContextPool.Rent();
+        var accelerator = context.Accelerator;
+        using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
+        buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
+        NttGpuMath.SquareDevice(accelerator, buffer.View, modulus, primitiveRoot);
+        buffer.View.CopyToCPU(ref gpu[0], gpu.Length);
 
-            for (int i = 0; i < cpu.Length; i++)
-            {
-                gpu[i].Should().Be(cpu[i]);
-            }
-        }
-        finally
+        for (int i = 0; i < cpu.Length; i++)
         {
-            context.Dispose();
+            gpu[i].Should().Be(cpu[i]);
         }
     }
 

PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs

@@ -11,21 +11,14 @@ public class GpuKernelPoolTests
     [Trait("Category", "Fast")]
     public void EnsureSmallCyclesOnDevice_returns_same_view_on_subsequent_calls()
     {
-        GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
-        {
-            var accelerator = lease.Accelerator;
+        using var lease = GpuContextPool.RentPreferred(preferCpu: true);
+        var accelerator = lease.Accelerator;
 
-            var first = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
-            var second = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
+        var first = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
+        var second = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
 
-            second.Equals(first).Should().BeTrue();
-            second.Length.Should().Be(MersenneDivisorCycles.SmallDivisorsMax + 1);
-        }
-        finally
-        {
-            lease.Dispose();
-        }
+        second.Equals(first).Should().BeTrue();
+        second.Length.Should().Be(MersenneDivisorCycles.SmallDivisorsMax + 1);
     }
 
     [Fact]

PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs

@@ -26,8 +26,7 @@ public class MersenneNumberDivisorByDivisorTesterTests
             markComposite: () => compositeMarks++,
             clearComposite: () => clearedMarks++,
             printResult: (prime, searched, detailed, passed) => results.Add((prime, searched, detailed, passed)),
-            threadCount: 1,
-            primeTestLimit: null);
+            threadCount: 1);
 
         tester.ConfiguredMaxPrime.Should().Be(13UL);
         compositeMarks.Should().Be(2);
@@ -60,7 +59,7 @@ public class MersenneNumberDivisorByDivisorTesterTests
             return DivisorLimit;
         }
 
-        public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
+        public bool IsPrime(ulong prime, out bool divisorsExhausted)
         {
             divisorsExhausted = true;
             return prime % 7UL != 0UL && prime % 11UL != 0UL;

PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs

@@ -7,88 +7,78 @@ namespace PerfectNumbers.Core.Tests;
 
 public class MersenneNumberDivisorCpuTesterTests
 {
-    [Fact]
-    [Trait("Category", "Fast")]
-    public void ByDivisor_tester_tracks_divisors_across_primes()
-    {
-        var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(11UL);
+        [Fact]
+        [Trait("Category", "Fast")]
+        public void ByDivisor_tester_tracks_divisors_across_primes()
+        {
+                var tester = new MersenneNumberDivisorByDivisorCpuTester();
+                tester.ConfigureFromMaxPrime(11UL);
 
-        tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue();
-        divisorsExhausted.Should().BeTrue();
+                tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue();
+                divisorsExhausted.Should().BeTrue();
 
-        tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue();
-        divisorsExhausted.Should().BeTrue();
+                tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue();
+                divisorsExhausted.Should().BeTrue();
 
-        tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse();
-        divisorsExhausted.Should().BeTrue();
-    }
+                tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse();
+                divisorsExhausted.Should().BeTrue();
+        }
 
-    [Fact]
-    [Trait("Category", "Fast")]
-    public void ByDivisor_prime_test_limit_disables_divisor_exhaustion()
-    {
-        var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(11UL);
+        [Fact]
+        [Trait("Category", "Fast")]
+        public void ByDivisor_session_checks_divisors_across_primes()
+        {
+                var tester = new MersenneNumberDivisorByDivisorCpuTester();
+                tester.ConfigureFromMaxPrime(13UL);
 
-        tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
-        divisorsExhausted.Should().BeFalse();
-    }
+                using var session = tester.CreateDivisorSession();
+                ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
+                byte[] hits = new byte[primes.Length];
 
-    [Fact]
-    [Trait("Category", "Fast")]
-    public void ByDivisor_session_checks_divisors_across_primes()
-    {
-        var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(13UL);
+                ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
+                session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
+                hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
 
-        using var session = tester.CreateDivisorSession();
-        ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
-        byte[] hits = new byte[primes.Length];
+                Array.Fill(hits, (byte)0);
+                ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
+                session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
+                hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
+        }
 
-        ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
-        session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
-        hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
+        [Fact]
+        [Trait("Category", "Fast")]
+        public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
+        {
+                var tester = new MersenneNumberDivisorByDivisorCpuTester();
+                tester.ConfigureFromMaxPrime(43UL);
 
-        Array.Fill(hits, (byte)0);
-        ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
-        session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
-        hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
-    }
+                using var session = tester.CreateDivisorSession();
+                ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
+                byte[] hits = new byte[exponents.Length];
 
-    [Fact]
-    [Trait("Category", "Fast")]
-    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
-    {
-        var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(43UL);
+                MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
+                ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
+                session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
 
-        using var session = tester.CreateDivisorSession();
-        ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
-        byte[] hits = new byte[exponents.Length];
+                hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+        }
 
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
-        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
-        session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
+        [Fact]
+        [Trait("Category", "Fast")]
+        public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
+        {
+                var tester = new MersenneNumberDivisorByDivisorCpuTester();
+                tester.ConfigureFromMaxPrime(61UL);
 
-        hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
-    }
+                using var session = tester.CreateDivisorSession();
+                ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
+                byte[] hits = new byte[exponents.Length];
 
-    [Fact]
-    [Trait("Category", "Fast")]
-    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
-    {
-        var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(61UL);
+                MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
+                ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
+                session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
 
-        using var session = tester.CreateDivisorSession();
-        ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
-        byte[] hits = new byte[exponents.Length];
-
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
-        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
-        session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
-
-        hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
-    }
+                hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+        }
 }
+

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -86,17 +86,6 @@ public class MersenneNumberDivisorGpuTesterTests
         exhausted.Should().BeTrue();
     }
 
-    [Fact]
-    [Trait("Category", "Fast")]
-    public void ByDivisor_gpu_prime_test_limit_stops_scan_immediately()
-    {
-        var tester = new MersenneNumberDivisorByDivisorGpuTester();
-        tester.ConfigureFromMaxPrime(43UL);
-
-        tester.IsPrime(31UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
-        divisorsExhausted.Should().BeFalse();
-    }
-
     [Fact]
     [Trait("Category", "Fast")]
     public void ByDivisor_tester_requires_configuration()

PerfectNumbers.Core.Tests/PerfectNumbers.Core.Tests.csproj

@@ -22,7 +22,4 @@
   <ItemGroup>
     <ProjectReference Include="../PerfectNumbers.Core/PerfectNumbers.Core.csproj" />
   </ItemGroup>
-  <ItemGroup>
-    <Compile Remove=".temp.cs" />
-  </ItemGroup>
 </Project>

PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs

@@ -7,9 +7,6 @@ namespace PerfectNumbers.Core.Tests;
 
 public class PrimeOrderCalculatorTests
 {
-    // EvenPerfectBitScanner never exercises the GPU heuristic device or the small-prime guards that the legacy tests cover,
-    // so those scenarios are commented out because they are unreachable in the production pipeline.
-    /*
     public static IEnumerable<object[]> SmallPrimes()
     {
         yield return new object[] { 3UL };
@@ -161,10 +158,9 @@ public class PrimeOrderCalculatorTests
         UInt128 prime = UInt128.Parse("18446744073709641691");
 
         UInt128 result = PrimeOrderCalculator.Calculate(
-                prime,
-                previousOrder: null,
-                PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
-                PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
+            prime,
+            previousOrder: null,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
 
         result.Should().Be((UInt128)1229782938247309446UL);
     }
@@ -200,5 +196,4 @@ public class PrimeOrderCalculatorTests
 
         return order;
     }
-*/
 }

PerfectNumbers.Core/ByDivisorDeltasDevice.cs

@@ -1,7 +0,0 @@
-namespace PerfectNumbers.Core;
-
-public enum ByDivisorDeltasDevice
-{
-    Cpu,
-    Gpu
-}

PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs

@@ -1,7 +0,0 @@
-namespace PerfectNumbers.Core;
-
-public enum ByDivisorMontgomeryDevice
-{
-    Cpu = 0,
-    Gpu = 1,
-}

PerfectNumbers.Core/DivisorCycleCache.cs

@@ -63,9 +63,7 @@ public sealed class DivisorCycleCache
         int[]? rentedMissing = null;
         Span<int> missingBuffer = length <= StackBufferThreshold
             ? stackalloc int[length]
-            : length < PerfectNumberConstants.PooledArrayThreshold
-                ? new Span<int>(new int[length])
-                : new Span<int>(rentedMissing = ThreadLocalArrayPool<int>.Shared.Rent(length), 0, length);
+            : new Span<int>(rentedMissing = ArrayPool<int>.Shared.Rent(length), 0, length);
 
         int missingCount = 0;
 
@@ -108,7 +106,7 @@ public sealed class DivisorCycleCache
 
         if (rentedMissing is not null)
         {
-            ThreadLocalArrayPool<int>.Shared.Return(rentedMissing, clearArray: false);
+            ArrayPool<int>.Shared.Return(rentedMissing, clearArray: false);
         }
     }
 
@@ -130,16 +128,12 @@ public sealed class DivisorCycleCache
         ulong[]? rentedDivisors = null;
         Span<ulong> workDivisors = count <= StackBufferThreshold
             ? stackalloc ulong[count]
-            : count < PerfectNumberConstants.PooledArrayThreshold
-                ? new Span<ulong>(new ulong[count])
-                : new Span<ulong>(rentedDivisors = ThreadLocalArrayPool<ulong>.Shared.Rent(count), 0, count);
+            : new Span<ulong>(rentedDivisors = ArrayPool<ulong>.Shared.Rent(count), 0, count);
 
         ulong[]? rentedResults = null;
         Span<ulong> workResults = count <= StackBufferThreshold
             ? stackalloc ulong[count]
-            : count < PerfectNumberConstants.PooledArrayThreshold
-                ? new Span<ulong>(new ulong[count])
-                : new Span<ulong>(rentedResults = ThreadLocalArrayPool<ulong>.Shared.Rent(count), 0, count);
+            : new Span<ulong>(rentedResults = ArrayPool<ulong>.Shared.Rent(count), 0, count);
 
         for (int i = 0; i < count; i++)
         {
@@ -155,12 +149,8 @@ public sealed class DivisorCycleCache
 
         if (rentedDivisors is not null)
         {
-            ThreadLocalArrayPool<ulong>.Shared.Return(rentedDivisors, clearArray: false);
-        }
-
-        if (rentedResults is not null)
-        {
-            ThreadLocalArrayPool<ulong>.Shared.Return(rentedResults, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(rentedDivisors, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(rentedResults!, clearArray: false);
         }
     }
 
@@ -182,30 +172,22 @@ public sealed class DivisorCycleCache
         ulong[]? rentedPow = null;
         Span<ulong> powSpan = length <= StackBufferThreshold
             ? stackalloc ulong[length]
-            : length < PerfectNumberConstants.PooledArrayThreshold
-                ? new Span<ulong>(new ulong[length])
-                : new Span<ulong>(rentedPow = ThreadLocalArrayPool<ulong>.Shared.Rent(length), 0, length);
+            : new Span<ulong>(rentedPow = ArrayPool<ulong>.Shared.Rent(length), 0, length);
 
         ulong[]? rentedOrder = null;
         Span<ulong> orderSpan = length <= StackBufferThreshold
             ? stackalloc ulong[length]
-            : length < PerfectNumberConstants.PooledArrayThreshold
-                ? new Span<ulong>(new ulong[length])
-                : new Span<ulong>(rentedOrder = ThreadLocalArrayPool<ulong>.Shared.Rent(length), 0, length);
+            : new Span<ulong>(rentedOrder = ArrayPool<ulong>.Shared.Rent(length), 0, length);
 
         ulong[]? rentedResult = null;
         Span<ulong> resultSpan = length <= StackBufferThreshold
             ? stackalloc ulong[length]
-            : length < PerfectNumberConstants.PooledArrayThreshold
-                ? new Span<ulong>(new ulong[length])
-                : new Span<ulong>(rentedResult = ThreadLocalArrayPool<ulong>.Shared.Rent(length), 0, length);
+            : new Span<ulong>(rentedResult = ArrayPool<ulong>.Shared.Rent(length), 0, length);
 
         byte[]? rentedStatus = null;
         Span<byte> statusSpan = length <= StackBufferThreshold
             ? stackalloc byte[length]
-            : length < PerfectNumberConstants.PooledArrayThreshold
-                ? new Span<byte>(new byte[length])
-                : new Span<byte>(rentedStatus = ThreadLocalArrayPool<byte>.Shared.Rent(length), 0, length);
+            : new Span<byte>(rentedStatus = ArrayPool<byte>.Shared.Rent(length), 0, length);
 
         ref ulong divisorRef = ref MemoryMarshal.GetReference(divisors);
         divisorBuffer.View.CopyFromCPU(ref divisorRef, length);
@@ -252,22 +234,10 @@ public sealed class DivisorCycleCache
         resultSpan.CopyTo(destination);
         if (rentedPow is not null)
         {
-            ThreadLocalArrayPool<ulong>.Shared.Return(rentedPow, clearArray: false);
-        }
-
-        if (rentedOrder is not null)
-        {
-            ThreadLocalArrayPool<ulong>.Shared.Return(rentedOrder, clearArray: false);
-        }
-
-        if (rentedResult is not null)
-        {
-            ThreadLocalArrayPool<ulong>.Shared.Return(rentedResult, clearArray: false);
-        }
-
-        if (rentedStatus is not null)
-        {
-            ThreadLocalArrayPool<byte>.Shared.Return(rentedStatus, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(rentedPow, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(rentedOrder!, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(rentedResult!, clearArray: false);
+            ArrayPool<byte>.Shared.Return(rentedStatus!, clearArray: false);
         }
 
         execution.Dispose();

PerfectNumbers.Core/ExponentRemainderStepper.cs

@@ -1,31 +1,20 @@
 namespace PerfectNumbers.Core;
 
-internal struct ExponentRemainderStepper
+internal struct ExponentRemainderStepper(in MontgomeryDivisorData divisor)
 {
-    private readonly MontgomeryDivisorData _divisor;
-    private readonly ulong _modulus;
-    private readonly ulong _nPrime;
-    private readonly ulong _montgomeryOne;
-    private ulong _previousExponent;
-    private ulong _currentMontgomery;
-    private bool _hasState;
+    private readonly MontgomeryDivisorData _divisor = divisor;
+    private readonly ulong _modulus = divisor.Modulus;
+    private readonly ulong _nPrime = divisor.NPrime;
+    private readonly ulong _montgomeryOne = divisor.MontgomeryOne;
+    private ulong _previousExponent = 0UL;
+    private ulong _currentMontgomery = divisor.MontgomeryOne;
+    private bool _hasState = false;
 
-    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 IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL;
 
-    public bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL;
+    public readonly bool HasState => _hasState;
 
-    public bool HasState => _hasState;
-
-    public ulong PreviousExponent => _previousExponent;
+    public readonly ulong PreviousExponent => _previousExponent;
 
     public void Reset()
     {

PerfectNumbers.Core/FactorEntry.cs

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

PerfectNumbers.Core/Gpu/GpuContextPool.cs

@@ -6,159 +6,156 @@ namespace PerfectNumbers.Core.Gpu;
 
 public static class GpuContextPool
 {
-    private static readonly bool PoolingEnabled = true;
-    // TODO: Introduce accelerator-specific warmup so pooled contexts precompile the ProcessEightBitWindows kernels and load
-    // divisor-cycle data before the first scan begins.
-    // Default device preference for generic GPU kernels (prime scans, NTT, etc.)
-    public static bool ForceCpu { get; set; } = false;
+        private static readonly bool PoolingEnabled = true;
+        // TODO: Introduce accelerator-specific warmup so pooled contexts precompile the ProcessEightBitWindows kernels and load
+        // divisor-cycle data before the first scan begins.
+        // Default device preference for generic GPU kernels (prime scans, NTT, etc.)
+        public static bool ForceCpu { get; set; } = false;
 
-    internal sealed class PooledContext
-    {
-        private bool disposed;
-
-        public Context Context { get; }
-        public Accelerator Accelerator { get; }
-        public bool IsCpu { get; }
-        public object ExecutionLock { get; } = new();
-
-        public PooledContext(bool preferCpu)
+        internal sealed class PooledContext
         {
-            Context = Context.CreateDefault();
-            Accelerator = Context.GetPreferredDevice(preferCpu).CreateAccelerator(Context);
-            IsCpu = Accelerator.AcceleratorType == AcceleratorType.CPU;
-            // NOTE: Avoid loading/compiling any kernel here to prevent implicit
-            // CL stream/queue creation during accelerator construction.
-            // Some OpenCL drivers are fragile when a queue is created immediately
-            // after device init. Real kernels will be JIT-loaded at first use.
+                private bool disposed;
+
+                public Context Context { get; }
+                public Accelerator Accelerator { get; }
+                public bool IsCpu { get; }
+                public object ExecutionLock { get; } = new();
+
+		public PooledContext(bool preferCpu)
+		{
+			Context = Context.CreateDefault();
+			Accelerator = Context.GetPreferredDevice(preferCpu).CreateAccelerator(Context);
+			IsCpu = Accelerator.AcceleratorType == AcceleratorType.CPU;
+			// NOTE: Avoid loading/compiling any kernel here to prevent implicit
+			// CL stream/queue creation during accelerator construction.
+			// Some OpenCL drivers are fragile when a queue is created immediately
+			// after device init. Real kernels will be JIT-loaded at first use.
+		}
+
+		public void Dispose()
+		{
+			if (disposed)
+			{
+				return;
+			}
+
+			// Ensure all cached GPU buffers for this accelerator are released.
+			NttGpuMath.ClearCaches(Accelerator);
+			// Release PrimeTester GPU state for this accelerator (kernel/device primes)
+			// Clear any per-accelerator cached resources to avoid releasing
+			// after the accelerator is destroyed.
+			PrimeTester.ClearGpuCaches(Accelerator);
+
+			Accelerator.Dispose();
+			Context.Dispose();
+			disposed = true;
+		}
+	}
+
+	private static readonly ConcurrentQueue<PooledContext> CpuPool = new();
+	private static readonly ConcurrentQueue<PooledContext> GpuPool = new();
+	private static readonly object CreationLock = new();
+
+	public static GpuContextLease Rent()
+	{
+		return RentPreferred(ForceCpu);
+	}
+
+	// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
+	public static GpuContextLease RentPreferred(bool preferCpu)
+	{
+		if (PoolingEnabled)
+		{
+			if (preferCpu && CpuPool.TryDequeue(out var cpu))
+			{
+				return new GpuContextLease(cpu, CreationLock);
+			}
+
+			if (!preferCpu && GpuPool.TryDequeue(out var gpu))
+			{
+				return new GpuContextLease(gpu, CreationLock);
+			}
+		}
+
+		// Serialize accelerator creation to avoid concurrent CL device/stream setup across threads.
+		lock (CreationLock)
+		{
+			return new GpuContextLease(new PooledContext(preferCpu), CreationLock);
+		}
+	}
+
+	public static void DisposeAll()
+	{
+		if (!PoolingEnabled)
+		{
+			return;
+		}
+
+		while (CpuPool.TryDequeue(out var cpu))
+		{
+			cpu.Dispose();
+		}
+
+		while (GpuPool.TryDequeue(out var gpu))
+		{
+			gpu.Dispose();
+		}
+	}
+
+	private static void Return(PooledContext ctx)
+	{
+		if (PoolingEnabled)
+		{
+			ctx.Accelerator.Synchronize();
+			if (ctx.IsCpu)
+			{
+				CpuPool.Enqueue(ctx);
+			}
+			else
+			{
+				GpuPool.Enqueue(ctx);
+			}
+		}
+
+                else
+                {
+                        lock (CreationLock)
+                        {
+                                ctx.Dispose();
+                        }
+                }
         }
 
-        public void Dispose()
-        {
-            if (disposed)
-            {
-                return;
-            }
-
-            // Ensure all cached GPU buffers for this accelerator are released.
-            NttGpuMath.ClearCaches(Accelerator);
-            // Release PrimeTester GPU state for this accelerator (kernel/device primes)
-            // Clear any per-accelerator cached resources to avoid releasing
-            // after the accelerator is destroyed.
-            PrimeTester.ClearGpuCaches(Accelerator);
-
-            Accelerator.Dispose();
-            Context.Dispose();
-            disposed = true;
-        }
-    }
-
-    private static readonly ConcurrentQueue<PooledContext> CpuPool = new();
-    private static readonly ConcurrentQueue<PooledContext> GpuPool = new();
-    private static readonly object CreationLock = new();
-
-    public static GpuContextLease Rent()
-    {
-        return RentPreferred(ForceCpu);
-    }
-
-    // Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
-    public static GpuContextLease RentPreferred(bool preferCpu)
-    {
-        if (PoolingEnabled)
-        {
-            if (preferCpu && CpuPool.TryDequeue(out var cpu))
-            {
-                return new GpuContextLease(cpu, CreationLock);
-            }
-
-            if (!preferCpu && GpuPool.TryDequeue(out var gpu))
-            {
-                return new GpuContextLease(gpu, CreationLock);
-            }
-        }
-
-        // Serialize accelerator creation to avoid concurrent CL device/stream setup across threads.
-        lock (CreationLock)
-        {
-            return new GpuContextLease(new PooledContext(preferCpu), CreationLock);
-        }
-    }
-
-    public static void DisposeAll()
-    {
-        if (!PoolingEnabled)
-        {
-            return;
-        }
-
-        while (CpuPool.TryDequeue(out var cpu))
-        {
-            cpu.Dispose();
-        }
-
-        while (GpuPool.TryDequeue(out var gpu))
-        {
-            gpu.Dispose();
-        }
-    }
-
-    private static void Return(PooledContext ctx)
-    {
-        if (PoolingEnabled)
-        {
-            ctx.Accelerator.Synchronize();
-            if (ctx.IsCpu)
-            {
-                CpuPool.Enqueue(ctx);
-            }
-            else
-            {
-                GpuPool.Enqueue(ctx);
-            }
-        }
-
-        else
-        {
-            lock (CreationLock)
-            {
-                ctx.Dispose();
-            }
-        }
-    }
-
-    public struct GpuContextLease
+    public struct GpuContextLease : IDisposable
     {
         private readonly PooledContext _ctx;
-        private bool disposed;
+                private bool disposed;
 
-        internal GpuContextLease(PooledContext ctx, object kernelInitLock)
+                internal GpuContextLease(PooledContext ctx, object kernelInitLock)
         {
             _ctx = ctx;
-            KernelInitLock = kernelInitLock;
-        }
+                        KernelInitLock = kernelInitLock;
+                }
 
         public Context Context => _ctx.Context;
 
         public Accelerator Accelerator => _ctx.Accelerator;
 
-        public object KernelInitLock { get; }
+                public object KernelInitLock { get; }
 
-        public object ExecutionLock => _ctx.ExecutionLock;
+                public object ExecutionLock => _ctx.ExecutionLock;
 
-        internal void ReturnToPool()
-        {
-            if (disposed)
-            {
-                return;
-            }
 
-            Return(_ctx);
-            disposed = true;
-        }
+                public void Dispose()
+                {
+                        if (disposed)
+			{
+				return;
+			}
 
-        public void Dispose()
-        {
-            ReturnToPool();
+			Return(_ctx);
+			disposed = true;
         }
     }
 }
+

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

@@ -37,7 +37,6 @@ public sealed class KernelContainer
             ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? IncrementalOrder;
     public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
         ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2ModOrder;
-    public Action<AcceleratorStream, Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int>? MontgomeryOddPowers;
 
     // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
     public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallCycles;
@@ -110,15 +109,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
         get
         {
             var accel = Accelerator; // avoid capturing 'this' in lambda
-            // lock (_gpu!.Value.KernelInitLock)
-            // {
+            lock (_gpu!.Value.KernelInitLock)
+            {
                 return KernelContainer.InitOnce(ref _kernels.Order, () =>
                 {
                     var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernelScan);
                     var kernel = KernelUtil.GetKernel(loaded);
                     return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>>();
                 });
-            // }
+            }
         }
     }
 
@@ -127,15 +126,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
         get
         {
             var accel = Accelerator;
-            // lock (_gpu!.Value.KernelInitLock)
-            // {
+            lock (_gpu!.Value.KernelInitLock)
+            {
                 return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () =>
                 {
                     var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernelScan);
                     var kernel = KernelUtil.GetKernel(loaded);
                     return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>>();
                 });
-            // }
+            }
         }
     }
 
@@ -144,15 +143,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
         get
         {
             var accel = Accelerator;
-            // lock (_gpu!.Value.KernelInitLock)
-            // {
+            lock (_gpu!.Value.KernelInitLock)
+            {
                 return KernelContainer.InitOnce(ref _kernels.Incremental, () =>
                 {
                     var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernelScan);
                     var kernel = KernelUtil.GetKernel(loaded);
                     return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>>();
                 });
-            // }
+            }
         }
     }
 
@@ -161,15 +160,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
         get
         {
             var accel = Accelerator;
-            // lock (_gpu!.Value.KernelInitLock)
-            // {
+            lock (_gpu!.Value.KernelInitLock)
+            {
                 return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () =>
                 {
                     var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalOrderKernelScan);
                     var kernel = KernelUtil.GetKernel(loaded);
                     return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
                 });
-            // }
+            }
         }
     }
 
@@ -178,66 +177,15 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
         get
         {
             var accel = Accelerator;
-            // lock (_gpu!.Value.KernelInitLock)
-            // {
+            lock (_gpu!.Value.KernelInitLock)
+            {
                 return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () =>
                 {
                     var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModOrderKernelScan);
                     var kernel = KernelUtil.GetKernel(loaded);
                     return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>>();
                 });
-            // }
-        }
-    }
-
-    public readonly Action<AcceleratorStream, Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int> MontgomeryOddPowerKernel
-    {
-        get
-        {
-            var accel = Accelerator;
-            // lock (_gpu!.Value.KernelInitLock)
-            // {
-                return KernelContainer.InitOnce(ref _kernels.MontgomeryOddPowers, () =>
-                {
-                    var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int>(MontgomeryOddPowerKernelScan);
-                    var kernel = KernelUtil.GetKernel(loaded);
-                    return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView1D<MontgomeryDivisorData, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int>>();
-                });
-            // }
-        }
-    }
-
-    private static void MontgomeryOddPowerKernelScan(
-        Index1D index,
-        ArrayView1D<MontgomeryDivisorData, Stride1D.Dense> divisors,
-        ArrayView1D<int, Stride1D.Dense> counts,
-        ArrayView1D<ulong, Stride1D.Dense> results,
-        int stride)
-    {
-        int tableIndex = index.X;
-        int baseOffset = tableIndex * stride;
-        MontgomeryDivisorData divisor = divisors[tableIndex];
-        int oddPowerCount = counts[tableIndex];
-        if (oddPowerCount <= 0)
-        {
-            return;
-        }
-
-        ulong montgomeryTwo = divisor.MontgomeryTwo;
-        results[baseOffset] = montgomeryTwo;
-        if (oddPowerCount == 1)
-        {
-            return;
-        }
-
-        ulong modulus = divisor.Modulus;
-        ulong nPrime = divisor.NPrime;
-        ulong square = divisor.MontgomeryTwoSquared;
-        ulong current = montgomeryTwo;
-        for (int i = 1; i < oddPowerCount; i++)
-        {
-            current = ULongExtensions.MontgomeryMultiply(current, square, modulus, nPrime);
-            results[baseOffset + i] = current;
+            }
         }
     }
 
@@ -396,10 +344,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
                 _stream?.Dispose();
                 _stream = null;
 
-                if (_gpu.HasValue)
-                {
-                    _gpu.Value.ReturnToPool();
-                }
+                _gpu?.Dispose();
                 _gpu = null;
 
                 _executionLock = null;
@@ -732,25 +677,25 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
         // GC.SuppressFinalize(this);
     }
 
-    public readonly struct ExecutionScope
+    public readonly struct ExecutionScope : IDisposable
     {
         private readonly object? _lock;
 
         public ExecutionScope(object? sync)
         {
             _lock = sync;
-            // if (_lock is not null)
-            // {
-            //     Monitor.Enter(_lock);
-            // }
+            if (_lock is not null)
+            {
+                Monitor.Enter(_lock);
+            }
         }
 
         public void Dispose()
         {
-            // if (_lock is not null)
-            // {
-            //     Monitor.Exit(_lock);
-            // }
+            if (_lock is not null)
+            {
+                Monitor.Exit(_lock);
+            }
         }
     }
 }
@@ -759,8 +704,6 @@ public class GpuKernelPool
 {
     private static readonly ConcurrentDictionary<Accelerator, KernelContainer> KernelCache = new(); // TODO: Replace this concurrent map with a simple accelerator-indexed lookup once kernel launchers are prewarmed during startup so we can drop the thread-safe wrapper entirely.
 
-    public static bool CpuCycleUsesGpu { get; set; }
-
     private static KernelContainer GetKernels(Accelerator accelerator)
     {
         return KernelCache.GetOrAdd(accelerator, _ => new KernelContainer());
@@ -777,8 +720,8 @@ public class GpuKernelPool
         }
 
         // Ensure single upload per accelerator even if multiple threads race here.
-        // lock (kernels) // TODO: Remove this lock by pre-uploading the immutable small-cycle snapshot during initialization; once no mutation happens at runtime, the pool must expose a simple reference without synchronization.
-        // {
+        lock (kernels) // TODO: Remove this lock by pre-uploading the immutable small-cycle snapshot during initialization; once no mutation happens at runtime, the pool must expose a simple reference without synchronization.
+        {
             if (kernels.SmallCycles is { } existing)
             {
                 return existing.View;
@@ -789,7 +732,7 @@ public class GpuKernelPool
             device.View.CopyFromCPU(host);
             kernels.SmallCycles = device;
             return device.View;
-        // }
+        }
     }
 
     public static ResiduePrimeViews EnsureSmallPrimesOnDevice(Accelerator accelerator)
@@ -803,8 +746,8 @@ public class GpuKernelPool
             return new ResiduePrimeViews(lastOne.View, lastSeven.View, lastOnePow2.View, lastSevenPow2.View);
         }
 
-        // lock (kernels) // TODO: Inline these small-prime uploads into startup initialization alongside the small-cycle snapshot so we can drop runtime locking and keep the GPU pool free of synchronization.
-        // {
+        lock (kernels) // TODO: Inline these small-prime uploads into startup initialization alongside the small-cycle snapshot so we can drop runtime locking and keep the GPU pool free of synchronization.
+        {
             if (kernels.SmallPrimesLastOne is { } existingLastOne &&
                 kernels.SmallPrimesLastSeven is { } existingLastSeven &&
                 kernels.SmallPrimesPow2LastOne is { } existingLastOnePow2 &&
@@ -833,7 +776,7 @@ public class GpuKernelPool
             kernels.SmallPrimesPow2LastSeven = deviceLastSevenPow2;
 
             return new ResiduePrimeViews(deviceLastOne.View, deviceLastSeven.View, deviceLastOnePow2.View, deviceLastSevenPow2.View);
-        // }
+        }
     }
 
     public static GpuKernelLease GetKernel(bool useGpuOrder)

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -7,6 +7,7 @@ using System.Numerics;
 using ILGPU;
 using ILGPU.Runtime;
 using MontgomeryDivisorData = PerfectNumbers.Core.MontgomeryDivisorData;
+using MontgomeryDivisorDataCache = PerfectNumbers.Core.MontgomeryDivisorDataCache;
 using CycleRemainderStepper = PerfectNumbers.Core.CycleRemainderStepper;
 
 namespace PerfectNumbers.Core.Gpu;
@@ -35,7 +36,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
 
     private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeExponent(Accelerator accelerator) =>
-        _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeExponentKernel));
+        _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeMontgomeryExponentKernel));
 
     private Action<Index1D, ArrayView<ulong>, ArrayView<byte>, byte> GetRemainderDeltaKernel(Accelerator accelerator) =>
         _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<byte>, byte>(ComputeRemainderDeltasKernel));
@@ -71,7 +72,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
     }
 
-    public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
+    public bool IsPrime(ulong prime, out bool divisorsExhausted)
     {
         ulong allowedMax;
         int batchCapacity;
@@ -108,7 +109,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             composite = CheckDivisors(
             prime,
             allowedMax,
-            timeLimit,
             accelerator,
             kernel,
             resources.DivisorsBuffer,
@@ -176,7 +176,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
     private static bool CheckDivisors(
         ulong prime,
         ulong allowedMax,
-        TimeSpan? timeLimit,
         Accelerator accelerator,
         Action<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>> kernel,
         MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> divisorsBuffer,
@@ -204,15 +203,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         bool processedAll = false;
         processedCount = 0UL;
         lastProcessed = 0UL;
-        coveredRange = false;
-
-        PrimeTestTimeLimit limitGuard;
-        if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard))
-        {
-            return false;
-        }
-
-        bool enforceLimit = limitGuard.IsActive;
 
         Span<MontgomeryDivisorData> divisorDataSpan;
         Span<ulong> divisorSpan;
@@ -288,12 +278,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         {
             while (currentK <= maxK)
             {
-                if (enforceLimit && limitGuard.HasExpired())
-                {
-                    processedAll = false;
-                    return false;
-                }
-
                 int chunkCount = Math.Min(chunkCapacity, batchCapacity);
                 ulong remainingK = maxK - currentK + 1UL;
                 if ((ulong)chunkCount > remainingK)
@@ -412,8 +396,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                     ulong divisorValue = filteredDivisorsSpan[i];
                     divisorSpan[i] = divisorValue;
                     exponentSpan[i] = prime;
-                    // Divisors never repeat while scanning the monotonically increasing sequence, so regenerate
-                    // the Montgomery data for each entry rather than storing cache lines that no future batch would hit.
                     divisorDataSpan[i] = MontgomeryDivisorData.FromModulus(divisorValue);
                 }
 
@@ -569,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulong exponent = exponents[index];
-        hits[index] = exponent.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -777,7 +759,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         mask[globalIndex] = 1;
     }
 
-    private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
+    private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
     {
         ulong modulus = divisor.Modulus;
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -788,7 +770,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         ulong exponent = exponents[index];
 
-        results[index] = exponent.Pow2ModBinaryGpu(modulus);
+        results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
     }
 
     public ulong DivisorLimit

PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs

@@ -1,262 +0,0 @@
-using System;
-using System.Buffers;
-using System.Runtime.InteropServices;
-using ILGPU;
-using ILGPU.Runtime;
-
-namespace PerfectNumbers.Core.Gpu;
-
-internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable
-{
-    private const int RemainderCount = 6;
-
-    private readonly GpuContextPool.GpuContextLease _lease;
-    private readonly Accelerator _accelerator;
-    private readonly Action<Index1D, CandidateFilterArgs, ArrayView<ulong>, ArrayView<byte>> _kernel;
-    private MemoryBuffer1D<ulong, Stride1D.Dense> _candidateBuffer;
-    private MemoryBuffer1D<byte, Stride1D.Dense> _maskBuffer;
-    private ulong[] _hostCandidates;
-    private byte[] _hostMask;
-    private int _capacity;
-    private bool _disposed;
-
-    internal MersenneNumberDivisorCandidateGpuEvaluator(int batchSize)
-    {
-        _lease = GpuContextPool.RentPreferred(preferCpu: false);
-        _accelerator = _lease.Accelerator;
-        _kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, CandidateFilterArgs, ArrayView<ulong>, ArrayView<byte>>(EvaluateCandidatesKernel);
-        _capacity = Math.Max(1, batchSize);
-        _candidateBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _maskBuffer = _accelerator.Allocate1D<byte>(_capacity);
-        _hostCandidates = ArrayPool<ulong>.Shared.Rent(_capacity);
-        _hostMask = ArrayPool<byte>.Shared.Rent(_capacity);
-    }
-
-    public void Dispose()
-    {
-        if (_disposed)
-        {
-            return;
-        }
-
-        _disposed = true;
-        _candidateBuffer.Dispose();
-        _maskBuffer.Dispose();
-        ArrayPool<ulong>.Shared.Return(_hostCandidates, clearArray: false);
-        ArrayPool<byte>.Shared.Return(_hostMask, clearArray: false);
-        _lease.Dispose();
-    }
-
-    internal void EvaluateCandidates(
-        ulong startDivisor,
-        ulong step,
-        ulong limit,
-        ReadOnlySpan<byte> remainders,
-        ReadOnlySpan<byte> steps,
-        bool lastIsSeven,
-        Span<ulong> candidates,
-        Span<byte> mask)
-    {
-        if (_disposed)
-        {
-            throw new ObjectDisposedException(nameof(MersenneNumberDivisorCandidateGpuEvaluator));
-        }
-
-        if (remainders.Length != RemainderCount || steps.Length != RemainderCount)
-        {
-            throw new ArgumentException("Remainder and step tables must contain exactly six entries.");
-        }
-
-        int count = candidates.Length;
-        if (count == 0)
-        {
-            return;
-        }
-
-        EnsureCapacity(count);
-
-        Span<ulong> hostCandidateSpan = _hostCandidates.AsSpan(0, count);
-        Span<byte> hostMaskSpan = _hostMask.AsSpan(0, count);
-
-        ArrayView1D<ulong, Stride1D.Dense> candidateView = _candidateBuffer.View.SubView(0, count);
-        ArrayView1D<byte, Stride1D.Dense> maskView = _maskBuffer.View.SubView(0, count);
-
-        var args = new CandidateFilterArgs(
-            startDivisor,
-            step,
-            limit,
-            remainders[0],
-            remainders[1],
-            remainders[2],
-            remainders[3],
-            remainders[4],
-            remainders[5],
-            steps[0],
-            steps[1],
-            steps[2],
-            steps[3],
-            steps[4],
-            steps[5],
-            lastIsSeven ? (byte)1 : (byte)0);
-
-        _kernel(count, args, candidateView, maskView);
-
-        candidateView.CopyToCPU(ref MemoryMarshal.GetReference(hostCandidateSpan), count);
-        maskView.CopyToCPU(ref MemoryMarshal.GetReference(hostMaskSpan), count);
-
-        hostCandidateSpan.CopyTo(candidates);
-        hostMaskSpan.CopyTo(mask);
-    }
-
-    private void EnsureCapacity(int required)
-    {
-        if (required <= _capacity)
-        {
-            return;
-        }
-
-        int newCapacity = _capacity == 0 ? 1 : _capacity;
-        while (newCapacity < required)
-        {
-            newCapacity <<= 1;
-        }
-
-        _candidateBuffer.Dispose();
-        _maskBuffer.Dispose();
-        ArrayPool<ulong>.Shared.Return(_hostCandidates, clearArray: false);
-        ArrayPool<byte>.Shared.Return(_hostMask, clearArray: false);
-
-        _capacity = newCapacity;
-        _candidateBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _maskBuffer = _accelerator.Allocate1D<byte>(_capacity);
-        _hostCandidates = ArrayPool<ulong>.Shared.Rent(_capacity);
-        _hostMask = ArrayPool<byte>.Shared.Rent(_capacity);
-    }
-
-    private static void EvaluateCandidatesKernel(
-        Index1D index,
-        CandidateFilterArgs args,
-        ArrayView<ulong> candidates,
-        ArrayView<byte> mask)
-    {
-        int idx = index;
-        if (idx >= mask.Length)
-        {
-            return;
-        }
-
-        ulong offset = (ulong)idx;
-        ulong candidate = args.StartDivisor + args.Step * offset;
-        candidates[idx] = candidate;
-
-        if (candidate > args.Limit || (args.Step != 0UL && candidate < args.StartDivisor))
-        {
-            mask[idx] = 0;
-            return;
-        }
-
-        int remainderValue10 = ComputeRemainder(args.Remainder10, args.Step10, idx, 10);
-        int remainderValue8 = ComputeRemainder(args.Remainder8, args.Step8, idx, 8);
-        int remainderValue5 = ComputeRemainder(args.Remainder5, args.Step5, idx, 5);
-        int remainderValue3 = ComputeRemainder(args.Remainder3, args.Step3, idx, 3);
-        int remainderValue7 = ComputeRemainder(args.Remainder7, args.Step7, idx, 7);
-        int remainderValue11 = ComputeRemainder(args.Remainder11, args.Step11, idx, 11);
-
-        bool lastIsSeven = args.LastIsSevenFlag != 0;
-        bool admissible = lastIsSeven
-            ? (remainderValue10 == 3 || remainderValue10 == 7 || remainderValue10 == 9)
-            : (remainderValue10 == 1 || remainderValue10 == 3 || remainderValue10 == 9);
-
-        if (!admissible)
-        {
-            mask[idx] = 0;
-            return;
-        }
-
-        if ((remainderValue8 != 1 && remainderValue8 != 7)
-            || remainderValue3 == 0
-            || remainderValue5 == 0
-            || remainderValue7 == 0
-            || remainderValue11 == 0)
-        {
-            mask[idx] = 0;
-            return;
-        }
-
-        mask[idx] = 1;
-    }
-
-    private static int ComputeRemainder(byte remainder, byte step, int index, int modulus)
-    {
-        if (modulus == 0)
-        {
-            return 0;
-        }
-
-        long value = remainder;
-        value += (long)step * index;
-        int result = (int)(value % modulus);
-        if (result < 0)
-        {
-            result += modulus;
-        }
-
-        return result;
-    }
-
-    private readonly struct CandidateFilterArgs
-    {
-        public CandidateFilterArgs(
-            ulong startDivisor,
-            ulong step,
-            ulong limit,
-            byte remainder10,
-            byte remainder8,
-            byte remainder5,
-            byte remainder3,
-            byte remainder7,
-            byte remainder11,
-            byte step10,
-            byte step8,
-            byte step5,
-            byte step3,
-            byte step7,
-            byte step11,
-            byte lastIsSevenFlag)
-        {
-            StartDivisor = startDivisor;
-            Step = step;
-            Limit = limit;
-            Remainder10 = remainder10;
-            Remainder8 = remainder8;
-            Remainder5 = remainder5;
-            Remainder3 = remainder3;
-            Remainder7 = remainder7;
-            Remainder11 = remainder11;
-            Step10 = step10;
-            Step8 = step8;
-            Step5 = step5;
-            Step3 = step3;
-            Step7 = step7;
-            Step11 = step11;
-            LastIsSevenFlag = lastIsSevenFlag;
-        }
-
-        public ulong StartDivisor { get; }
-        public ulong Step { get; }
-        public ulong Limit { get; }
-        public byte Remainder10 { get; }
-        public byte Remainder8 { get; }
-        public byte Remainder5 { get; }
-        public byte Remainder3 { get; }
-        public byte Remainder7 { get; }
-        public byte Remainder11 { get; }
-        public byte Step10 { get; }
-        public byte Step8 { get; }
-        public byte Step5 { get; }
-        public byte Step3 { get; }
-        public byte Step7 { get; }
-        public byte Step11 { get; }
-        public byte LastIsSevenFlag { get; }
-    }
-}

PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs

@@ -1,179 +0,0 @@
-using System;
-using System.Buffers;
-using System.Runtime.InteropServices;
-using PerfectNumbers.Core;
-using ILGPU;
-using ILGPU.Runtime;
-
-namespace PerfectNumbers.Core.Gpu;
-
-internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
-{
-    private readonly GpuContextPool.GpuContextLease _lease;
-    private readonly Accelerator _accelerator;
-    private readonly Action<Index1D, ArrayView<ulong>, ArrayView<MontgomeryDivisorData>> _kernel;
-    private MemoryBuffer1D<ulong, Stride1D.Dense> _modulusBuffer;
-    private MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> _resultBuffer;
-    private ulong[] _hostModulusBuffer;
-    private MontgomeryDivisorData[] _hostResultBuffer;
-    private int _capacity;
-    private bool _disposed;
-
-    internal MersenneNumberDivisorMontgomeryGpuBuilder(int batchSize)
-    {
-        _lease = GpuContextPool.RentPreferred(preferCpu: false);
-        _accelerator = _lease.Accelerator;
-        _kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<MontgomeryDivisorData>>(BuildKernel);
-        _capacity = Math.Max(1, batchSize);
-        _modulusBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _resultBuffer = _accelerator.Allocate1D<MontgomeryDivisorData>(_capacity);
-        _hostModulusBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
-        _hostResultBuffer = ArrayPool<MontgomeryDivisorData>.Shared.Rent(_capacity);
-    }
-
-    public void Dispose()
-    {
-        if (_disposed)
-        {
-            return;
-        }
-
-        _disposed = true;
-        _modulusBuffer.Dispose();
-        _resultBuffer.Dispose();
-        ArrayPool<ulong>.Shared.Return(_hostModulusBuffer, clearArray: false);
-        ArrayPool<MontgomeryDivisorData>.Shared.Return(_hostResultBuffer, clearArray: false);
-        _lease.Dispose();
-    }
-
-    internal void Build(ReadOnlySpan<ulong> moduli, Span<MontgomeryDivisorData> destination)
-    {
-        if (_disposed)
-        {
-            throw new ObjectDisposedException(nameof(MersenneNumberDivisorMontgomeryGpuBuilder));
-        }
-
-        int length = moduli.Length;
-        if (length == 0)
-        {
-            return;
-        }
-
-        EnsureCapacity(length);
-
-        Span<ulong> hostModulusSpan = _hostModulusBuffer.AsSpan(0, length);
-        moduli.CopyTo(hostModulusSpan);
-
-        ArrayView1D<ulong, Stride1D.Dense> modulusView = _modulusBuffer.View.SubView(0, length);
-        modulusView.CopyFromCPU(ref MemoryMarshal.GetReference(hostModulusSpan), length);
-
-        ArrayView1D<MontgomeryDivisorData, Stride1D.Dense> resultView = _resultBuffer.View.SubView(0, length);
-        _kernel(length, modulusView, resultView);
-
-        Span<MontgomeryDivisorData> hostResultSpan = _hostResultBuffer.AsSpan(0, length);
-        resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostResultSpan), length);
-
-        hostResultSpan.CopyTo(destination);
-    }
-
-    private void EnsureCapacity(int required)
-    {
-        if (required <= _capacity)
-        {
-            return;
-        }
-
-        int newCapacity = _capacity == 0 ? 1 : _capacity;
-        while (newCapacity < required)
-        {
-            newCapacity <<= 1;
-        }
-
-        _modulusBuffer.Dispose();
-        _resultBuffer.Dispose();
-        ArrayPool<ulong>.Shared.Return(_hostModulusBuffer, clearArray: false);
-        ArrayPool<MontgomeryDivisorData>.Shared.Return(_hostResultBuffer, clearArray: false);
-
-        _capacity = newCapacity;
-        _modulusBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _resultBuffer = _accelerator.Allocate1D<MontgomeryDivisorData>(_capacity);
-        _hostModulusBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
-        _hostResultBuffer = ArrayPool<MontgomeryDivisorData>.Shared.Rent(_capacity);
-    }
-
-    private static void BuildKernel(Index1D index, ArrayView<ulong> moduli, ArrayView<MontgomeryDivisorData> destination)
-    {
-        int idx = index;
-        if (idx >= moduli.Length)
-        {
-            return;
-        }
-
-        ulong modulus = moduli[idx];
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            destination[idx] = default;
-            return;
-        }
-
-        ulong nPrime = ComputeMontgomeryNPrime(modulus);
-        ulong montgomeryOne = ComputeMontgomeryResidue64(modulus);
-        ulong montgomeryTwo = montgomeryOne.MulMod64GpuCompatible(2UL, modulus);
-        ulong montgomeryTwoSquared = MontgomeryMultiplyGpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
-
-        destination[idx] = new MontgomeryDivisorData(
-            modulus,
-            nPrime,
-            montgomeryOne,
-            montgomeryTwo,
-            montgomeryTwoSquared);
-    }
-
-    private static ulong ComputeMontgomeryResidue64(ulong modulus)
-    {
-        ulong remainder = ulong.MaxValue % modulus;
-        remainder++;
-        if (remainder >= modulus)
-        {
-            remainder -= modulus;
-        }
-
-        return remainder;
-    }
-
-    private static ulong MontgomeryMultiplyGpu(ulong a, ulong b, ulong modulus, ulong nPrime)
-    {
-        GpuUInt128 product = new(a);
-        product.Mul64(new GpuUInt128(b));
-        ulong tLow = product.Low;
-        ulong tHigh = product.High;
-
-        ulong m = unchecked(tLow * nPrime);
-        GpuUInt128 mProduct = new(m);
-        mProduct.Mul64(new GpuUInt128(modulus));
-        ulong mLow = mProduct.Low;
-        ulong mHigh = mProduct.High;
-
-        ulong carry = unchecked((tLow + mLow) < tLow ? 1UL : 0UL);
-        ulong result = tHigh + mHigh + carry;
-        if (result >= modulus)
-        {
-            result -= modulus;
-        }
-
-        return result;
-    }
-
-    private static ulong ComputeMontgomeryNPrime(ulong modulus)
-    {
-        ulong inv = modulus;
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        return unchecked(0UL - inv);
-    }
-}

PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs

@@ -1,140 +0,0 @@
-using System;
-using System.Buffers;
-using System.Runtime.InteropServices;
-using ILGPU;
-using ILGPU.Runtime;
-
-namespace PerfectNumbers.Core.Gpu;
-
-internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable
-{
-    private const int RemainderCount = 6;
-
-    private readonly GpuContextPool.GpuContextLease _lease;
-    private readonly Accelerator _accelerator;
-    private readonly Action<Index1D, int, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>> _kernel;
-    private MemoryBuffer1D<byte, Stride1D.Dense> _remainderBuffer;
-    private MemoryBuffer1D<byte, Stride1D.Dense> _stepBuffer;
-    private MemoryBuffer1D<byte, Stride1D.Dense> _resultBuffer;
-    private byte[] _hostRemainders;
-    private byte[] _hostSteps;
-    private readonly int _length;
-    private int _capacity;
-    private bool _disposed;
-
-    internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan<byte> steps)
-    {
-        if (steps.Length != RemainderCount)
-        {
-            throw new ArgumentException("Step table must contain exactly six entries.");
-        }
-
-        _lease = GpuContextPool.RentPreferred(preferCpu: false);
-        _accelerator = _lease.Accelerator;
-        _kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, int, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>>(AdvanceRemaindersKernel);
-
-        _length = steps.Length;
-        _capacity = Math.Max(1, _length);
-        _remainderBuffer = _accelerator.Allocate1D<byte>(_capacity);
-        _stepBuffer = _accelerator.Allocate1D<byte>(_capacity);
-        _resultBuffer = _accelerator.Allocate1D<byte>(_capacity);
-
-        _hostRemainders = ArrayPool<byte>.Shared.Rent(_capacity);
-        _hostSteps = ArrayPool<byte>.Shared.Rent(_capacity);
-
-        steps.CopyTo(_hostSteps.AsSpan(0, _length));
-
-        ArrayView1D<byte, Stride1D.Dense> stepView = _stepBuffer.View.SubView(0, _length);
-        stepView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostSteps.AsSpan(0, _length)), _length);
-    }
-
-    public void Dispose()
-    {
-        if (_disposed)
-        {
-            return;
-        }
-
-        _disposed = true;
-        _remainderBuffer.Dispose();
-        _stepBuffer.Dispose();
-        _resultBuffer.Dispose();
-        ArrayPool<byte>.Shared.Return(_hostRemainders, clearArray: false);
-        ArrayPool<byte>.Shared.Return(_hostSteps, clearArray: false);
-        _lease.Dispose();
-    }
-
-    internal void Advance(int advanceCount, Span<byte> remainders)
-    {
-        if (_disposed)
-        {
-            throw new ObjectDisposedException(nameof(MersenneNumberDivisorRemainderGpuStepper));
-        }
-
-        if (advanceCount <= 0 || remainders.Length == 0)
-        {
-            return;
-        }
-
-        if (remainders.Length != _length)
-        {
-            throw new ArgumentException("Remainder buffer length must match the configured length.");
-        }
-
-        Span<byte> hostRemainderSpan = _hostRemainders.AsSpan(0, _length);
-        remainders.CopyTo(hostRemainderSpan);
-
-        ArrayView1D<byte, Stride1D.Dense> remainderView = _remainderBuffer.View.SubView(0, _length);
-        remainderView.CopyFromCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length);
-
-        ArrayView1D<byte, Stride1D.Dense> stepView = _stepBuffer.View.SubView(0, _length);
-        ArrayView1D<byte, Stride1D.Dense> resultView = _resultBuffer.View.SubView(0, _length);
-
-        _kernel(_length, advanceCount, remainderView, stepView, resultView);
-
-        resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length);
-        hostRemainderSpan.CopyTo(remainders);
-    }
-
-    private static void AdvanceRemaindersKernel(Index1D index, int advanceCount, ArrayView<byte> remainders, ArrayView<byte> steps, ArrayView<byte> results)
-    {
-        int idx = index;
-        if (idx >= remainders.Length)
-        {
-            return;
-        }
-
-        int modulus = GetModulus(idx);
-        if (modulus == 0)
-        {
-            results[idx] = remainders[idx];
-            return;
-        }
-
-        int remainder = remainders[idx] % modulus;
-        int step = steps[idx] % modulus;
-
-        long delta = (long)step * advanceCount;
-        int adjusted = (int)((remainder + delta) % modulus);
-        if (adjusted < 0)
-        {
-            adjusted += modulus;
-        }
-
-        results[idx] = (byte)adjusted;
-    }
-
-    private static int GetModulus(int index)
-    {
-        return index switch
-        {
-            0 => 10,
-            1 => 8,
-            2 => 5,
-            3 => 3,
-            4 => 7,
-            5 => 11,
-            _ => 0,
-        };
-    }
-}

PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs

@@ -1,121 +0,0 @@
-using System;
-using System.Buffers;
-using System.Runtime.InteropServices;
-using PerfectNumbers.Core;
-using ILGPU;
-using ILGPU.Runtime;
-
-namespace PerfectNumbers.Core.Gpu;
-
-internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
-{
-    private readonly GpuContextPool.GpuContextLease _lease;
-    private readonly Accelerator _accelerator;
-    private readonly Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel;
-    private MemoryBuffer1D<ulong, Stride1D.Dense> _exponentBuffer;
-    private MemoryBuffer1D<ulong, Stride1D.Dense> _resultBuffer;
-    private ulong[] _hostBuffer;
-    private int _capacity;
-    private bool _disposed;
-
-    internal MersenneNumberDivisorResidueGpuEvaluator(int batchSize)
-    {
-        _lease = GpuContextPool.RentPreferred(preferCpu: false);
-        _accelerator = _lease.Accelerator;
-        _kernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeResiduesKernel);
-        _capacity = Math.Max(1, batchSize);
-        _exponentBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _resultBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
-    }
-
-    public void Dispose()
-    {
-        if (_disposed)
-        {
-            return;
-        }
-
-        _disposed = true;
-        _exponentBuffer.Dispose();
-        _resultBuffer.Dispose();
-        ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
-        _lease.Dispose();
-    }
-
-    internal void ComputeResidues(ReadOnlySpan<ulong> exponents, in MontgomeryDivisorData divisorData, Span<ulong> destination)
-    {
-        if (_disposed)
-        {
-            throw new ObjectDisposedException(nameof(MersenneNumberDivisorResidueGpuEvaluator));
-        }
-
-        int length = exponents.Length;
-        if (length == 0)
-        {
-            return;
-        }
-
-        EnsureCapacity(length);
-
-        Span<ulong> hostSpan = _hostBuffer.AsSpan(0, length);
-        exponents.CopyTo(hostSpan);
-
-        ArrayView1D<ulong, Stride1D.Dense> exponentView = _exponentBuffer.View.SubView(0, length);
-        exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length);
-
-        ArrayView1D<ulong, Stride1D.Dense> resultView = _resultBuffer.View.SubView(0, length);
-        _kernel(length, divisorData, exponentView, resultView);
-        resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostSpan), length);
-
-        hostSpan.CopyTo(destination);
-    }
-
-    private void EnsureCapacity(int required)
-    {
-        if (required <= _capacity)
-        {
-            return;
-        }
-
-        int newCapacity = _capacity;
-        if (newCapacity == 0)
-        {
-            newCapacity = 1;
-        }
-
-        while (newCapacity < required)
-        {
-            newCapacity <<= 1;
-        }
-
-        _exponentBuffer.Dispose();
-        _resultBuffer.Dispose();
-        ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
-
-        _capacity = newCapacity;
-        _exponentBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _resultBuffer = _accelerator.Allocate1D<ulong>(_capacity);
-        _hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
-    }
-
-    private static void ComputeResiduesKernel(Index1D index, MontgomeryDivisorData divisorData, ArrayView<ulong> exponents, ArrayView<ulong> results)
-    {
-        int idx = index;
-        int length = (int)exponents.Length;
-        if (idx >= length)
-        {
-            return;
-        }
-
-        ulong exponent = exponents[idx];
-        ulong modulus = divisorData.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            results[idx] = 0UL;
-            return;
-        }
-
-        results[idx] = exponent.Pow2ModBinaryGpu(modulus);
-    }
-}

PerfectNumbers.Core/Gpu/MontgomeryOddPowerGpu.cs

@@ -1,94 +0,0 @@
-using System;
-using System.Runtime.InteropServices;
-using ILGPU.Runtime;
-
-namespace PerfectNumbers.Core.Gpu;
-
-internal static class MontgomeryOddPowerGpu
-{
-    public static bool TryCompute(in MontgomeryDivisorData divisor, int oddPowerCount, Span<ulong> destination)
-    {
-        if (oddPowerCount <= 0)
-        {
-            return true;
-        }
-
-        Span<MontgomeryDivisorData> divisorSpan = stackalloc MontgomeryDivisorData[1];
-        divisorSpan[0] = divisor;
-        Span<int> countSpan = stackalloc int[1];
-        countSpan[0] = oddPowerCount;
-        Span<ulong> resultSpan = stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount];
-
-        if (!TryComputeBatch(divisorSpan, countSpan, resultSpan, PerfectNumberConstants.MaxOddPowersCount))
-        {
-            return false;
-        }
-
-        resultSpan[..oddPowerCount].CopyTo(destination);
-        return true;
-    }
-
-    public static bool TryComputeBatch(
-        ReadOnlySpan<MontgomeryDivisorData> divisors,
-        ReadOnlySpan<int> oddPowerCounts,
-        Span<ulong> results,
-        int stride)
-    {
-        int divisorCount = divisors.Length;
-        if (divisorCount == 0)
-        {
-            return true;
-        }
-
-        if (oddPowerCounts.Length < divisorCount)
-        {
-            throw new ArgumentException(
-                "oddPowerCounts span must contain at least as many elements as divisors span.",
-                nameof(oddPowerCounts));
-        }
-
-        if (results.Length < divisorCount * stride)
-        {
-            throw new ArgumentException(
-                "The result span is shorter than the required batched odd power storage.",
-                nameof(results));
-        }
-
-        if (GpuContextPool.ForceCpu)
-        {
-            return false;
-        }
-
-        var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
-        var execution = lease.EnterExecutionScope();
-        var accelerator = lease.Accelerator;
-        var stream = lease.Stream;
-
-        using var divisorBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(divisorCount);
-        using var countBuffer = accelerator.Allocate1D<int>(divisorCount);
-        using var resultBuffer = accelerator.Allocate1D<ulong>(divisorCount * stride);
-
-        try
-        {
-            divisorBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(divisors), divisorCount);
-            countBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(oddPowerCounts), divisorCount);
-            resultBuffer.MemSetToZero();
-
-            var kernel = lease.MontgomeryOddPowerKernel;
-            kernel(stream, divisorCount, divisorBuffer.View, countBuffer.View, resultBuffer.View, stride);
-            stream.Synchronize();
-
-            resultBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(results), divisorCount * stride);
-            return true;
-        }
-        catch
-        {
-            return false;
-        }
-        finally
-        {
-            execution.Dispose();
-            lease.Dispose();
-        }
-    }
-}

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -811,7 +811,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             ulong reduced = phi / factor;
-            if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL)
+            if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL)
             {
                 return false;
             }
@@ -1595,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
     {
-        return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL;
+        return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
     }
 
     private static ulong CalculateByDoublingKernel(ulong prime)
@@ -1715,10 +1715,10 @@ internal static partial class PrimeOrderGpuHeuristics
         GpuUInt128[]? rentedExponents = null;
         GpuUInt128[]? rentedResults = null;
         var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
-        var execution = lease.EnterExecutionScope();
 
         try
         {
+            var execution = lease.EnterExecutionScope();
             Accelerator accelerator = lease.Accelerator;
             AcceleratorStream stream = lease.Stream;
             var kernel = GetPow2ModWideKernel(accelerator);
@@ -1755,6 +1755,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
             exponentBuffer.Dispose();
             remainderBuffer.Dispose();
+            execution.Dispose();
             return true;
         }
         catch (Exception)
@@ -1774,7 +1775,6 @@ internal static partial class PrimeOrderGpuHeuristics
                 ArrayPool<GpuUInt128>.Shared.Return(rentedResults, clearArray: false);
             }
 
-            execution.Dispose();
             lease.Dispose();
         }
     }
@@ -1811,57 +1811,7 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong exponent = exponents[index];
-        ulong modulus = divisor.Modulus;
-		// ulong nPrime = divisor.NPrime;
-		// ulong result = divisor.MontgomeryOne;
-		// ulong baseValue = divisor.MontgomeryTwo;
-		// ulong remainingExponent = exponent;
-
-		// if (remainingExponent == 0UL)
-		// {
-		//     remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
-		//     return;
-		// }
-
-		// while (remainingExponent != 0UL)
-		// {
-		//     if ((remainingExponent & 1UL) != 0UL)
-		//     {
-		//         result = result.MontgomeryMultiply(baseValue, modulus, nPrime);
-		//     }
-
-		//     remainingExponent >>= 1;
-		//     if (remainingExponent == 0UL)
-		//     {
-		//         break;
-		//     }
-
-		//     // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term.
-		//     baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime);
-		// }
-
-		// remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
-
-		// Plain modular exponentiation without Montgomery reduction kept for reference.
-		ulong plainResult = 1UL;// % modulus;
-		ulong plainBase = 2UL;// % modulus;
-		while (exponent != 0UL)
-		{
-			if ((exponent & 1UL) != 0UL)
-			{
-				plainResult = (plainResult * plainBase) % modulus;
-			}
-
-			exponent >>= 1;
-			if (exponent == 0UL)
-			{
-				break;
-			}
-
-			plainBase = (plainBase * plainBase) % modulus;
-		}
-		
-        remainders[index] = plainResult;
+        remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)

PerfectNumbers.Core/HeuristicFailureLog.cs

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

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

PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs

@@ -12,7 +12,7 @@ public interface IMersenneNumberDivisorByDivisorTester
 
     ulong GetAllowedMaxDivisor(ulong prime);
 
-    bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null);
+    bool IsPrime(ulong prime, out bool divisorsExhausted);
 
     IDivisorScanSession CreateDivisorSession();
 

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -1,12 +1,8 @@
-using System;
 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;
@@ -19,12 +15,6 @@ public class MersenneDivisorCycles
     // Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown
     private ulong[]? _smallCycles;
 
-    [ThreadStatic]
-    private static Dictionary<ulong, int>? s_factorCountsPool;
-
-    [ThreadStatic]
-    private static Dictionary<ulong, int>? s_factorScratchPool;
-
     public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles();
 
     public MersenneDivisorCycles()
@@ -61,15 +51,15 @@ public class MersenneDivisorCycles
         return snapshot;
     }
 
-    public static bool CycleEqualsExponent(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
+    public static bool CycleEqualsExponentCpu(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
     {
         Dictionary<ulong, FactorCacheEntry>? cache = null;
-        return CycleEqualsExponent(divisor, exponent, divisorData, ref cache);
+        return CycleEqualsExponentCpu(divisor, exponent, divisorData, ref cache);
     }
 
-    public static bool CycleEqualsExponent(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary<ulong, FactorCacheEntry>? factorCache)
+    public static bool CycleEqualsExponentCpu(ulong divisor, ulong exponent, in MontgomeryDivisorData divisorData, ref Dictionary<ulong, FactorCacheEntry>? factorCache)
     {
-        if (CycleEqualsExponentForMersenneCandidate(divisor, divisorData, exponent))
+        if (CycleEqualsExponentForMersenneCandidateCpu(divisor, divisorData, exponent))
         {
             return true;
         }
@@ -104,7 +94,7 @@ public class MersenneDivisorCycles
             factorCache = cache;
         }
 
-        if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL)
+        if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL)
         {
             return cycleLength == exponent;
         }
@@ -112,7 +102,7 @@ public class MersenneDivisorCycles
         return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
     }
 
-    public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
+    public static bool CycleEqualsExponentForMersenneCandidateCpu(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
     {
         if (!IsValidMersenneDivisorCandidate(divisor, exponent))
         {
@@ -271,7 +261,7 @@ public class MersenneDivisorCycles
         return order;
     }
 
-    public static bool TryCalculateCycleLengthForExponent(
+    public static bool TryCalculateCycleLengthForExponentCpu(
         ulong divisor,        
         ulong exponent,
         in MontgomeryDivisorData divisorData,
@@ -280,85 +270,62 @@ public class MersenneDivisorCycles
     {
         cycleLength = 0UL;
 
-        // The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
-        // so this power-of-two fast-path never triggers on that call path.
-        // The historical fast-path is left here for reference but commented out because it is unreachable
-        // on the EvenPerfectBitScanner pipeline.
-        // if ((divisor & (divisor - 1UL)) == 0UL)
-        // {
-        //     cycleLength = 1UL;
-        //     return true;
-        // }
+        if ((divisor & (divisor - 1UL)) == 0UL)
+        {
+            cycleLength = 1UL;
+            return true;
+        }
 
-        // EvenPerfectBitScanner filters exponents <= 3 before invoking this helper and each
-        // divisor candidate is at least 2 * exponent + 1, so this guard never matches.
-        // The fallback calculation is therefore unreachable on the active call path.
-        // if (divisor <= 3UL || exponent <= 1UL)
-        // {
-        //     cycleLength = CalculateCycleLength(divisor, divisorData);
-        //     return true;
-        // }
+        if (divisor <= 3UL || exponent <= 1UL)
+        {
+            cycleLength = CalculateCycleLength(divisor, divisorData);
+            return true;
+        }
 
-        ulong computedOrder = PrimeOrderCalculator.Calculate(
+        ulong order = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
             divisorData,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
             PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-        if (computedOrder != 0UL)
+        if (order != 0UL)
         {
-            cycleLength = computedOrder;
+            cycleLength = order;
             return true;
         }
 
         ulong phi = divisor - 1UL;
-        // Divisors reach this helper only when they exceed the cached table, which means
-        // divisor >= 2 * exponent + 1 and phi cannot underflow to zero here.
-        // The defensive underflow guard is kept for posterity but disabled for the EvenPerfectBitScanner path.
-        // if (phi == 0UL)
-        // {
-        //     return false;
-        // }
+        if (phi == 0UL)
+        {
+            return false;
+        }
 
         int twoCount = BitOperations.TrailingZeroCount(phi);
         ulong reducedPhi = phi >> twoCount;
-        // Every divisor candidate generated by the by-divisor scan satisfies q ≡ 1 (mod 2 * exponent),
-        // ensuring reducedPhi keeps a factor of exponent and remains non-zero in this context.
-        // The divisibility check is therefore unreachable for the current caller and is commented out.
-        // if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
-        // {
-        //     return false;
-        // }
+        if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
+        {
+            return false;
+        }
 
         ulong k = reducedPhi / exponent;
-        // reducedPhi always contains at least one factor of exponent, so dividing by exponent leaves k >= 1.
-        // if (k == 0UL)
-        // {
-        //     return false;
-        // }
+        if (k == 0UL)
+        {
+            return false;
+        }
 
-        Dictionary<ulong, int> factorCounts = AcquireFactorCountDictionary(8);
-        try
+        Dictionary<ulong, int> factorCounts = new(8);
+        if (twoCount > 0)
         {
-            // The divisor candidates are odd, so phi is even and this branch always records at least one factor of two.
-            // if (twoCount > 0)
-            // {
-            //     factorCounts[2UL] = twoCount;
-            // }
             factorCounts[2UL] = twoCount;
-
-            if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
-            {
-                return false;
-            }
-
-            cycleLength = ReduceOrder(divisorData, phi, factorCounts, exponent);
-            return true;
         }
-        finally
+
+        if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
         {
-            ReleaseFactorCountDictionary(factorCounts);
+            return false;
         }
+
+        cycleLength = ReduceOrderCpu(divisorData, phi, factorCounts);
+        return true;
     }
 
     private static bool AccumulateFactors(
@@ -372,119 +339,61 @@ public class MersenneDivisorCycles
             return true;
         }
 
-        FactorCacheEntry cachedFactorization;
-        PrimeFactorPartialResult? lease = null;
-        bool fromCache;
-
-        try
+        if (!TryGetFactorization(value, cache, cacheResult, out FactorCacheEntry factorization))
         {
-            if (!TryGetFactorization(value, cache, cacheResult, out cachedFactorization, out lease, out fromCache))
-            {
-                return false;
-            }
-
-            if (fromCache)
-            {
-                ulong[] primes = cachedFactorization.Primes;
-                byte[] exponents = cachedFactorization.Exponents;
-                int length = cachedFactorization.Count;
-
-                for (int i = 0; i < length; i++)
-                {
-                    AddFactor(counts, primes[i], exponents[i]);
-                }
-
-                return true;
-            }
-
-            if (lease is null || lease.Count == 0)
-            {
-                return true;
-            }
-
-            ReadOnlySpan<ulong> primesSpan = lease.Primes;
-            ReadOnlySpan<byte> exponentsSpan = lease.Exponents;
-            int spanLength = lease.Count;
-
-            for (int i = 0; i < spanLength; i++)
-            {
-                AddFactor(counts, primesSpan[i], exponentsSpan[i]);
-            }
-
-            return true;
+            return false;
         }
-        finally
+
+        ulong[] primes = factorization.Primes;
+        byte[] exponents = factorization.Exponents;
+        int length = factorization.Count;
+
+        for (int i = 0; i < length; i++)
         {
-            lease?.Dispose();
+            AddFactor(counts, primes[i], exponents[i]);
         }
+
+        return true;
     }
 
     private static bool TryGetFactorization(
         ulong value,
         Dictionary<ulong, FactorCacheEntry>? cache,
         bool cacheResult,
-        out FactorCacheEntry factorization,
-        out PrimeFactorPartialResult? lease,
-        out bool fromCache)
+        out FactorCacheEntry factorization)
     {
-        factorization = default;
-        lease = null;
-        fromCache = false;
-
-        if (value <= 1UL)
-        {
-            return true;
-        }
-
         if (cache is not null && cache.TryGetValue(value, out factorization))
         {
-            fromCache = true;
             return true;
         }
 
-        Dictionary<ulong, int> scratch = AcquireFactorScratchDictionary(8);
-        try
+        Dictionary<ulong, int> scratch = new(8);
+        if (!TryFactorIntoCountsInternal(value, scratch))
         {
-            if (!TryFactorIntoCountsInternal(value, scratch))
-            {
-                factorization = default;
-                return false;
-            }
-
-            int count = scratch.Count;
-            if (count == 0)
-            {
-                return true;
-            }
-
-            PrimeFactorPartialResult result = PrimeFactorPartialResult.Rent(count);
-            lease = result;
-
-            Span<ulong> primes = result.PrimeWriteSpan;
-            Span<byte> exponents = result.ExponentWriteSpan;
-            int index = 0;
-            foreach (KeyValuePair<ulong, int> entry in scratch)
-            {
-                primes[index] = entry.Key;
-                exponents[index] = checked((byte)entry.Value);
-                index++;
-            }
-
-            result.CommitCount(index);
-            result.Sort();
-
-            if (cache is not null && cacheResult)
-            {
-                FactorCacheEntry cacheEntry = result.ToCacheEntry();
-                cache[value] = cacheEntry;
-            }
-
-            return true;
+            factorization = default!;
+            return false;
         }
-        finally
+
+        int count = scratch.Count;
+        ulong[] primes = new ulong[count];
+        byte[] exponents = new byte[count];
+        int index = 0;
+        foreach (KeyValuePair<ulong, int> entry in scratch)
         {
-            ReleaseFactorScratchDictionary(scratch);
+            primes[index] = entry.Key;
+            exponents[index] = checked((byte)entry.Value);
+            index++;
         }
+
+        Array.Sort(primes, exponents);
+        factorization = new FactorCacheEntry(primes, exponents, count);
+
+        if (cache is not null && cacheResult)
+        {
+            cache[value] = factorization;
+        }
+
+        return true;
     }
 
     private static bool TryFactorIntoCountsInternal(ulong value, Dictionary<ulong, int> counts)
@@ -559,22 +468,15 @@ public class MersenneDivisorCycles
         }
     }
 
-    private static ulong ReduceOrder(
-        in MontgomeryDivisorData divisorData,
-        ulong initialOrder,
-        Dictionary<ulong, int> factorCounts,
-        ulong targetOrder)
+    private static ulong ReduceOrderCpu(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
     {
-        if (factorCounts.Count == 0 || initialOrder == targetOrder)
+        if (factorCounts.Count == 0)
         {
             return initialOrder;
         }
 
         int count = factorCounts.Count;
-        bool rentFromPool = count >= PerfectNumberConstants.PooledArrayThreshold;
-        ulong[] primes = rentFromPool
-            ? ThreadLocalArrayPool<ulong>.Shared.Rent(count)
-            : new ulong[count];
+        ulong[] primes = ArrayPool<ulong>.Shared.Rent(count);
         try
         {
             int index = 0;
@@ -603,13 +505,6 @@ public class MersenneDivisorCycles
                     if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL)
                     {
                         order = candidate;
-                        if (order == targetOrder)
-                        {
-                            // Early-exit once the order matches the tested exponent;
-                            // further factor processing would not change the outcome for the by-divisor scan.
-                            return order;
-                        }
-
                         continue;
                     }
 
@@ -621,10 +516,7 @@ public class MersenneDivisorCycles
         }
         finally
         {
-            if (rentFromPool)
-            {
-                ThreadLocalArrayPool<ulong>.Shared.Return(primes, clearArray: false);
-            }
+            ArrayPool<ulong>.Shared.Return(primes, clearArray: false);
         }
     }
 
@@ -635,20 +527,31 @@ public class MersenneDivisorCycles
             return 2UL;
         }
 
+        PollardRhoMontgomeryReducer reducer = PollardRhoMontgomeryReducer.Create(n);
+        ulong diff = 0UL;
+        ulong d = 0UL;
+
         while (true)
         {
-            ulong c = (NextRandomUInt64() % (n - 1UL)) + 1UL;
-            ulong x = (NextRandomUInt64() % (n - 2UL)) + 2UL;
-            ulong y = x;
-            ulong d = 1UL;
+            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)
             {
-                x = AdvancePolynomial(x, c, n);
-                y = AdvancePolynomial(y, c, n);
-                y = AdvancePolynomial(y, c, n);
+                xMontgomery = reducer.Advance(xMontgomery, cMontgomery);
+                xStandard = reducer.AdvanceStandard(xStandard, c);
 
-                ulong diff = x > y ? x - y : y - x;
+                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);
             }
 
@@ -659,13 +562,6 @@ 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)
@@ -697,60 +593,6 @@ public class MersenneDivisorCycles
         }
     }
 
-    // Thread-local dictionary pools keep the factoring helpers reusable so the GC stays out of the by-divisor CPU loop.
-    private static Dictionary<ulong, int> AcquireFactorCountDictionary(int capacityHint)
-    {
-        Dictionary<ulong, int>? dictionary = s_factorCountsPool;
-        if (dictionary is null)
-        {
-            return new Dictionary<ulong, int>(capacityHint);
-        }
-
-        s_factorCountsPool = null;
-        dictionary.Clear();
-        dictionary.EnsureCapacity(capacityHint);
-        return dictionary;
-    }
-
-    private static void ReleaseFactorCountDictionary(Dictionary<ulong, int> dictionary)
-    {
-        dictionary.Clear();
-        if (s_factorCountsPool is null)
-        {
-            s_factorCountsPool = dictionary;
-        }
-    }
-
-    private static Dictionary<ulong, int> AcquireFactorScratchDictionary(int capacityHint)
-    {
-        Dictionary<ulong, int>? dictionary = s_factorScratchPool;
-        if (dictionary is null)
-        {
-            return new Dictionary<ulong, int>(capacityHint);
-        }
-
-        s_factorScratchPool = null;
-        dictionary.Clear();
-        dictionary.EnsureCapacity(capacityHint);
-        return dictionary;
-    }
-
-    private static void ReleaseFactorScratchDictionary(Dictionary<ulong, int> dictionary)
-    {
-        dictionary.Clear();
-        if (s_factorScratchPool is null)
-        {
-            s_factorScratchPool = dictionary;
-        }
-    }
-
-    private static ulong NextRandomUInt64()
-    {
-        Span<byte> buffer = stackalloc byte[8];
-        RandomNumberGenerator.Fill(buffer);
-        return BinaryPrimitives.ReadUInt64LittleEndian(buffer);
-    }
-
     public readonly struct FactorCacheEntry
     {
         public FactorCacheEntry(ulong[] primes, byte[] exponents, int count)
@@ -810,11 +652,7 @@ public class MersenneDivisorCycles
             tasks[taskIndex] = Task.Run(() =>
             {
                 ulong rangeLength = localEnd - localStart + 1UL;
-                int localLength = checked((int)rangeLength);
-                bool rentLocalCycles = localLength >= PerfectNumberConstants.PooledArrayThreshold;
-                (ulong divisor, ulong cycle)[] localCycles = rentLocalCycles
-                    ? ThreadLocalArrayPool<(ulong, ulong)>.Shared.Rent(localLength)
-                    : new (ulong divisor, ulong cycle)[localLength];
+                (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength));
                 int localCycleIndex = 0;
 
                 try
@@ -862,10 +700,7 @@ public class MersenneDivisorCycles
                 }
                 finally
                 {
-                    if (rentLocalCycles)
-                    {
-                        ThreadLocalArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false);
-                    }
+                    ArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false);
                 }
             });
         }
@@ -880,7 +715,7 @@ public class MersenneDivisorCycles
         ulong start = 3UL;
         var locker = new object();
 
-        var pool = ThreadLocalArrayPool<ulong>.Shared;
+        var pool = ArrayPool<ulong>.Shared;
         ulong batchSizeUL = (ulong)batchSize, d, end;
         int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx;
         ulong[] divisors, outCycles, validDivisors;
@@ -1085,34 +920,29 @@ public class MersenneDivisorCycles
 
     internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength)
     {
-        // The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
-        // so this power-of-two fast-path never triggers on that call path.
-        // The early exit is retained for reference but commented out because it is unreachable on the scanner path.
-        // if ((divisor & (divisor - 1UL)) == 0UL)
-        // {
-        //     cycleLength = 1UL;
-        //     return true;
-        // }
+        if ((divisor & (divisor - 1UL)) == 0UL)
+        {
+            cycleLength = 1UL;
+            return true;
+        }
 
-        // By the divisor-by-divisor pipeline, candidates have the shape 2 * exponent * k + 1 with exponent >= 5,
-        // so the divisor never drops below 11 on this call path.
-        // if (divisor <= 3UL)
-        // {
-        //     cycleLength = CalculateCycleLengthFallback(divisor);
-        //     return true;
-        // }
+        if (divisor <= 3UL)
+        {
+            cycleLength = CalculateCycleLengthFallback(divisor);
+            return true;
+        }
 
         if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
         {
-            ulong computedOrder = PrimeOrderCalculator.Calculate(
+            ulong order = PrimeOrderCalculator.Calculate(
                     divisor,
                     previousOrder: null,
                     divisorData,
                     PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
                     PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-            if (computedOrder != 0UL)
+            if (order != 0UL)
             {
-                cycleLength = computedOrder;
+                cycleLength = order;
                 return true;
             }
         }

PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs

@@ -18,8 +18,7 @@ public static class MersenneNumberDivisorByDivisorTester
             Action markComposite,
             Action clearComposite,
             Action<ulong, bool, bool, bool> printResult,
-            int threadCount,
-            TimeSpan? primeTestLimit)
+            int threadCount)
     {
         if (candidates.Count == 0)
         {
@@ -27,8 +26,6 @@ public static class MersenneNumberDivisorByDivisorTester
             return;
         }
 
-        _ = threadCount;
-
         candidates.Sort();
         Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
         int candidateCount = candidateSpan.Length;
@@ -39,34 +36,65 @@ public static class MersenneNumberDivisorByDivisorTester
         if (previousResults is not null && previousResults.Count > 0)
         {
             int recordedCount = previousResults.Count;
-            ulong[] recordedCandidates = ArrayPool<ulong>.Shared.Rent(recordedCount);
-            try
+            bool preferSortedMerge = recordedCount <= 4096 || (recordedCount * 4L) <= candidateCount;
+            if (preferSortedMerge)
             {
-                Span<ulong> recordedSpan = recordedCandidates.AsSpan(0, recordedCount);
-                int recordedIndex = 0;
-                foreach (ulong value in previousResults.Keys)
+                ulong[] recordedCandidates = ArrayPool<ulong>.Shared.Rent(recordedCount);
+                try
                 {
-                    recordedSpan[recordedIndex++] = value;
+                    Span<ulong> recordedSpan = recordedCandidates.AsSpan(0, recordedCount);
+                    int recordedIndex = 0;
+                    foreach (ulong value in previousResults.Keys)
+                    {
+                        recordedSpan[recordedIndex++] = value;
+                    }
+
+                    recordedSpan.Sort();
+
+                    int writeIndex = 0;
+                    int skipIndex = 0;
+
+                    for (int readIndex = 0; readIndex < candidateCount; readIndex++)
+                    {
+                        ulong candidate = candidateSpan[readIndex];
+
+                        while (skipIndex < recordedCount && recordedSpan[skipIndex] < candidate)
+                        {
+                            skipIndex++;
+                        }
+
+                        if (skipIndex < recordedCount && recordedSpan[skipIndex] == candidate)
+                        {
+                            skippedByPreviousResults++;
+                            skipIndex++;
+                            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;
+                    }
                 }
-
-                recordedSpan.Sort();
-
+                finally
+                {
+                    ArrayPool<ulong>.Shared.Return(recordedCandidates, clearArray: true);
+                }
+            }
+            else
+            {
                 int writeIndex = 0;
-                int skipIndex = 0;
-
                 for (int readIndex = 0; readIndex < candidateCount; readIndex++)
                 {
                     ulong candidate = candidateSpan[readIndex];
 
-                    while (skipIndex < recordedCount && recordedSpan[skipIndex] < candidate)
-                    {
-                        skipIndex++;
-                    }
-
-                    if (skipIndex < recordedCount && recordedSpan[skipIndex] == candidate)
+                    if (previousResults.ContainsKey(candidate))
                     {
                         skippedByPreviousResults++;
-                        skipIndex++;
                         continue;
                     }
 
@@ -76,14 +104,10 @@ public static class MersenneNumberDivisorByDivisorTester
                 if (writeIndex < candidateCount)
                 {
                     candidates.RemoveRange(writeIndex, candidateCount - writeIndex);
-                    candidateSpan = CollectionsMarshal.AsSpan(candidates);
+                    candidateSpan = CollectionsMarshal.AsSpan(candidates); // candidateSpan now tracks the reduced list.
                     candidateCount = candidateSpan.Length;
                 }
             }
-            finally
-            {
-                ArrayPool<ulong>.Shared.Return(recordedCandidates, clearArray: true);
-            }
         }
 
         if (skippedByPreviousResults > 0)
@@ -113,12 +137,9 @@ public static class MersenneNumberDivisorByDivisorTester
         ulong maxPrime = 0UL;
         int primeWriteIndex = 0;
 
-		// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
-		// We don't need to additionally check it.
-
-        // using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
-        // bool hasPrime = primeEnumerator.MoveNext();
-        // ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
+        using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
+        bool hasPrime = primeEnumerator.MoveNext();
+        ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
 
         for (int index = startIndex; index < candidateCount; index++)
         {
@@ -131,29 +152,25 @@ public static class MersenneNumberDivisorByDivisorTester
                 continue;
             }
 
-			// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
-			// We don't need to additionally check it.
-			
-			// while (hasPrime && currentPrime < candidate)
-			// {
-			//     hasPrime = primeEnumerator.MoveNext();
-			//     if (!hasPrime)
-			//     {
-			//         break;
-			//     }
+            while (hasPrime && currentPrime < candidate)
+            {
+                hasPrime = primeEnumerator.MoveNext();
+                if (!hasPrime)
+                {
+                    break;
+                }
 
-			//     currentPrime = primeEnumerator.Current;
-			// }
+                currentPrime = primeEnumerator.Current;
+            }
 
-			// if (!hasPrime || currentPrime != candidate)
-			// if (!Number.IsPrime(candidate))
-			// {
-			//     markComposite();
-			//     printResult(candidate, false, false, false);
-			//     continue;
-			// }
+            if (!hasPrime || currentPrime != candidate)
+            {
+                markComposite();
+                printResult(candidate, false, false, false);
+                continue;
+            }
 
-			primesSpan[primeWriteIndex++] = candidate;
+            primesSpan[primeWriteIndex++] = candidate;
 
             if (candidate > maxPrime)
             {
@@ -244,7 +261,7 @@ public static class MersenneNumberDivisorByDivisorTester
 
         void ProcessPrime(ulong prime)
         {
-            bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted, primeTestLimit);
+            bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted);
 
             if (!isPrime)
             {
@@ -257,52 +274,46 @@ public static class MersenneNumberDivisorByDivisorTester
             printResult(prime, true, divisorsExhausted, true);
         }
 
-        if (workerCount > filteredPrimes.Count)
-        {
-            workerCount = filteredPrimes.Count;
-        }
-
-        if (workerCount <= 1)
+        if (workerCount == 1)
         {
             foreach (ulong prime in filteredPrimes)
             {
                 ProcessPrime(prime);
             }
-
-            return;
         }
-
-        (int Start, int End)[] workerRanges = new (int, int)[workerCount];
-        int totalPrimeCount = filteredPrimes.Count;
-        int baseGroupSize = totalPrimeCount / workerCount;
-        int remaining = totalPrimeCount % workerCount;
-        int groupStartIndex = 0;
-
-        for (int workerIndex = 0; workerIndex < workerCount; workerIndex++)
+        else
         {
-            int groupSize = baseGroupSize;
-            if (workerIndex < remaining)
+            ParallelOptions options = new()
             {
-                groupSize++;
-            }
+                MaxDegreeOfParallelism = workerCount
+            };
 
-            startIndex = groupStartIndex;
-            int endIndex = startIndex + groupSize;
-            workerRanges[workerIndex] = (startIndex, endIndex);
-            groupStartIndex = endIndex;
+            primeCount = filteredPrimes.Count; // Reusing primeCount to track the filtered count.
+            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]);
+                        }
+                    });
         }
-
-        ParallelOptions parallelOptions = new()
-        {
-            MaxDegreeOfParallelism = workerCount
-        };
-
-        Parallel.ForEach(workerRanges, parallelOptions, range =>
-        {
-            for (int candidateIndex = range.Start; candidateIndex < range.End; candidateIndex++)
-            {
-                ProcessPrime(filteredPrimes[candidateIndex]);
-            }
-        });
     }
 }

PerfectNumbers.Core/MontgomeryDivisorData.cs

@@ -45,4 +45,32 @@ public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong
             montgomeryTwo,
             montgomeryTwoSquared);
     }
+
+}
+
+internal static class MontgomeryDivisorDataCache
+{
+    [ThreadStatic]
+    private static bool s_hasCachedValue;
+
+    [ThreadStatic]
+    private static ulong s_cachedModulus;
+
+    [ThreadStatic]
+    private static MontgomeryDivisorData s_cachedData;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static MontgomeryDivisorData Get(ulong modulus)
+    {
+        if (s_hasCachedValue && s_cachedModulus == modulus)
+        {
+            return s_cachedData;
+        }
+
+        MontgomeryDivisorData computed = MontgomeryDivisorData.FromModulus(modulus);
+        s_cachedModulus = modulus;
+        s_cachedData = computed;
+        s_hasCachedValue = true;
+        return computed;
+    }
 }

PerfectNumbers.Core/MutableUInt128.cs

@@ -1,243 +0,0 @@
-using System;
-using System.Numerics;
-using System.Runtime.CompilerServices;
-
-namespace PerfectNumbers.Core;
-
-/// <summary>
-/// Mutable helper for 128-bit arithmetic derived from the .NET runtime UInt128 implementation.
-/// Portions adapted from https://github.com/dotnet/runtime/blob/v8.0.0/src/libraries/System.Private.CoreLib/src/System/UInt128.cs.
-/// </summary>
-public struct MutableUInt128
-{
-    public ulong High;
-    public ulong Low;
-
-    public MutableUInt128(ulong value)
-    {
-        Low = value;
-    }
-
-    public MutableUInt128(ulong high, ulong low)
-    {
-        High = high;
-        Low = low;
-    }
-
-    public MutableUInt128(UInt128 value)
-    {
-        High = (ulong)(value >> 64);
-        Low = (ulong)value;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Set(ulong value)
-    {
-        High = 0UL;
-        Low = value;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Set(UInt128 value)
-    {
-        High = (ulong)(value >> 64);
-        Low = (ulong)value;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Set(ulong high, ulong low)
-    {
-        High = high;
-        Low = low;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public readonly UInt128 ToUInt128() => ((UInt128)High << 64) | Low;
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Add(ulong value)
-    {
-        ulong low = Low + value;
-        if (low < value)
-        {
-            High++;
-        }
-
-        Low = low;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Add(MutableUInt128 value)
-    {
-        ulong temp = Low + value.Low;
-        Low = temp;
-        temp = temp < Low ? 1UL : 0UL;
-		temp += High;
-		temp += value.High;
-        High = temp;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Sub(ulong value)
-    {
-        ulong temp = Low;
-        ulong result = temp - value;
-        temp = result > temp ? 1UL : 0UL;
-
-        Low = result;
-        High -= temp;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Sub(MutableUInt128 value)
-    {
-        ulong temp = Low;
-        ulong result = temp - value.Low;
-        temp = temp < value.Low ? 1UL : 0UL;
-
-        Low = result;
-        High = High - value.High - temp;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void Mul(ulong value)
-    {
-        ulong currentHigh = High;
-        ulong currentLow = Low;
-
-        ulong temp = Mul64(currentLow, value, out currentLow);
-        Mul64(currentHigh, value, out currentHigh);
-        temp += currentHigh;
-
-        Low = currentLow;
-        High = temp;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Mul64(ulong left, ulong right, out ulong low)
-    {
-        const ulong Mask32 = 0xFFFF_FFFFUL;
-
-        ulong leftHigh = left >> 32;
-        ulong rightLow = right & Mask32;
-
-        ulong leftLow = left & Mask32;
-        ulong lowLow = leftLow * rightLow;
-        ulong temp = right >> 32;
-        leftLow *= temp;
-        rightLow *= leftHigh;
-        leftHigh *= temp;
-
-        temp = (lowLow >> 32) + (leftLow & Mask32) + (rightLow & Mask32);
-        leftHigh += (leftLow >> 32) + (rightLow >> 32) + (temp >> 32);
-
-        low = ((temp & Mask32) << 32) | (lowLow & Mask32);
-        return leftHigh;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public void MulAdd(ulong multiplier, ulong addend)
-    {
-        const ulong Mask32 = 0xFFFF_FFFFUL;
-
-        ulong multiplierLow = multiplier & Mask32;
-        ulong multiplierHigh = multiplier >> 32;
-
-        ulong product11 = Low;
-        ulong product01 = product11 & Mask32;
-        product11 >>= 32;
-
-        ulong product00 = product01 * multiplierLow;
-        product01 *= multiplierHigh;
-        ulong product10 = product11 * multiplierLow;
-        product11 *= multiplierHigh;
-
-        ulong high = (product00 >> 32) + (product01 & Mask32) + (product10 & Mask32);
-        product11 = product11 + (product01 >> 32) + (product10 >> 32) + (high >> 32);
-        product10 = ((high & Mask32) << 32) | (product00 & Mask32);
-
-        high = High;
-        product00 = high & Mask32;
-        high >>= 32;
-
-        product01 = product00 * multiplierLow;
-        product00 *= multiplierHigh;
-        high *= multiplierLow;
-
-        high = (product01 >> 32) + (product00 & Mask32) + (high & Mask32);
-		product01 = ((high & Mask32) << 32) | (product01 & Mask32);
-
-		// Calculate low. We're assigning addend to local variable for a tiny bit better performance.
-		product00 = addend;
-		high = product10 + product00;
-        Low = high;
-
-		// Calculate high including carry
-		high = high < product00 ? 1UL : 0UL;
-		high += product11 + product01;
-		
-        // Discard overflow above 128 bits to match Multiply's truncation semantics.
-        High = high;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public readonly ulong Mod(ulong modulus)
-	{
-		// This case will never happen in production code.
-		// if (modulus == 0UL)
-		// {
-		//     throw new DivideByZeroException();
-		// }
-
-		ulong high = High;
-		if (high == 0UL)
-        {
-            return Low % modulus;
-        }
-
-        MutableUInt128 remainder = default;
-        int highBitLength = 64 - BitOperations.LeadingZeroCount(high);
-        if (highBitLength > 0)
-        {
-            ProcessWord(ref remainder, high, highBitLength, modulus);
-        }
-
-        ProcessWord(ref remainder, Low, 64, modulus);
-        return remainder.Low;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static void ProcessWord(ref MutableUInt128 remainder, ulong word, int bitCount, ulong modulus)
-    {
-        for (int bit = bitCount - 1; bit >= 0; bit--)
-        {
-            remainder.ShiftLeft();
-            if (((word >> bit) & 1UL) != 0)
-            {
-                remainder.Add(1UL);
-            }
-
-            remainder.Reduce(modulus);
-        }
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private void ShiftLeft()
-    {
-        High = (High << 1) | (Low >> 63);
-        Low <<= 1;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private void Reduce(ulong modulus)
-    {
-        if (High != 0UL || Low >= modulus)
-        {
-            Sub(modulus);
-            if (High != 0UL || Low >= modulus)
-            {
-                Sub(modulus);
-            }
-        }
-    }
-}

PerfectNumbers.Core/PartialFactorResult.cs

@@ -1,4 +1,4 @@
-using System;
+using System.Buffers;
 
 namespace PerfectNumbers.Core
 {
@@ -7,88 +7,70 @@ namespace PerfectNumbers.Core
         [ThreadStatic]
         private static PartialFactorResult? s_poolHead;
 
-        public static readonly PartialFactorResult Empty = new()
-        {
-            Cofactor = 1UL,
-            FullyFactored = true,
-        };
-
         public FactorEntry[]? Factors;
+        private PartialFactorResult? _next;
         public ulong Cofactor;
         public bool FullyFactored;
         public int Count;
 
-        private PartialFactorResult? _next;
-        private bool _factorsFromPool;
-
         private PartialFactorResult()
         {
         }
 
-        public static PartialFactorResult Rent(in FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count)
-        {
-            PartialFactorResult? instance = s_poolHead;
-            if (instance is null)
-            {
-                instance = new PartialFactorResult();
-            }
-            else
-            {
-                s_poolHead = instance._next;
-            }
+		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;
+			}
 
-            FactorEntry[]? factorArray = factors;
-            instance._next = null;
-            instance.Factors = factorArray;
-            instance.Cofactor = cofactor;
-            instance.FullyFactored = fullyFactored;
-            instance.Count = count;
-            instance._factorsFromPool = factorArray is not null && factorArray.Length >= PerfectNumberConstants.PooledArrayThreshold;
-            return instance;
-        }
+			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[] single = new FactorEntry[1];
-                single[0] = new FactorEntry(prime);
-                return Rent(single, 1UL, true, 1);
+                FactorEntry[] local = ArrayPool<FactorEntry>.Shared.Rent(1);
+                local[0] = new FactorEntry(prime, 1);
+                return Rent(local, 1UL, true, 1);
             }
 
-            int newCount = Count + 1;
-            FactorEntry[] extended = newCount >= PerfectNumberConstants.PooledArrayThreshold
-                ? ThreadLocalArrayPool<FactorEntry>.Shared.Rent(newCount)
-                : new FactorEntry[newCount];
-
+            FactorEntry[] extended = ArrayPool<FactorEntry>.Shared.Rent(Count + 1);
             Array.Copy(source, 0, extended, 0, Count);
-            extended[Count] = new FactorEntry(prime);
-            Span<FactorEntry> span = extended.AsSpan(0, newCount);
+            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, newCount);
+            return Rent(extended, 1UL, true, Count + 1);
         }
 
         public void Dispose()
-        {
-            if (this == Empty)
-            {
-                return;
-            }
+		{
+			if (this != Empty)
+			{
+				FactorEntry[]? factors = Factors;
+				if (factors is not null)
+				{
+					Factors = null;
+					ArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
+				}
 
-            FactorEntry[]? factors = Factors;
-            if (factors is not null)
-            {
-                Factors = null;
-                if (_factorsFromPool)
-                {
-                    ThreadLocalArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
-                }
-            }
-
-            _factorsFromPool = false;
-            _next = s_poolHead;
-            s_poolHead = this;
+				_next = s_poolHead;
+				s_poolHead = this;
+			}
         }
     }
-}
+}

PerfectNumbers.Core/PerfectNumbers.Core.csproj

@@ -22,7 +22,4 @@
     <InternalsVisibleTo Include="EvenPerfectBitScanner.Benchmarks" />
   </ItemGroup>
 
-  <ItemGroup>
-    <Compile Remove=".temp.cs" />
-  </ItemGroup>
 </Project>

PerfectNumbers.Core/PollardRhoMontgomeryReducer.cs

@@ -0,0 +1,119 @@
+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/PrimeFactorPartialResult.cs

@@ -1,150 +0,0 @@
-using System;
-using System.Buffers;
-
-namespace PerfectNumbers.Core
-{
-    internal sealed class PrimeFactorPartialResult : IDisposable
-    {
-        [ThreadStatic]
-        private static PrimeFactorPartialResult? s_poolHead;
-
-        private PrimeFactorPartialResult? _next;
-        private ulong[]? _primes;
-        private byte[]? _exponents;
-        private bool _arraysFromPool;
-        private int _capacity;
-
-        private PrimeFactorPartialResult()
-        {
-        }
-
-        public static PrimeFactorPartialResult Rent(int capacity)
-        {
-            PrimeFactorPartialResult? instance = s_poolHead;
-            if (instance is null)
-            {
-                instance = new PrimeFactorPartialResult();
-            }
-            else
-            {
-                s_poolHead = instance._next;
-            }
-
-            instance._next = null;
-            instance.EnsureCapacity(capacity);
-            instance.Count = 0;
-            return instance;
-        }
-
-        public int Count { get; private set; }
-
-        internal ReadOnlySpan<ulong> Primes => _primes is null ? ReadOnlySpan<ulong>.Empty : _primes.AsSpan(0, Count);
-
-        internal ReadOnlySpan<byte> Exponents => _exponents is null ? ReadOnlySpan<byte>.Empty : _exponents.AsSpan(0, Count);
-
-        internal Span<ulong> PrimeWriteSpan => _primes is null ? Span<ulong>.Empty : _primes.AsSpan(0, _capacity);
-
-        internal Span<byte> ExponentWriteSpan => _exponents is null ? Span<byte>.Empty : _exponents.AsSpan(0, _capacity);
-
-        internal void CommitCount(int count)
-        {
-            Count = count;
-        }
-
-        internal void Sort()
-        {
-            if (Count == 0 || _primes is null || _exponents is null)
-            {
-                return;
-            }
-
-            Array.Sort(_primes, _exponents, 0, Count);
-        }
-
-        internal MersenneDivisorCycles.FactorCacheEntry ToCacheEntry()
-        {
-            if (Count == 0 || _primes is null || _exponents is null)
-            {
-                return default;
-            }
-
-            ulong[] primes = new ulong[Count];
-            byte[] exponents = new byte[Count];
-            Array.Copy(_primes, 0, primes, 0, Count);
-            Array.Copy(_exponents, 0, exponents, 0, Count);
-            return new MersenneDivisorCycles.FactorCacheEntry(primes, exponents, Count);
-        }
-
-        public void Dispose()
-        {
-            if (_arraysFromPool)
-            {
-                if (_primes is not null)
-                {
-                    ThreadLocalArrayPool<ulong>.Shared.Return(_primes, clearArray: false);
-                }
-
-                if (_exponents is not null)
-                {
-                    ThreadLocalArrayPool<byte>.Shared.Return(_exponents, clearArray: false);
-                }
-
-                _arraysFromPool = false;
-                _primes = null;
-                _exponents = null;
-            }
-
-            _capacity = 0;
-            Count = 0;
-
-            _next = s_poolHead;
-            s_poolHead = this;
-        }
-
-        private void EnsureCapacity(int capacity)
-        {
-            if (capacity <= 0)
-            {
-                _capacity = 0;
-                return;
-            }
-
-            if (_primes is not null && _primes.Length >= capacity)
-            {
-                _capacity = capacity;
-                return;
-            }
-
-            if (_arraysFromPool)
-            {
-                if (_primes is not null)
-                {
-                    ThreadLocalArrayPool<ulong>.Shared.Return(_primes, clearArray: false);
-                }
-
-                if (_exponents is not null)
-                {
-                    ThreadLocalArrayPool<byte>.Shared.Return(_exponents, clearArray: false);
-                }
-
-                _arraysFromPool = false;
-                _primes = null;
-                _exponents = null;
-            }
-
-            if (capacity >= PerfectNumberConstants.PooledArrayThreshold)
-            {
-                _primes = ThreadLocalArrayPool<ulong>.Shared.Rent(capacity);
-                _exponents = ThreadLocalArrayPool<byte>.Shared.Rent(capacity);
-                _arraysFromPool = true;
-            }
-            else
-            {
-                _primes = new ulong[capacity];
-                _exponents = new byte[capacity];
-            }
-
-            _capacity = capacity;
-        }
-    }
-}

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -9,6 +9,9 @@ namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
+    public static UInt128 Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
+        => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
+
     public static UInt128 Calculate(
         UInt128 prime,
         UInt128? previousOrder,
@@ -162,7 +165,7 @@ internal static partial class PrimeOrderCalculator
                 return order;
             }
 
-            tempArray = ThreadLocalArrayPool<FactorEntry128>.Shared.Rent(length + 1);
+            tempArray = ArrayPool<FactorEntry128>.Shared.Rent(length + 1);
             Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
             factorSpan.CopyTo(buffer);
 
@@ -200,7 +203,7 @@ internal static partial class PrimeOrderCalculator
         {
             if (tempArray is not null)
             {
-                ThreadLocalArrayPool<FactorEntry128>.Shared.Return(tempArray, clearArray: false);
+                ArrayPool<FactorEntry128>.Shared.Return(tempArray, clearArray: false);
             }
         }
     }
@@ -664,7 +667,7 @@ internal static partial class PrimeOrderCalculator
             return new PartialFactorResult128(null, value, false, 0);
         }
 
-        FactorEntry128[]? rented = counts.Count > 0 ? ThreadLocalArrayPool<FactorEntry128>.Shared.Rent(counts.Count) : null;
+        FactorEntry128[]? rented = counts.Count > 0 ? ArrayPool<FactorEntry128>.Shared.Rent(counts.Count) : null;
         int index = 0;
         if (rented is not null)
         {
@@ -682,7 +685,7 @@ internal static partial class PrimeOrderCalculator
             span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
             resultArray = new FactorEntry128[index];
             span.CopyTo(resultArray);
-            ThreadLocalArrayPool<FactorEntry128>.Shared.Return(rented, clearArray: false);
+            ArrayPool<FactorEntry128>.Shared.Return(rented, clearArray: false);
         }
 
         bool fullyFactored = cofactor == UInt128.One;
@@ -1062,20 +1065,15 @@ internal static partial class PrimeOrderCalculator
 
     private readonly struct FactorEntry128
     {
-		public FactorEntry128(UInt128 value)
-		{
-			Value = value;
-		}
+        public FactorEntry128(UInt128 value, int exponent)
+        {
+            Value = value;
+            Exponent = exponent;
+        }
 
-		public FactorEntry128(UInt128 value, int exponent)
-		{
-			Value = value;
-			Exponent = exponent;
-		}
+        public UInt128 Value { get; }
 
-        public readonly UInt128 Value;
-
-		public readonly int Exponent = 1;
+        public int Exponent { get; }
     }
 
     private readonly struct PartialFactorResult128
@@ -1088,13 +1086,13 @@ internal static partial class PrimeOrderCalculator
             Count = count;
         }
 
-		public readonly FactorEntry128[]? Factors;
+        public FactorEntry128[]? Factors { get; }
 
-		public readonly UInt128 Cofactor;
+        public UInt128 Cofactor { get; }
 
-		public readonly bool FullyFactored;
+        public bool FullyFactored { get; }
 
-		public readonly int Count;
+        public int Count { get; }
 
         public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0);
 
@@ -1102,13 +1100,14 @@ internal static partial class PrimeOrderCalculator
         {
             if (Factors is null)
             {
-                FactorEntry128[] local = [new FactorEntry128(prime)];
-				return new PartialFactorResult128(local, UInt128.One, true, 1);
+                FactorEntry128[] local = new FactorEntry128[1];
+                local[0] = new FactorEntry128(prime, 1);
+                return new PartialFactorResult128(local, UInt128.One, true, 1);
             }
 
             FactorEntry128[] extended = new FactorEntry128[Count + 1];
             Array.Copy(Factors, extended, Count);
-            extended[Count] = new FactorEntry128(prime);
+            extended[Count] = new FactorEntry128(prime, 1);
             return new PartialFactorResult128(extended, UInt128.One, true, Count + 1);
         }
     }

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -17,14 +17,14 @@ internal static partial class PrimeOrderCalculator
     }
 
     internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode)
-    {
-        public static PrimeOrderSearchConfig HeuristicDefault => new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 256, maxPowChecks: 24, PrimeOrderMode.Heuristic);
+	{
+		public static PrimeOrderSearchConfig HeuristicDefault => new(smallFactorLimit: 100_000, pollardRhoMilliseconds: 128, maxPowChecks: 24, PrimeOrderMode.Heuristic);
         public static PrimeOrderSearchConfig StrictDefault => new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict);
 
-        public readonly uint SmallFactorLimit = smallFactorLimit;
-        public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds;
-        public readonly int MaxPowChecks = maxPowChecks;
-        public readonly PrimeOrderMode Mode = mode;
+		public readonly uint SmallFactorLimit = smallFactorLimit;
+		public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds;
+		public readonly int MaxPowChecks = maxPowChecks;
+		public readonly PrimeOrderMode Mode = mode;		
     }
 
     internal enum PrimeOrderHeuristicDevice
@@ -39,12 +39,8 @@ internal static partial class PrimeOrderCalculator
     [ThreadStatic]
     private static bool s_allowGpuPow2;
 
-    [ThreadStatic]
-    private static bool s_allowGpuCycle;
-
     private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
 
-    private static bool IsGpuCycleAllowed => s_pow2ModeInitialized && s_allowGpuCycle;
 
     [ThreadStatic]
     private static PrimeOrderHeuristicDevice s_deviceMode;
@@ -61,53 +57,31 @@ internal static partial class PrimeOrderCalculator
     [ThreadStatic]
     private static Dictionary<ulong, int>? s_partialFactorCounts;
 
-    [ThreadStatic]
-    private static ulong s_cachedDivisorCycleModulus;
-
-    [ThreadStatic]
-    private static ulong s_cachedDivisorCycleLength;
-
-    [ThreadStatic]
-    private static bool s_hasCachedDivisorCycle;
-
-    [ThreadStatic]
-    private static int s_divisorCycleSuppressionDepth;
-
-    [ThreadStatic]
-    private static MutableUInt128 s_pollardRhoPolynomialBuffer;
-
     private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu;
 
-    private readonly struct Pow2ModeScope
-    {
-        private readonly bool _previousInitialized;
-        private readonly bool _previousAllowPow2;
-        private readonly bool _previousAllowCycle;
-        private readonly PrimeOrderHeuristicDevice _previousDevice;
+	private readonly struct Pow2ModeScope
+	{
+		private readonly bool _previousInitialized;
+		private readonly bool _previousAllow;
+		private readonly PrimeOrderHeuristicDevice _previousDevice;
 
-        public Pow2ModeScope(PrimeOrderHeuristicDevice device)
-        {
-            _previousInitialized = s_pow2ModeInitialized;
-            _previousAllowPow2 = s_allowGpuPow2;
-            _previousAllowCycle = s_allowGpuCycle;
-            _previousDevice = s_deviceMode;
-            s_pow2ModeInitialized = true;
-            bool gpuAvailable = !GpuContextPool.ForceCpu;
-            bool allowCpuPow2 = gpuAvailable && device == PrimeOrderHeuristicDevice.Cpu;
-            bool allowCpuCycle = gpuAvailable && device == PrimeOrderHeuristicDevice.Cpu && GpuKernelPool.CpuCycleUsesGpu;
-            s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu || allowCpuPow2;
-            s_allowGpuCycle = device == PrimeOrderHeuristicDevice.Gpu || allowCpuCycle;
-            s_deviceMode = device;
-        }
+		public Pow2ModeScope(PrimeOrderHeuristicDevice device)
+		{
+			_previousInitialized = s_pow2ModeInitialized;
+			_previousAllow = s_allowGpuPow2;
+			_previousDevice = s_deviceMode;
+			s_pow2ModeInitialized = true;
+			s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu;
+			s_deviceMode = device;
+		}
 
-        public void Dispose()
-        {
-            s_allowGpuPow2 = _previousAllowPow2;
-            s_allowGpuCycle = _previousAllowCycle;
-            s_pow2ModeInitialized = _previousInitialized;
-            s_deviceMode = _previousDevice;
-        }
-    }
+		public void Dispose()
+		{
+			s_allowGpuPow2 = _previousAllow;
+			s_pow2ModeInitialized = _previousInitialized;
+			s_deviceMode = _previousDevice;
+		}
+	}
 
     private readonly struct DebugLoggingScope
     {
@@ -125,45 +99,10 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool IsDivisorCycleSuppressed => s_divisorCycleSuppressionDepth > 0;
-
-    private readonly struct DivisorCycleSuppressionScope
-    {
-        private readonly bool _suppress;
-
-        private DivisorCycleSuppressionScope(bool suppress)
-        {
-            _suppress = suppress;
-        }
-
-        public static readonly DivisorCycleSuppressionScope Suppressed = new(true);
-
-        public static readonly DivisorCycleSuppressionScope Unsuppressed = new(false);
-
-        public void Enter()
-        {
-            if (_suppress)
-            {
-                s_divisorCycleSuppressionDepth++;
-            }
-        }
-
-        public void Exit()
-        {
-            if (_suppress)
-            {
-                s_divisorCycleSuppressionDepth--;
-            }
-        }
-    }
-
     private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device);
 
     private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled);
 
-    private static DivisorCycleSuppressionScope SuppressDivisorCycleUsage(bool suppress) =>
-        suppress ? DivisorCycleSuppressionScope.Suppressed : DivisorCycleSuppressionScope.Unsuppressed;
-
     [Conditional("DEBUG")]
     private static void DebugLog(string message)
     {
@@ -182,6 +121,13 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
+    public static ulong Calculate(
+        ulong prime,
+        ulong? previousOrder,
+        in MontgomeryDivisorData divisorData,
+        in PrimeOrderSearchConfig config)
+        => Calculate(prime, previousOrder, divisorData, config, PrimeOrderHeuristicDevice.Cpu);
+
     public static ulong Calculate(
         ulong prime,
         ulong? previousOrder,
@@ -197,16 +143,13 @@ internal static partial class PrimeOrderCalculator
 
     private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
     {
-        // EvenPerfectBitScanner only reaches this helper with divisors of the form 2 * exponent * k + 1 where exponent >= 5,
-        // so the small-prime return path is unreachable in production.
-        // if (prime <= 3UL)
-        // {
-        //     return prime == 3UL ? 2UL : 1UL;
-        // }
+        if (prime <= 3UL)
+        {
+            return prime == 3UL ? 2UL : 1UL;
+        }
 
         ulong phi = prime - 1UL;
 
-        // The by-divisor scanner currently defaults to the CPU heuristic device, but --mersenne-device may request GPU heuristics.
         if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder))
         {
             return gpuOrder;
@@ -214,19 +157,17 @@ internal static partial class PrimeOrderCalculator
 
         PartialFactorResult phiFactors = PartialFactor(phi, config);
 
-        ulong result;
-        // The by-divisor scan always yields phi = prime - 1 with at least the factor 2 captured during small-prime sieving,
-        // so PartialFactor never returns a null factor list on this call path.
-        // if (phiFactors.Factors is null)
-        // {
-        //     result = CalculateByFactorizationCpu(prime, divisorData, config);
-        //     phiFactors.Dispose();
-        //     return result;
-        // }
+		ulong result;
+		if (phiFactors.Factors is null)
+		{
+			result = FinishStrictlyCpu(prime, divisorData);
+			phiFactors.Dispose();
+			return result;
+		}
 
-        result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors);
-        phiFactors.Dispose();
-        return result;
+		result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors);
+		phiFactors.Dispose();
+		return result;
     }
 
     private static ulong RunHeuristicPipelineCpu(
@@ -237,52 +178,35 @@ internal static partial class PrimeOrderCalculator
         ulong phi,
         PartialFactorResult phiFactors)
     {
-        var cycleScope = SuppressDivisorCycleUsage(config.Mode == PrimeOrderMode.Heuristic && !IsGpuHeuristicDevice);
-        cycleScope.Enter();
-        try
+        if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData))
         {
-            if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData))
-            {
-                return phi;
-            }
+            return phi;
+        }
 
-            ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData);
-            candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData);
+        ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData);
+        candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData);
 
-            PartialFactorResult orderFactors = PartialFactor(candidateOrder, config);
-            if (orderFactors.Factors is null)
-            {
-                orderFactors.Dispose();
-                return candidateOrder;
-            }
-
-            // Reuse the computed factorization for both the confirmation check and the heuristic pass
-            // to avoid factoring the same order twice on the by-divisor CPU pipeline.
-            if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, ref orderFactors))
-            {
-                orderFactors.Dispose();
-                return candidateOrder;
-            }
-
-            if (config.Mode == PrimeOrderMode.Strict)
-            {
-                orderFactors.Dispose();
-                return CalculateByFactorizationCpu(prime, divisorData, config);
-            }
-
-            if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, ref orderFactors, phiFactors, out ulong order))
-            {
-                orderFactors.Dispose();
-                return order;
-            }
-
-            orderFactors.Dispose();
+        if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config))
+        {
             return candidateOrder;
         }
-        finally
+
+        if (config.Mode == PrimeOrderMode.Strict)
         {
-            cycleScope.Exit();
+            return FinishStrictlyCpu(prime, divisorData);
         }
+
+        if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order))
+        {
+            return order;
+        }
+
+        return candidateOrder;
+    }
+
+    private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData)
+    {
+        return CalculateByFactorizationCpu(prime, divisorData);
     }
 
     private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
@@ -317,52 +241,31 @@ internal static partial class PrimeOrderCalculator
         return order;
     }
 
-    [ThreadStatic]
-    private static ArrayPool<FactorEntry>? _factorPool;
-    private static ArrayPool<FactorEntry> FactorPool => _factorPool ??= ArrayPool<FactorEntry>.Create();
-
-    private const int FactorEntryStackThreshold = 32;
-
-    [ThreadStatic]
-    private static ArrayPool<ulong>? _ulongPool;
-    private static ArrayPool<ulong> ULongPool => _ulongPool ??= ArrayPool<ulong>.Create();
-
     private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
     {
-        ArrayPool<FactorEntry> factorPool = FactorPool;
-
         FactorEntry[]? tempArray = null;
-		ReadOnlySpan<FactorEntry> factorSpan = factors.Factors;
-		int length = factors.Count;
-		int capacity = length + 1;
-
-		Span<FactorEntry> buffer = capacity <= FactorEntryStackThreshold ? stackalloc FactorEntry[capacity] : default;
-		if (capacity > FactorEntryStackThreshold)
-		{
-			tempArray = factorPool.Rent(capacity);
-			buffer = tempArray.AsSpan(0, capacity);
-		}
-
         try
         {
+            ReadOnlySpan<FactorEntry> factorSpan = factors.Factors;
+            int length = factors.Count;
+            tempArray = ArrayPool<FactorEntry>.Shared.Rent(length + 1);
+            Span<FactorEntry> buffer = tempArray;
             factorSpan.CopyTo(buffer);
 
             bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor);
             // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None);
             if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime)
             {
-                buffer[length] = new FactorEntry(factors.Cofactor);
+                buffer[length] = new FactorEntry(factors.Cofactor, 1);
                 length++;
             }
 
-            buffer = buffer.Slice(0, length); // Reuse the same span restricted to populated entries.
-            buffer.Sort(static (a, b) => a.Value.CompareTo(b.Value));
+            buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value));
 
-            for (int i = 0; i < buffer.Length; i++)
+            for (int i = 0; i < length; i++)
             {
-                FactorEntry entry = buffer[i];
-                ulong primeFactor = entry.Value;
-                int exponent = entry.Exponent;
+                ulong primeFactor = buffer[i].Value;
+                int exponent = buffer[i].Exponent;
                 for (int iteration = 0; iteration < exponent; iteration++)
                 {
                     if ((order % primeFactor) != 0UL)
@@ -387,18 +290,17 @@ internal static partial class PrimeOrderCalculator
         {
             if (tempArray is not null)
             {
-                factorPool.Return(tempArray, clearArray: false);
+                ArrayPool<FactorEntry>.Shared.Return(tempArray, clearArray: false);
             }
         }
     }
 
-    private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, ref PartialFactorResult orderFactors)
+    private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
     {
-        // The candidate order is always a positive divisor of phi in this flow, so zero never appears here.
-        // if (order == 0UL)
-        // {
-        //     return false;
-        // }
+        if (order == 0UL)
+        {
+            return false;
+        }
 
         // 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)");
@@ -407,57 +309,66 @@ internal static partial class PrimeOrderCalculator
             return false;
         }
 
-        // PartialFactor returns a null factor span when it cannot peel off any divisors, so keep the guard for composite candidate orders.
-        if (orderFactors.Factors is null)
-        {
-            return false;
-        }
+        // DebugLog("Partial factoring order");
 
-        if (!orderFactors.FullyFactored)
+        // TODO: Do we do partial factoring of order multiple times?
+        PartialFactorResult factorization = PartialFactor(order, config);
+        try
         {
-            // A leftover cofactor of 0 or 1 indicates the partial factorization produced no useful prime divisor.
-            if (orderFactors.Cofactor <= 1UL)
+            if (factorization.Factors is null)
             {
                 return false;
             }
 
-            // DebugLog("Cofactor > 1. Testing primality of cofactor");
-            bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor);
-            // bool isPrime = PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None);
-            if (!isPrime)
+            if (!factorization.FullyFactored)
             {
-                return false;
-            }
-
-            // DebugLog("Adding cofactor as prime factor");
-            PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
-            orderFactors.Dispose();
-            orderFactors = extended;
-        }
-
-        ReadOnlySpan<FactorEntry> span = orderFactors.Factors!;
-        // DebugLog("Verifying prime-power reductions");
-        int length = orderFactors.Count;
-        for (int i = 0; i < length; i++)
-        {
-            ulong primeFactor = span[i].Value;
-            ulong reduced = order;
-            for (int iteration = 0; iteration < span[i].Exponent; iteration++)
-            {
-                if ((reduced % primeFactor) != 0UL)
-                {
-                    break;
-                }
-
-                reduced /= primeFactor;
-                if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+                if (factorization.Cofactor <= 1UL)
                 {
                     return false;
                 }
-            }
-        }
 
-        return true;
+                // DebugLog("Cofactor > 1. Testing primality of cofactor");
+                bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
+                // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
+                if (!isPrime)
+                {
+                    return false;
+                }
+
+                // DebugLog("Adding cofactor as prime factor");
+                PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
+                factorization.Dispose();
+                factorization = extended;
+            }
+
+            ReadOnlySpan<FactorEntry> span = factorization.Factors;
+            // DebugLog("Verifying prime-power reductions");
+            int length = factorization.Count;
+            for (int i = 0; i < length; i++)
+            {
+                ulong primeFactor = span[i].Value;
+                ulong reduced = order;
+                for (int iteration = 0; iteration < span[i].Exponent; iteration++)
+                {
+                    if ((reduced % primeFactor) != 0UL)
+                    {
+                        break;
+                    }
+
+                    reduced /= primeFactor;
+                    if (Pow2EqualsOneCpu(reduced, prime, divisorData))
+                    {
+                        return false;
+                    }
+                }
+            }
+
+            return true;
+        }
+        finally
+        {
+            factorization.Dispose();
+        }
     }
 
     private static bool TryHeuristicFinishCpu(
@@ -466,164 +377,169 @@ internal static partial class PrimeOrderCalculator
         ulong? previousOrder,
         in MontgomeryDivisorData divisorData,
         in PrimeOrderSearchConfig config,
-        ref PartialFactorResult orderFactors,
         PartialFactorResult phiFactors,
         out ulong result)
     {
         result = 0UL;
-        // The candidate order stays at least 2 for odd divisors above 3, so this guard never triggers in the EvenPerfectBitScanner flow.
-        // if (order <= 1UL)
-        // {
-        //     return false;
-        // }
-
-        // PartialFactor returns a null factor span when it cannot peel off any divisors, so the heuristic pass must bail in that case.
-        if (orderFactors.Factors is null)
+        if (order <= 1UL)
         {
             return false;
         }
 
-        if (!orderFactors.FullyFactored)
+        // TODO: Do we do partial factoring of order multiple times?
+        // DebugLog("Trying heuristic. Partial factoring order");
+        PartialFactorResult orderFactors = PartialFactor(order, config);
+        try
         {
-            // The partial factoring helper reports a cofactor of 0 or 1 when nothing usable remains, so short-circuit the heuristic stage.
-            if (orderFactors.Cofactor <= 1UL)
+            if (orderFactors.Factors is null)
             {
                 return false;
             }
 
-            var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor);
-            // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None))
-            if (!isPrime)
+            if (!orderFactors.FullyFactored)
             {
-                return false;
-            }
-
-            PartialFactorResult extended = 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");
-        BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
-        if (candidates.Count == 0)
-        {
-            return false;
-        }
-
-        // DebugLog("Sorting candidates");
-        SortCandidates(prime, previousOrder, candidates);
-
-        int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
-        int powUsed = 0;
-        int candidateCount = candidates.Count;
-        bool allowGpuBatch = true;
-        Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
-
-        // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
-        int index = 0;
-        const int MaxGpuBatchSize = 1024;
-        const int StackGpuBatchSize = 1024;
-        Span<ulong> stackGpuRemainders = stackalloc ulong[StackGpuBatchSize];
-        ArrayPool<ulong> uLongPool = ULongPool;
-        while (index < candidateCount && powUsed < powBudget)
-        {
-            int remaining = candidateCount - index;
-            int budgetRemaining = powBudget - powUsed;
-            int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize));
-            if (batchSize <= 0)
-            {
-                break;
-            }
-
-            ReadOnlySpan<ulong> batch = candidateSpan.Slice(index, batchSize);
-            ulong[]? gpuPool = null;
-            Span<ulong> pooledGpuRemainders = default;
-            bool gpuSuccess = false;
-            bool gpuStackRemainders = false;
-            GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
-
-            if (allowGpuBatch && IsGpuPow2Allowed)
-            {
-                if (batchSize <= StackGpuBatchSize)
+                if (orderFactors.Cofactor <= 1UL)
                 {
-                    Span<ulong> localRemainders = stackGpuRemainders.Slice(0, batchSize);
-                    status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData);
-                    if (status == GpuPow2ModStatus.Success)
-                    {
-                        gpuSuccess = true;
-                        gpuStackRemainders = true;
-                    }
+                    return false;
                 }
-                else
+
+                var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor);
+                // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None))
+                if (!isPrime)
                 {
-                    gpuPool = uLongPool.Rent(batchSize);
-                    Span<ulong> pooledRemainders = gpuPool.AsSpan(0, batchSize);
-                    status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData);
-                    if (status == GpuPow2ModStatus.Success)
+                    return false;
+                }
+
+                PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
+                orderFactors.Dispose();
+                orderFactors = extended;
+            }
+
+            int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4;
+            List<ulong> candidates = new(capacity);
+            FactorEntry[] factorArray = orderFactors.Factors!;
+            // DebugLog("Building candidates list");
+            BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
+            if (candidates.Count == 0)
+            {
+                return false;
+            }
+
+            // DebugLog("Sorting candidates");
+            SortCandidates(prime, previousOrder, candidates);
+
+            int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
+            int powUsed = 0;
+            int candidateCount = candidates.Count;
+            bool allowGpuBatch = true;
+            Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
+
+            // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
+            int index = 0;
+            const int MaxGpuBatchSize = 1024;
+            const int StackGpuBatchSize = 1024;
+            Span<ulong> stackGpuRemainders = stackalloc ulong[StackGpuBatchSize];
+            while (index < candidateCount && powUsed < powBudget)
+            {
+                int remaining = candidateCount - index;
+                int budgetRemaining = powBudget - powUsed;
+                int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize));
+                if (batchSize <= 0)
+                {
+                    break;
+                }
+
+                ReadOnlySpan<ulong> batch = candidateSpan.Slice(index, batchSize);
+                ulong[]? gpuPool = null;
+                Span<ulong> pooledGpuRemainders = default;
+                bool gpuSuccess = false;
+                bool gpuStackRemainders = false;
+                GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
+
+                if (allowGpuBatch && IsGpuPow2Allowed)
+                {
+                    if (batchSize <= StackGpuBatchSize)
                     {
-                        pooledGpuRemainders = pooledRemainders;
-                        gpuSuccess = true;
+                        Span<ulong> localRemainders = stackGpuRemainders.Slice(0, batchSize);
+                        status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData);
+                        if (status == GpuPow2ModStatus.Success)
+                        {
+                            gpuSuccess = true;
+                            gpuStackRemainders = true;
+                        }
                     }
                     else
                     {
-                        uLongPool.Return(gpuPool, clearArray: false);
-                        gpuPool = null;
+                        gpuPool = ArrayPool<ulong>.Shared.Rent(batchSize);
+                        Span<ulong> pooledRemainders = gpuPool.AsSpan(0, batchSize);
+                        status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData);
+                        if (status == GpuPow2ModStatus.Success)
+                        {
+                            pooledGpuRemainders = pooledRemainders;
+                            gpuSuccess = true;
+                        }
+                        else
+                        {
+                            ArrayPool<ulong>.Shared.Return(gpuPool, clearArray: false);
+                            gpuPool = null;
+                        }
+                    }
+
+                    if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable))
+                    {
+                        allowGpuBatch = false;
                     }
                 }
 
-                if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable))
+                for (int i = 0; i < batchSize && powUsed < powBudget; i++)
                 {
-                    allowGpuBatch = false;
-                }
-            }
+                    ulong candidate = batch[i];
+                    powUsed++;
 
-            for (int i = 0; i < batchSize && powUsed < powBudget; i++)
-            {
-                ulong candidate = batch[i];
-                powUsed++;
+                    bool equalsOne;
+                    if (gpuSuccess)
+                    {
+                        ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i];
+                        equalsOne = remainderValue == 1UL;
+                    }
+                    else
+                    {
+                        equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData);
+                    }
+                    if (!equalsOne)
+                    {
+                        continue;
+                    }
 
-                bool equalsOne;
-                if (gpuSuccess)
-                {
-                    ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i];
-                    equalsOne = remainderValue == 1UL;
-                }
-                else
-                {
-                    equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData);
-                }
-                if (!equalsOne)
-                {
-                    continue;
-                }
+                    if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget))
+                    {
+                        continue;
+                    }
 
-                if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget))
-                {
-                    continue;
+                    if (gpuPool is not null)
+                    {
+                        ArrayPool<ulong>.Shared.Return(gpuPool, clearArray: false);
+                    }
+
+                    result = candidate;
+                    return true;
                 }
 
                 if (gpuPool is not null)
                 {
-                    uLongPool.Return(gpuPool, clearArray: false);
+                    ArrayPool<ulong>.Shared.Return(gpuPool, clearArray: false);
                 }
 
-                result = candidate;
-                return true;
+                index += batchSize;
             }
 
-            if (gpuPool is not null)
-            {
-                uLongPool.Return(gpuPool, clearArray: false);
-            }
-
-            index += batchSize;
+            // DebugLog("No candidate confirmed");
+            return false;
+        }
+        finally
+        {
+            orderFactors.Dispose();
         }
-
-        // DebugLog("No candidate confirmed");
-        return false;
     }
 
     private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates)
@@ -857,36 +773,6 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool TryGetDivisorCycle(ulong prime, in MontgomeryDivisorData divisorData, out ulong cycleLength)
-    {
-        if (prime <= 1UL || IsDivisorCycleSuppressed)
-        {
-            cycleLength = 0UL;
-            return false;
-        }
-
-        if (s_hasCachedDivisorCycle && s_cachedDivisorCycleModulus == prime)
-        {
-            cycleLength = s_cachedDivisorCycleLength;
-            return cycleLength != 0UL;
-        }
-
-        var suppression = SuppressDivisorCycleUsage(true);
-        suppression.Enter();
-        try
-        {
-            cycleLength = MersenneDivisorCycles.CalculateCycleLength(prime, divisorData);
-            s_cachedDivisorCycleModulus = prime;
-            s_cachedDivisorCycleLength = cycleLength;
-            s_hasCachedDivisorCycle = true;
-            return cycleLength != 0UL;
-        }
-        finally
-        {
-            suppression.Exit();
-        }
-    }
-
     private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
     {
         if (IsGpuPow2Allowed)
@@ -898,11 +784,6 @@ internal static partial class PrimeOrderCalculator
             }
         }
 
-        if (TryGetDivisorCycle(prime, divisorData, out ulong cycleLength) && cycleLength != 0UL)
-        {
-            return exponent.Pow2MontgomeryModWithCycleCpu(cycleLength, divisorData) == 1UL;
-        }
-
         return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
     }
 
@@ -1072,10 +953,7 @@ internal static partial class PrimeOrderCalculator
                     return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
                 }
 
-                int dictionaryCount = counts.Count;
-                FactorEntry[] factors = dictionaryCount >= PerfectNumberConstants.PooledArrayThreshold
-                    ? ThreadLocalArrayPool<FactorEntry>.Shared.Rent(dictionaryCount)
-                    : new FactorEntry[dictionaryCount];
+                FactorEntry[] factors = ArrayPool<FactorEntry>.Shared.Rent(counts.Count);
                 int index = 0;
                 foreach (KeyValuePair<ulong, int> entry in counts)
                 {
@@ -1109,9 +987,7 @@ internal static partial class PrimeOrderCalculator
                 return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
             }
 
-            FactorEntry[] array = actualCount >= PerfectNumberConstants.PooledArrayThreshold
-                ? ThreadLocalArrayPool<FactorEntry>.Shared.Rent(actualCount)
-                : new FactorEntry[actualCount];
+            FactorEntry[] array = ArrayPool<FactorEntry>.Shared.Rent(actualCount);
             int arrayIndex = 0;
             for (int i = 0; i < factorCount; i++)
             {
@@ -1288,6 +1164,10 @@ internal static partial class PrimeOrderCalculator
             return true;
         }
 
+        PollardRhoMontgomeryReducer reducer = PollardRhoMontgomeryReducer.Create(n);
+        ulong diff = 0UL;
+        ulong d = 0UL;
+
         while (true)
         {
             if (stopwatch.ElapsedTicks > budgetTicks)
@@ -1296,9 +1176,12 @@ internal static partial class PrimeOrderCalculator
             }
 
             ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
-            ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
-            ulong y = x;
-            ulong d = 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)
             {
@@ -1307,10 +1190,15 @@ internal static partial class PrimeOrderCalculator
                     return false;
                 }
 
-                x = AdvancePolynomial(x, c, n);
-                y = AdvancePolynomial(y, c, n);
-                y = AdvancePolynomial(y, c, n);
-                ulong diff = x > y ? x - y : y - x;
+                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);
             }
 
@@ -1322,16 +1210,6 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus)
-    {
-        ref MutableUInt128 polynomial = ref s_pollardRhoPolynomialBuffer;
-        polynomial.Set(x);
-        polynomial.Mul(x);
-        polynomial.Add(c);
-        return polynomial.Mod(modulus);
-    }
-
     private static ulong BinaryGcd(ulong a, ulong b)
     {
         if (a == 0UL)
@@ -1453,19 +1331,8 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
+    private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData)
     {
-        if (!IsGpuHeuristicDevice && IsGpuCycleAllowed)
-        {
-            PrimeOrderSearchConfig strictConfig = config.Mode == PrimeOrderMode.Strict
-                ? config
-                : new PrimeOrderSearchConfig(config.SmallFactorLimit, config.PollardRhoMilliseconds, config.MaxPowChecks, PrimeOrderMode.Strict);
-            if (PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder: null, strictConfig, divisorData, out ulong gpuOrder))
-            {
-                return gpuOrder;
-            }
-        }
-
         ulong phi = prime - 1UL;
         Dictionary<ulong, int> counts = new(capacity: 8);
         FactorCompletely(phi, counts);
@@ -1529,19 +1396,31 @@ internal static partial class PrimeOrderCalculator
             return 2UL;
         }
 
+        PollardRhoMontgomeryReducer reducer = PollardRhoMontgomeryReducer.Create(n);
+        ulong diff = 0UL;
+        ulong d = 0UL;
+
         while (true)
         {
             ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
-            ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
-            ulong y = x;
-            ulong d = 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)
             {
-                x = AdvancePolynomial(x, c, n);
-                y = AdvancePolynomial(y, c, n);
-                y = AdvancePolynomial(y, c, n);
-                ulong diff = x > y ? x - y : y - x;
+                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);
             }
 

PerfectNumbers.Core/PrimeTestTimeLimit.cs

@@ -1,72 +0,0 @@
-using System;
-using System.Diagnostics;
-
-namespace PerfectNumbers.Core;
-
-internal readonly struct PrimeTestTimeLimit
-{
-    private readonly long _deadlineTimestamp;
-
-    private PrimeTestTimeLimit(bool isActive, long deadlineTimestamp)
-    {
-        IsActive = isActive;
-        _deadlineTimestamp = deadlineTimestamp;
-    }
-
-    internal bool IsActive { get; }
-
-    internal static bool TryCreate(TimeSpan? timeLimit, out PrimeTestTimeLimit limit)
-    {
-        if (!timeLimit.HasValue)
-        {
-            limit = new PrimeTestTimeLimit(isActive: false, deadlineTimestamp: long.MaxValue);
-            return true;
-        }
-
-        TimeSpan limitValue = timeLimit.GetValueOrDefault();
-        if (limitValue <= TimeSpan.Zero)
-        {
-            limit = new PrimeTestTimeLimit(isActive: true, deadlineTimestamp: long.MinValue);
-            return false;
-        }
-
-        double seconds = limitValue.TotalSeconds;
-        if (double.IsInfinity(seconds))
-        {
-            limit = new PrimeTestTimeLimit(isActive: true, deadlineTimestamp: long.MaxValue);
-            return true;
-        }
-
-        double scaled = seconds * Stopwatch.Frequency;
-        long limitTicks = scaled >= long.MaxValue ? long.MaxValue : (long)scaled;
-        long startTimestamp = Stopwatch.GetTimestamp();
-        long deadline = startTimestamp <= long.MaxValue - limitTicks ? startTimestamp + limitTicks : long.MaxValue;
-        if (deadline < startTimestamp)
-        {
-            deadline = long.MaxValue;
-        }
-
-        limit = new PrimeTestTimeLimit(isActive: true, deadlineTimestamp: deadline);
-        return true;
-    }
-
-    internal bool HasExpired()
-    {
-        if (!IsActive)
-        {
-            return false;
-        }
-
-        if (_deadlineTimestamp == long.MaxValue)
-        {
-            return false;
-        }
-
-        if (_deadlineTimestamp == long.MinValue)
-        {
-            return true;
-        }
-
-        return Stopwatch.GetTimestamp() >= _deadlineTimestamp;
-    }
-}

PerfectNumbers.Core/ReciprocalMath.cs

@@ -0,0 +1,125 @@
+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/ThreadLocalArrayPool.cs

@@ -1,27 +0,0 @@
-using System;
-using System.Buffers;
-using System.Runtime.CompilerServices;
-
-namespace PerfectNumbers.Core;
-
-internal static class ThreadLocalArrayPool<T>
-{
-    [ThreadStatic]
-    private static ArrayPool<T>? s_threadPool;
-
-    public static ArrayPool<T> Shared
-    {
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        get
-        {
-            ArrayPool<T>? pool = s_threadPool;
-            if (pool is null)
-            {
-                pool = ArrayPool<T>.Create();
-                s_threadPool = pool;
-            }
-
-            return pool;
-        }
-    }
-}

PerfectNumbers.Core/ULongExtensions.cs

@@ -310,37 +310,6 @@ public static class ULongExtensions
         return state.MulModWithNativeModulo(b, modulus);
     }
 
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus)
-    {
-        if (modulus <= 1UL)
-        {
-            return 0UL;
-        }
-
-        ulong remainingExponent = exponent;
-        ulong result = 1UL % modulus;
-        ulong baseValue = 2UL % modulus;
-
-        while (remainingExponent != 0UL)
-        {
-            if ((remainingExponent & 1UL) != 0UL)
-            {
-                result = result.MulMod64GpuCompatible(baseValue, modulus);
-            }
-
-            remainingExponent >>= 1;
-            if (remainingExponent == 0UL)
-            {
-                break;
-            }
-
-            baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus);
-        }
-
-        return result;
-    }
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime)
     {
@@ -377,63 +346,48 @@ public static class ULongExtensions
 
         ulong result = divisor.MontgomeryOne;
         ulong nPrime = divisor.NPrime;
-
-        Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : default;
-        ulong[]? pooledOddPowers = null;
-
-		if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold)
-        {
-            oddPowerStorage = new ulong[oddPowerCount];
-        }
-        else
-        {
-            pooledOddPowers = ArrayPool<ulong>.Shared.Rent(oddPowerCount);
-            oddPowerStorage = pooledOddPowers;
-        }
-
-        Span<ulong> oddPowers = oddPowerStorage[..oddPowerCount];
-        InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers);
+        ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount);
 
         int index = bitLength - 1;
-        while (index >= 0)
-        {
-            if (((exponent >> index) & 1UL) == 0UL)
-            {
-                result = result.MontgomeryMultiply(result, modulus, nPrime);
-                index--;
-                continue;
-            }
+		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;
-            }
+			int windowStart = index - windowSize + 1;
+			if (windowStart < 0)
+			{
+				windowStart = 0;
+			}
 
-            while (((exponent >> windowStart) & 1UL) == 0UL)
-            {
-                windowStart++;
-            }
+			while (((exponent >> windowStart) & 1UL) == 0UL)
+			{
+				windowStart++;
+			}
 
-            int windowLength = index - windowStart + 1;
-            for (int square = 0; square < windowLength; square++)
-            {
-                result = result.MontgomeryMultiply(result, modulus, nPrime);
-            }
+			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);
+			ulong mask = (1UL << windowLength) - 1UL;
+			ulong windowValue = (exponent >> windowStart) & mask;
+			int tableIndex = (int)((windowValue - 1UL) >> 1);
+			ulong multiplier = oddPowers[tableIndex];
+			result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
 
-            index = windowStart - 1;
-        }
-
-        if (pooledOddPowers is not null)
-        {
-            ArrayPool<ulong>.Shared.Return(pooledOddPowers);
-        }
+			index = windowStart - 1;
+		}
+		
+		if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold)
+		{
+			ArrayPool<ulong>.Shared.Return(oddPowers);
+		}
 
         if (keepMontgomery)
         {
@@ -579,34 +533,26 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> destination)
+    private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
     {
-        int oddPowerCount = destination.Length;
-        if (oddPowerCount == 0)
-        {
-            return;
-        }
-
-        bool computedOnGpu = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
-            && MontgomeryOddPowerGpu.TryCompute(divisor, oddPowerCount, destination);
-        if (computedOnGpu)
-        {
-            return;
-        }
-
-        destination[0] = divisor.MontgomeryTwo;
+		ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold ? new ulong[oddPowerCount] : ArrayPool<ulong>.Shared.Rent(oddPowerCount);
+		
+        oddPowers[0] = divisor.MontgomeryTwo;
         if (oddPowerCount == 1)
         {
-            return;
+            return oddPowers;
         }
 
-        ulong square = divisor.MontgomeryTwoSquared;
-        ulong current = destination[0];
+        ulong previous,
+			  square = divisor.MontgomeryTwoSquared;
+
         for (int i = 1; i < oddPowerCount; i++)
         {
-            current = current.MontgomeryMultiply(square, modulus, nPrime);
-            destination[i] = current;
+            previous = oddPowers[i - 1];
+            oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
         }
+
+        return oddPowers;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -636,11 +582,11 @@ public static class ULongExtensions
         return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false);
     }
 
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
-    {
-        return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false);
-    }
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+	{
+		return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false);
+	}
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
@@ -649,12 +595,12 @@ public static class ULongExtensions
         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 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

@@ -0,0 +1,196 @@
+# 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`.
+

gpu-montgomery-usage.md

@@ -1,18 +0,0 @@
-# GPU Modular Arithmetic in the `--mersenne=bydivisor` Path
-
-The GPU implementation now routes single-use exponent checks through standard modular arithmetic while keeping Montgomery reduction only where repeated multiplications reuse the same modulus.
-
-## Kernels using standard modular arithmetic
-
-- `CheckKernel` in `MersenneNumberDivisorByDivisorGpuTester` evaluates each divisor candidate with `Pow2ModBinaryGpu` and compares the result directly against 1. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L561-L645】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
-- `ComputeExponentKernel` in the same tester precomputes batched residues with the plain helper so CPU lookups receive standard-domain results. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L780-L791】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
-- `ComputeResiduesKernel` inside `MersenneNumberDivisorResidueGpuEvaluator` now mirrors that logic when building residue tables for mixed CPU/GPU scans. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs†L102-L114】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
-- The order-search helpers in `PrimeOrderGpuHeuristics` (`TrySpecialMaxKernel`, `ExponentLoweringKernel`, and `Pow2EqualsOneKernel`) also rely on `Pow2ModBinaryGpu` for their equality checks. 【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L798-L866】【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L1596-L1598】
-
-## Kernels still using Montgomery arithmetic
-
-- `MontgomeryOddPowerKernelScan` in `GpuKernelPool` continues to build windowed odd-power tables with `MontgomeryMultiply`, which amortizes the Montgomery setup over many multiplications per modulus. 【F:PerfectNumbers.Core/Gpu/GpuKernelPool.cs†L200-L237】【F:PerfectNumbers.Core/ULongExtensions.cs†L532-L559】
-
-## Shared helper
-
-- `Pow2ModBinaryGpu` in `ULongExtensions` implements the GPU-friendly binary exponentiation loop that powers the standard modular path across these kernels. 【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】

todo_list.txt

@@ -123,7 +123,6 @@ 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.
@@ -267,7 +266,6 @@ 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.

tools/HeuristicMultiplyShiftTiming/HeuristicMultiplyShiftTiming.csproj

@@ -8,7 +8,4 @@
   <ItemGroup>
     <ProjectReference Include="../../PerfectNumbers.Core/PerfectNumbers.Core.csproj" />
   </ItemGroup>
-  <ItemGroup>
-    <Compile Remove=".temp.cs" />
-  </ItemGroup>
-</Project>
+</Project>