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

Cache Pollard rho reciprocals for faster Montgomery steps

AGENTS.md

@@ -106,16 +106,11 @@ There is only one agents file - this one. Don't search for another one.
 
 Always use indentation with four space characters in source code files, even if a file currently uses tabs or appears to rely on eight-space indentation. Increase indentation level by exactly one indentation unit (four spaces) per nested block. When parameters or long expressions wrap, indent the continuation with two indentation units for clarity.
 Explicitly prefer reusing existing variables once their previous values are no longer required instead of declaring new locals of the same type. Call out the variable reuse in comments when its meaning changes so the intent stays clear and register pressure remains minimal.
-Do not introduce try/finally blocks solely to release pooled resources. Restructure the code so that every return path releases the resources without relying on finally unless the control flow would become unclear or unmaintainable.
 
 Keep the divisor cycle cache limited to a single block in memory (only the on-disk snapshot). When a lookup misses, compute the required cycle on the configured device and discard the transient result instead of scheduling background generation or caching additional blocks.
 
 All code, comments, commit messages, branch names and PR descriptions must be written in American English.
 
-- When modifying existing files, preserve all relevant TODO notes and explanatory comments unless they are no longer valid.
-- Prefer defining data holders as structs or classes with public fields over auto-properties when performance on the hot path is a concern.
-- Do not add or keep methods that simply wrap another method call without adding any logic; inline those calls or extend the wrapper with meaningful work.
-
 ## Performance implementation checklist
 
 - Prefer dedicated kernels over flag-controlled or parameterized kernels. Avoid conditional branches in kernels and tight loops by creating separate kernels when necessary.

EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core.Gpu;

EvenPerfectBitScanner.Benchmarks/GpuUInt128AddBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;

EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs

@@ -1,3 +1,5 @@
+using System;
+using System.Collections.Generic;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core.Gpu;

EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;

EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;

EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModByLimbBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core.Gpu;

EvenPerfectBitScanner.Benchmarks/GpuUInt128NativeModuloBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core.Gpu;

EvenPerfectBitScanner.Benchmarks/GpuUInt128ScalarMulModBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core.Gpu;

EvenPerfectBitScanner.Benchmarks/GpuUInt128SubModScalarBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core.Gpu;

EvenPerfectBitScanner.Benchmarks/HeuristicMultiplyShiftBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core;

EvenPerfectBitScanner.Benchmarks/MersenneDivisorCycleLengthGpuBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;

EvenPerfectBitScanner.Benchmarks/Mod11ULongBenchmarks.cs

@@ -1,6 +1,7 @@
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
+using PerfectNumbers.Core;
 
 namespace EvenPerfectBitScanner.Benchmarks;
 

EvenPerfectBitScanner.Benchmarks/Mul64Benchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;

EvenPerfectBitScanner.Benchmarks/MulHighBenchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;

EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs

@@ -1,3 +1,4 @@
+using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs

@@ -1,7 +1,9 @@
+using System;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using System.Numerics;
 using PerfectNumbers.Core;
+using PerfectNumbers.Core.Gpu;
 
 namespace EvenPerfectBitScanner.Benchmarks;
 

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs

@@ -1,4 +1,6 @@
+using System;
 using System.Buffers.Binary;
+using System.Threading;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core;
@@ -375,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/Pow2MontgomeryModUInt128Benchmarks.cs

@@ -1,4 +1,6 @@
+using System;
 using System.Buffers.Binary;
+using System.Numerics;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core;

EvenPerfectBitScanner.Benchmarks/Pow2OddPowerTableBenchmarks.cs

@@ -1,4 +1,8 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
 using System.Runtime.CompilerServices;
+using System.Threading.Tasks;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using PerfectNumbers.Core.Gpu;

EvenPerfectBitScanner.Benchmarks/Program.cs

@@ -1,3 +1,5 @@
+using BenchmarkDotNet.Configs;
+using BenchmarkDotNet.Running;
 using PerfectNumbers.Core.Gpu;
 
 namespace EvenPerfectBitScanner.Benchmarks;

EvenPerfectBitScanner/Program.cs

@@ -1,4 +1,5 @@
-using System.Diagnostics;
+using System.Buffers;
+using System.Collections.Generic;
 using System.Globalization;
 using System.Numerics;
 using System.Runtime.CompilerServices;
@@ -44,11 +45,6 @@ internal static class Program
     private static readonly Optimized PrimeIterator = new();
     private static ulong _byDivisorStartPrime;
 
-    private static bool _testMode;
-    private static int _testTargetPrimeCount;
-    private static int _testProcessedPrimeCount;
-    private const string TestTimeFileName = "even_perfect_test_time.txt";
-
     [ThreadStatic]
     private static bool _lastCompositeP;
 
@@ -85,8 +81,6 @@ internal static class Program
         int scanBatchSize = 2_097_152, sliceSize = 32;
         ulong residueKMax = 5_000_000UL;
         string filterFile = string.Empty;
-        UInt128 maxPrimeLimit = UInt128.MaxValue;
-        bool maxPrimeConfigured = false;
         string cyclesPath = DefaultCyclesPath;
         int cyclesBatchSize = 512;
         bool continueCyclesGeneration = false;
@@ -99,10 +93,6 @@ internal static class Program
         ulong remainder = 0UL;
         bool startPrimeProvided = false;
 
-        _testMode = false;
-        _testTargetPrimeCount = 0;
-        _testProcessedPrimeCount = 0;
-
         // NTT backend selection (GPU): reference vs staged
         for (; argIndex < args.Length; argIndex++)
         {
@@ -119,20 +109,6 @@ internal static class Program
                 currentP = ulong.Parse(arg.AsSpan(arg.IndexOf('=') + 1));
                 startPrimeProvided = true;
             }
-            else if (arg.StartsWith("--max-prime=", StringComparison.OrdinalIgnoreCase))
-            {
-                ReadOnlySpan<char> value = arg.AsSpan(arg.IndexOf('=') + 1);
-                if (UInt128.TryParse(value, NumberStyles.None, CultureInfo.InvariantCulture, out UInt128 parsedMaxPrime))
-                {
-                    maxPrimeLimit = parsedMaxPrime;
-                    maxPrimeConfigured = true;
-                }
-                else
-                {
-                    Console.WriteLine("Invalid value for --max-prime.");
-                    return;
-                }
-            }
             else if (arg.Equals("--increment=bit", StringComparison.OrdinalIgnoreCase))
             {
                 useBitTransform = true;
@@ -209,10 +185,6 @@ internal static class Program
             {
                 useModuloWorkaround = true;
             }
-            else if (arg.Equals("--test", StringComparison.OrdinalIgnoreCase))
-            {
-                _testMode = true;
-            }
             else if (arg.Equals("--use-order", StringComparison.OrdinalIgnoreCase))
             {
                 useOrder = true;
@@ -405,19 +377,6 @@ internal static class Program
             }
         }
 
-        int testPrimeCandidateLimit = 0;
-
-        if (_testMode)
-        {
-            currentP = 31UL;
-            remainder = currentP % 6UL;
-            filterFile = string.Empty;
-            maxPrimeLimit = UInt128.MaxValue;
-            maxPrimeConfigured = false;
-            startPrimeProvided = true;
-            testPrimeCandidateLimit = Math.Max(1, threadCount) * 3;
-        }
-
         if (useByDivisor)
         {
             _byDivisorStartPrime = startPrimeProvided ? currentP : 0UL;
@@ -427,7 +386,7 @@ internal static class Program
             _byDivisorStartPrime = 0UL;
         }
 
-        if (useByDivisor && !_testMode && string.IsNullOrEmpty(filterFile))
+        if (useByDivisor && string.IsNullOrEmpty(filterFile))
         {
             Console.WriteLine("--mersenne=bydivisor requires --filter-p=<path>.");
             return;
@@ -645,102 +604,51 @@ internal static class Program
             File.WriteAllText(ResultsFileName, $"p,searchedMersenne,detailedCheck,passedAllTests{Environment.NewLine}");
         }
 
-        bool useFilter = !_testMode && !string.IsNullOrEmpty(filterFile);
+        bool useFilter = !string.IsNullOrEmpty(filterFile);
         HashSet<ulong> filter = [];
         List<ulong> byDivisorCandidates = [];
         ulong maxP = 0UL;
         ulong[] localFilter = Array.Empty<ulong>();
         int filterCount = 0;
-
-        if (_useByDivisorMode)
+        if (useFilter)
         {
-            if (_testMode)
+            Console.WriteLine("Loading filter...");
+            if (_useByDivisorMode)
             {
-                byDivisorCandidates = BuildTestPrimeCandidates(testPrimeCandidateLimit);
+                byDivisorCandidates = LoadByDivisorCandidates(filterFile);
             }
             else
             {
-                Console.WriteLine("Loading filter...");
-                int skippedByLimit = 0;
-                byDivisorCandidates = LoadByDivisorCandidates(filterFile, maxPrimeLimit, maxPrimeConfigured, out skippedByLimit);
-                if (maxPrimeConfigured)
+                // TODO: Rent this filter buffer from ArrayPool<ulong> and reuse it across reload batches so
+                // we do not allocate fresh arrays while replaying large result filters.
+                localFilter = new ulong[1024];
+                filterCount = 0;
+                LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) =>
                 {
-                    if (skippedByLimit > 0)
+                    if (passedAllTests)
                     {
-                        Console.WriteLine($"Skipped {skippedByLimit} by-divisor candidate(s) above --max-prime limit.");
+                        localFilter[filterCount++] = p;
+                        if (p > maxP)
+                        {
+                            maxP = p;
+                        }
+
+                        if (filterCount == 1024)
+                        {
+                            filter.AddRange(localFilter[..filterCount]);
+                            filterCount = 0;
+                            Console.WriteLine($"Added {p}");
+                        }
                     }
+                });
 
-                    if (byDivisorCandidates.Count == 0)
-                    {
-                        Console.WriteLine("No by-divisor candidates fall within the --max-prime limit.");
-                        return;
-                    }
-                }
-            }
-        }
-        else if (useFilter)
-        {
-            Console.WriteLine("Loading filter...");
-            // TODO: Rent this filter buffer from ArrayPool<ulong> and reuse it across reload batches so
-            // we do not allocate fresh arrays while replaying large result filters.
-            localFilter = new ulong[1024];
-            filterCount = 0;
-            bool addedWithinLimit = false;
-            int skippedByLimit = 0;
-            LoadResultsFile(filterFile, (p, detailedCheck, passedAllTests) =>
-            {
-                if (!passedAllTests)
-                {
-                    return;
-                }
-
-                if (maxPrimeConfigured && (UInt128)p > maxPrimeLimit)
-                {
-                    skippedByLimit++;
-                    return;
-                }
-
-                localFilter[filterCount++] = p;
-                addedWithinLimit = true;
-                if (p > maxP)
-                {
-                    maxP = p;
-                }
-
-                if (filterCount == 1024)
+                if (filterCount > 0)
                 {
                     filter.AddRange(localFilter[..filterCount]);
-                    filterCount = 0;
-                    Console.WriteLine($"Added {p}");
-                }
-            });
-
-            if (filterCount > 0)
-            {
-                filter.AddRange(localFilter[..filterCount]);
-            }
-
-            if (maxPrimeConfigured)
-            {
-                if (!addedWithinLimit)
-                {
-                    Console.WriteLine("No filter candidates fall within the --max-prime limit.");
-                    return;
-                }
-
-                if (skippedByLimit > 0)
-                {
-                    Console.WriteLine($"Skipped {skippedByLimit} filter candidate(s) above --max-prime limit.");
                 }
             }
         }
 
-
-        if (_testMode)
-        {
-            _testTargetPrimeCount = testPrimeCandidateLimit;
-        }
-
         if (_primeCount >= 2)
         {
             _primeFoundAfterInit = true;
@@ -753,12 +661,6 @@ internal static class Program
         // Configure batch size for GPU primality sieve
         PrimeTester.GpuBatchSize = gpuPrimeBatch;
 
-        Stopwatch? stopwatch = null;
-        if (_testMode)
-        {
-            stopwatch = Stopwatch.StartNew();
-        }
-
         Console.WriteLine("Starting scan...");
 
         if (_useByDivisorMode)
@@ -766,12 +668,6 @@ internal static class Program
             RunByDivisorMode(byDivisorCandidates, threadCount);
             FlushBuffer();
             StringBuilderPool.Return(_outputBuilder!);
-            if (stopwatch is not null)
-            {
-                stopwatch.Stop();
-                ReportTestTime(stopwatch.Elapsed);
-            }
-
             return;
         }
 
@@ -786,14 +682,7 @@ internal static class Program
             }
         }
 
-        if (_testMode)
-        {
-            testPrimeCandidateLimit = Math.Max(1, threadCount) * 3;
-            _testTargetPrimeCount = testPrimeCandidateLimit;
-        }
-
         _state = ((long)currentP << 3) | (long)remainder;
-
         Task[] tasks = new Task[threadCount];
 
         int taskIndex = 0;
@@ -863,71 +752,13 @@ internal static class Program
         }
 
         Task.WaitAll(tasks);
-
-        TimeSpan testElapsed = default;
-        bool reportTestTime = false;
-        if (stopwatch is not null)
-        {
-            stopwatch.Stop();
-            testElapsed = stopwatch.Elapsed;
-            reportTestTime = true;
-        }
-
         FlushBuffer();
         StringBuilderPool.Return(_outputBuilder!);
-
-        if (reportTestTime)
-        {
-            ReportTestTime(testElapsed);
-        }
     }
 
-    private static List<ulong> BuildTestPrimeCandidates(int targetCount)
-    {
-        if (targetCount <= 0)
-        {
-            return [];
-        }
-
-        List<ulong> candidates = new(targetCount);
-        ulong previous = 29UL;
-        int remaining = targetCount;
-        ulong nextPrime = 0UL;
-
-        while (remaining > 0)
-        {
-            try
-            {
-                nextPrime = PrimeIterator.Next(in previous);
-            }
-            catch (InvalidOperationException)
-            {
-                break;
-            }
-
-            if (nextPrime < 31UL)
-            {
-                previous = nextPrime;
-                continue;
-            }
-
-            candidates.Add(nextPrime);
-            remaining--;
-            previous = nextPrime;
-        }
-
-        if (remaining > 0)
-        {
-            Console.WriteLine("Unable to populate the requested number of test primes before reaching the 64-bit limit.");
-        }
-
-        return candidates;
-    }
-
-    private static List<ulong> LoadByDivisorCandidates(string candidateFile, UInt128 maxPrimeLimit, bool maxPrimeConfigured, out int skippedByLimit)
+    private static List<ulong> LoadByDivisorCandidates(string candidateFile)
     {
         List<ulong> candidates = [];
-        skippedByLimit = 0;
         using FileStream readStream = new(candidateFile, FileMode.Open, FileAccess.Read, FileShare.ReadWrite);
         using StreamReader reader = new(readStream);
         string? line;
@@ -967,14 +798,7 @@ internal static class Program
                 // arrive so by-divisor candidate loading avoids per-value string allocations.
                 if (ulong.TryParse(span[start..index], NumberStyles.None, CultureInfo.InvariantCulture, out ulong parsed))
                 {
-                    if (!maxPrimeConfigured || (UInt128)parsed <= maxPrimeLimit)
-                    {
-                        candidates.Add(parsed);
-                    }
-                    else
-                    {
-                        skippedByLimit++;
-                    }
+                    candidates.Add(parsed);
                 }
             }
         }
@@ -1072,7 +896,6 @@ internal static class Program
         Console.WriteLine();
         Console.WriteLine("Options:");
         Console.WriteLine("  --prime=<value>        starting exponent (p)");
-        Console.WriteLine("  --max-prime=<value>    inclusive upper bound for primes from filter files");
         Console.WriteLine("  --increment=bit|add    exponent increment method");
         Console.WriteLine("  --threads=<value>      number of worker threads");
         Console.WriteLine("  --block-size=<value>   values processed per thread batch");
@@ -1166,21 +989,6 @@ internal static class Program
         }
     }
 
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static void OnPrimeCandidateConfirmed()
-    {
-        if (!_testMode)
-        {
-            return;
-        }
-
-        int processed = Interlocked.Increment(ref _testProcessedPrimeCount);
-        if (processed == _testTargetPrimeCount)
-        {
-            Volatile.Write(ref _limitReached, true);
-        }
-    }
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static bool TryAcquireConsoleSlot(bool primeCandidate, out bool primeFlag)
     {
@@ -1317,75 +1125,6 @@ internal static class Program
         }
     }
 
-    private static void ReportTestTime(TimeSpan elapsed)
-    {
-        double elapsedSeconds = elapsed.TotalSeconds;
-        double elapsedMilliseconds = elapsed.TotalMilliseconds;
-        Console.WriteLine($"Test elapsed time: {elapsedSeconds:F3} s");
-
-        string filePath = GetTestTimeFilePath();
-        bool hasPrevious = TryReadTestTime(filePath, out double previousMilliseconds);
-        string elapsedText = elapsedMilliseconds.ToString("F3", CultureInfo.InvariantCulture);
-
-        if (!hasPrevious)
-        {
-            WriteTestTime(filePath, elapsedMilliseconds);
-            Console.WriteLine($"FIRST TIME: {elapsedText} ms");
-            return;
-        }
-
-        string previousText = previousMilliseconds.ToString("F3", CultureInfo.InvariantCulture);
-        if (elapsedMilliseconds < previousMilliseconds)
-        {
-            WriteTestTime(filePath, elapsedMilliseconds);
-            Console.WriteLine($"BETTER TIME: {elapsedText} ms (previous {previousText} ms)");
-        }
-        else
-        {
-            Console.WriteLine($"WORSE TIME: {elapsedText} ms (best {previousText} ms)");
-        }
-    }
-
-    private static string GetTestTimeFilePath()
-    {
-        if (string.IsNullOrEmpty(_resultsDir))
-        {
-            return TestTimeFileName;
-        }
-
-        return Path.Combine(_resultsDir!, TestTimeFileName);
-    }
-
-    private static bool TryReadTestTime(string filePath, out double milliseconds)
-    {
-        if (!File.Exists(filePath))
-        {
-            milliseconds = 0.0;
-            return false;
-        }
-
-        string content = File.ReadAllText(filePath).Trim();
-        if (content.Length == 0)
-        {
-            milliseconds = 0.0;
-            return false;
-        }
-
-        return double.TryParse(content, NumberStyles.Float, CultureInfo.InvariantCulture, out milliseconds);
-    }
-
-    private static void WriteTestTime(string filePath, double milliseconds)
-    {
-        string? directory = Path.GetDirectoryName(filePath);
-        if (!string.IsNullOrEmpty(directory) && !Directory.Exists(directory))
-        {
-            Directory.CreateDirectory(directory);
-        }
-
-        string formatted = milliseconds.ToString("F6", CultureInfo.InvariantCulture);
-        File.WriteAllText(filePath, formatted);
-    }
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     internal static ulong CountOnes(ulong value)
     {
@@ -1616,8 +1355,6 @@ internal static class Program
                 _lastCompositeP = true;
                 return false;
             }
-
-            OnPrimeCandidateConfirmed();
         }
 
         searchedMersenne = true;

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

@@ -1,7 +1,10 @@
+using System;
+using System.Collections.Generic;
 using System.Reflection;
 using FluentAssertions;
 using ILGPU;
 using ILGPU.Runtime;
+using PerfectNumbers.Core;
 using PerfectNumbers.Core.Gpu;
 using Xunit;
 
@@ -196,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/Gpu/PrimeOrderGpuHeuristicsTests.cs

@@ -1,6 +1,8 @@
+using System;
 using System.Numerics;
 using System.Reflection;
 using FluentAssertions;
+using PerfectNumbers.Core;
 using PerfectNumbers.Core.Gpu;
 using Xunit;
 

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

@@ -1,4 +1,7 @@
+using System;
+using System.Collections.Generic;
 using FluentAssertions;
+using PerfectNumbers.Core;
 using Xunit;
 
 namespace PerfectNumbers.Core.Tests;
@@ -67,7 +70,7 @@ public class MersenneNumberDivisorByDivisorTesterTests
             return new DummySession();
         }
 
-        public void PrepareCandidates(in ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues)
+        public void PrepareCandidates(ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues)
         {
             for (int i = 0; i < primes.Length; i++)
             {
@@ -77,7 +80,7 @@ public class MersenneNumberDivisorByDivisorTesterTests
 
         private sealed class DummySession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
         {
-            public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan<ulong> primes, Span<byte> hits)
+            public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
             {
                 throw new NotSupportedException();
             }

PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs

@@ -1,3 +1,4 @@
+using System;
 using FluentAssertions;
 using PerfectNumbers.Core.Cpu;
 using Xunit;

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -1,5 +1,8 @@
+using System;
+using System.Collections.Generic;
 using System.Reflection;
 using FluentAssertions;
+using PerfectNumbers.Core;
 using PerfectNumbers.Core.Cpu;
 using PerfectNumbers.Core.Gpu;
 using Xunit;

PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs

@@ -1,3 +1,5 @@
+using System.Collections.Generic;
+using System.Numerics;
 using FluentAssertions;
 using Xunit;
 
@@ -156,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);
     }

PerfectNumbers.Core.Tests/ULongExtensionsMultiplyShiftTests.cs

@@ -1,5 +1,6 @@
 using System.Numerics;
 using FluentAssertions;
+using PerfectNumbers.Core;
 using Xunit;
 
 namespace PerfectNumbers.Core.Tests;

PerfectNumbers.Core/Cpu/FactorCacheLease.cs

@@ -1,30 +0,0 @@
-namespace PerfectNumbers.Core.Cpu;
-
-internal struct FactorCacheLease
-{
-    private Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? _cache;
-
-    public void EnsureInitialized(ref Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? cache)
-    {
-        Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? current = cache;
-        if (current is null)
-        {
-            current = ThreadStaticPools.RentMersenneFactorCacheDictionary();
-            cache = current;
-        }
-
-        _cache = current;
-    }
-
-    public void Dispose()
-    {
-        Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? cache = _cache;
-        if (cache is null)
-        {
-            return;
-        }
-
-        ThreadStaticPools.ReturnMersenneFactorCacheDictionary(cache);
-        _cache = null;
-    }
-}

PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs

@@ -1,82 +0,0 @@
-namespace PerfectNumbers.Core.Cpu;
-
-internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
-{
-    private bool _disposed;
-
-    public void Reset()
-    {
-        _disposed = false;
-    }
-
-    public void CheckDivisor(
-        ulong divisor,
-        in MontgomeryDivisorData divisorData,
-        ulong divisorCycle,
-        in ReadOnlySpan<ulong> primes,
-        Span<byte> hits)
-    {
-        int length = primes.Length;
-        // EvenPerfectBitScanner always provides at least one exponent on this path.
-        // Keep the guard commented out for tests that might reuse the session differently.
-        // if (length == 0)
-        // {
-        //     return;
-        // }
-
-        MontgomeryDivisorData cachedData = divisorData;
-        if (cachedData.Modulus != divisor)
-        {
-            cachedData = MontgomeryDivisorDataCache.Get(divisor);
-        }
-
-        if (divisorCycle == 0UL)
-        {
-            divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
-            if (divisorCycle == 0UL)
-            {
-                hits.Clear();
-                return;
-            }
-        }
-
-        // Keep these remainder steppers in place so future updates continue reusing the previously computed residues.
-        // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors.
-        var exponentStepper = new ExponentRemainderStepper(cachedData);
-        if (!exponentStepper.IsValidModulus)
-        {
-            hits.Clear();
-            return;
-        }
-
-        var cycleStepper = new CycleRemainderStepper(divisorCycle);
-
-        ulong remainder = cycleStepper.Initialize(primes[0]);
-        hits[0] = remainder == 0UL
-            ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0)
-            : (byte)0;
-
-        for (int i = 1; i < length; i++)
-        {
-            remainder = cycleStepper.ComputeNext(primes[i]);
-            if (remainder != 0UL)
-            {
-                hits[i] = 0;
-                continue;
-            }
-
-            hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0;
-        }
-    }
-
-    public void Dispose()
-    {
-        if (_disposed)
-        {
-            throw new ObjectDisposedException(nameof(MersenneCpuDivisorScanSession), "Divisor scan session disposed twice.");
-        }
-
-        _disposed = true;
-        ThreadStaticPools.ReturnMersenneCpuDivisorSession(this);
-    }
-}

PerfectNumbers.Core/CycleRemainderStepper.cs

@@ -1,8 +1,9 @@
-using System.Runtime.CompilerServices;
+using System;
+using System.Numerics;
 
 namespace PerfectNumbers.Core;
 
-public struct CycleRemainderStepper
+internal struct CycleRemainderStepper
 {
     private readonly ulong _cycleLength;
     private ulong _previousPrime;
@@ -22,16 +23,15 @@ public struct CycleRemainderStepper
         _hasState = false;
     }
 
-	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public void Reset()
-	{
-		// TODO: Inline this reset at the call sites so the hot loops reuse struct reinitialization
-		// measured fastest in MersenneDivisorCycleLengthGpuBenchmarks, avoiding the extra
-		// method call when scanners need to restart stepping.
-		_previousPrime = 0UL;
-		_currentRemainder = 0UL;
-		_hasState = false;
-	}
+    public void Reset()
+    {
+        // TODO: Inline this reset at the call sites so the hot loops reuse struct reinitialization
+        // measured fastest in MersenneDivisorCycleLengthGpuBenchmarks, avoiding the extra
+        // method call when scanners need to restart stepping.
+        _previousPrime = 0UL;
+        _currentRemainder = 0UL;
+        _hasState = false;
+    }
 
     public ulong Initialize(ulong prime)
     {

PerfectNumbers.Core/DeterministicRandom.cs

@@ -1,3 +1,6 @@
+using System;
+using System.Threading;
+
 namespace PerfectNumbers.Core;
 
 internal static class DeterministicRandom

PerfectNumbers.Core/DivisorCycleCache.cs

@@ -1,7 +1,10 @@
+using System;
 using System.Buffers;
 using System.Collections.Concurrent;
+using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+using System.Threading;
 using ILGPU;
 using ILGPU.Runtime;
 using PerfectNumbers.Core.Gpu;

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/FactorEntry128.cs

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

PerfectNumbers.Core/Gpu/GpuContextPool.cs

@@ -126,7 +126,7 @@ public static class GpuContextPool
                 }
         }
 
-    public struct GpuContextLease
+    public struct GpuContextLease : IDisposable
     {
         private readonly PooledContext _ctx;
                 private bool disposed;

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

@@ -677,7 +677,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne
         // GC.SuppressFinalize(this);
     }
 
-    public readonly struct ExecutionScope
+    public readonly struct ExecutionScope : IDisposable
     {
         private readonly object? _lock;
 

PerfectNumbers.Core/Gpu/GpuRational.cs

@@ -1,3 +1,5 @@
+using System;
+using System.Numerics;
 using System.Runtime.CompilerServices;
 
 namespace PerfectNumbers.Core.Gpu;

PerfectNumbers.Core/Gpu/GpuUInt128.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using ILGPU.Algorithms;

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -1,9 +1,14 @@
+using System;
 using System.Buffers;
 using System.Collections.Concurrent;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+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;
 
@@ -143,7 +148,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         return true;
     }
 
-    public void PrepareCandidates(in ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues)
+    public void PrepareCandidates(ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues)
     {
         if (allowedMaxValues.Length < primes.Length)
         {
@@ -546,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulong exponent = exponents[index];
-        hits[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -578,7 +583,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             _disposed = false;
         }
 
-        public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan<ulong> primes, Span<byte> hits)
+        public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
         {
             if (_disposed)
             {
@@ -765,7 +770,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         ulong exponent = exponents[index];
 
-        results[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: true);
+        results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
     }
 
     public ulong DivisorLimit

PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs

@@ -1,5 +1,8 @@
+using System;
 using System.Buffers;
+using System.Threading;
 using System.Runtime.InteropServices;
+using ILGPU;
 using ILGPU.Runtime;
 
 namespace PerfectNumbers.Core.Gpu;

PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs

@@ -1,5 +1,6 @@
 using System.Buffers;
 using System.Runtime.InteropServices;
+using System.Threading;
 using ILGPU.Runtime;
 
 namespace PerfectNumbers.Core.Gpu;

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -1,8 +1,10 @@
+using System;
 using System.Buffers;
 using System.Collections.Concurrent;
 using ILGPU;
 using ILGPU.Runtime;
 using System.Numerics;
+using PerfectNumbers.Core;
 using System.Runtime.InteropServices;
 using System.Runtime.CompilerServices;
 using ILGPU.Runtime.OpenCL;
@@ -809,7 +811,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             ulong reduced = phi / factor;
-            if (ULongExtensions.Pow2MontgomeryModWindowedKernel(reduced, divisor, keepMontgomery: false) == 1UL)
+            if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL)
             {
                 return false;
             }
@@ -1593,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
     {
-        return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
     }
 
     private static ulong CalculateByDoublingKernel(ulong prime)
@@ -1793,7 +1795,7 @@ internal static partial class PrimeOrderGpuHeuristics
             return 0UL;
         }
 
-        return ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisorData, keepMontgomery: false);
+        return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
@@ -1809,7 +1811,7 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong exponent = exponents[index];
-        remainders[index] = ULongExtensions.Pow2MontgomeryModWindowedKernel(exponent, divisor, keepMontgomery: false);
+        remainders[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)

PerfectNumbers.Core/Gpu/ReadOnlyGpuUInt128.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 

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

@@ -1,3 +1,5 @@
+using System;
+
 namespace PerfectNumbers.Core;
 
 public interface IMersenneNumberDivisorByDivisorTester
@@ -14,10 +16,10 @@ public interface IMersenneNumberDivisorByDivisorTester
 
     IDivisorScanSession CreateDivisorSession();
 
-    void PrepareCandidates(in ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues);
+    void PrepareCandidates(ReadOnlySpan<ulong> primes, Span<ulong> allowedMaxValues);
 
     public interface IDivisorScanSession : IDisposable
     {
-        void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan<ulong> primes, Span<byte> hits);
+        void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits);
     }
 }

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -1,8 +1,7 @@
 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 PerfectNumbers.Core.Gpu;
 using ILGPU;
@@ -15,6 +14,7 @@ public class MersenneDivisorCycles
     private List<(ulong divisor, ulong cycleLength)> _table = [];
     // Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown
     private ulong[]? _smallCycles;
+
     public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles();
 
     public MersenneDivisorCycles()
@@ -51,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;
         }
@@ -94,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;
         }
@@ -102,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))
         {
@@ -261,8 +261,8 @@ public class MersenneDivisorCycles
         return order;
     }
 
-    public static bool TryCalculateCycleLengthForExponent(
-        ulong divisor,
+    public static bool TryCalculateCycleLengthForExponentCpu(
+        ulong divisor,        
         ulong exponent,
         in MontgomeryDivisorData divisorData,
         Dictionary<ulong, FactorCacheEntry>? factorCache,
@@ -282,25 +282,23 @@ public class MersenneDivisorCycles
             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;
-        // A zero totient would imply divisor == 0, which never occurs on the tested path. Keep the guard documented without
-        // branching in the hot loop.
-        // if (phi == 0UL)
-        // {
-        //     return false;
-        // }
+        if (phi == 0UL)
+        {
+            return false;
+        }
 
         int twoCount = BitOperations.TrailingZeroCount(phi);
         ulong reducedPhi = phi >> twoCount;
@@ -310,32 +308,24 @@ public class MersenneDivisorCycles
         }
 
         ulong k = reducedPhi / exponent;
-        // k == 0 would indicate reducedPhi < exponent, which the EvenPerfectBitScanner call flow prevents. Document the
-        // assumption instead of branching.
-        // if (k == 0UL)
-        // {
-        //     return false;
-        // }
-
-        Dictionary<ulong, int> factorCounts = ThreadStaticPools.RentFactorCountDictionary();
-        bool success = false;
+        if (k == 0UL)
+        {
+            return false;
+        }
 
+        Dictionary<ulong, int> factorCounts = new(8);
         if (twoCount > 0)
         {
             factorCounts[2UL] = twoCount;
         }
 
-        bool exponentAccumulated = AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true);
-        bool factorsAccumulated = exponentAccumulated && AccumulateFactors(k, factorCache, factorCounts, cacheResult: false);
-
-        if (factorsAccumulated)
+        if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
         {
-            cycleLength = ReduceOrder(divisorData, phi, factorCounts);
-            success = true;
+            return false;
         }
 
-        ThreadStaticPools.ReturnFactorCountDictionary(factorCounts);
-        return success;
+        cycleLength = ReduceOrderCpu(divisorData, phi, factorCounts);
+        return true;
     }
 
     private static bool AccumulateFactors(
@@ -349,27 +339,18 @@ public class MersenneDivisorCycles
             return true;
         }
 
-        if (!TryGetFactorization(value, cache, cacheResult, out FactorCacheEntry? factorization, out bool releaseAfterUse) || factorization is null)
+        if (!TryGetFactorization(value, cache, cacheResult, out FactorCacheEntry factorization))
         {
             return false;
         }
 
-        FactorCacheEntry actualFactorization = factorization;
-        KeyValuePair<ulong, int>[]? pairs = actualFactorization.Pairs;
-        if (pairs is not null)
-        {
-            int length = actualFactorization.Count;
+        ulong[] primes = factorization.Primes;
+        byte[] exponents = factorization.Exponents;
+        int length = factorization.Count;
 
-            for (int i = 0; i < length; i++)
-            {
-                KeyValuePair<ulong, int> pair = pairs[i];
-                AddFactor(counts, pair.Key, pair.Value);
-            }
-        }
-
-        if (releaseAfterUse)
+        for (int i = 0; i < length; i++)
         {
-            actualFactorization.ReturnToPool();
+            AddFactor(counts, primes[i], exponents[i]);
         }
 
         return true;
@@ -379,57 +360,37 @@ public class MersenneDivisorCycles
         ulong value,
         Dictionary<ulong, FactorCacheEntry>? cache,
         bool cacheResult,
-        out FactorCacheEntry? factorization,
-        out bool releaseAfterUse)
+        out FactorCacheEntry factorization)
     {
-        releaseAfterUse = false;
-
-        if (cache is not null && cache.TryGetValue(value, out FactorCacheEntry? cachedEntry))
+        if (cache is not null && cache.TryGetValue(value, out factorization))
         {
-            factorization = cachedEntry;
             return true;
         }
 
-        Dictionary<ulong, int> scratch = ThreadStaticPools.RentFactorScratchDictionary();
-        bool factored = TryFactorIntoCountsInternal(value, scratch);
-        if (!factored)
+        Dictionary<ulong, int> scratch = new(8);
+        if (!TryFactorIntoCountsInternal(value, scratch))
         {
-            ThreadStaticPools.ReturnFactorScratchDictionary(scratch);
-            factorization = null;
+            factorization = default!;
             return false;
         }
 
         int count = scratch.Count;
-        FactorCacheEntry entry = ThreadStaticPools.RentFactorCacheEntry();
-        if (count == 0)
+        ulong[] primes = new ulong[count];
+        byte[] exponents = new byte[count];
+        int index = 0;
+        foreach (KeyValuePair<ulong, int> entry in scratch)
         {
-            entry.Initialize(null, 0, fromPool: false);
-            factorization = entry;
-            releaseAfterUse = cache is null || !cacheResult;
-            if (cache is not null && cacheResult)
-            {
-                cache[value] = entry;
-            }
-
-            ThreadStaticPools.ReturnFactorScratchDictionary(scratch);
-            return true;
+            primes[index] = entry.Key;
+            exponents[index] = checked((byte)entry.Value);
+            index++;
         }
 
-        KeyValuePair<ulong, int>[] pairs = ThreadStaticPools.FactorKeyValuePairPool.Rent(count);
-        ((ICollection<KeyValuePair<ulong, int>>)scratch).CopyTo(pairs, 0);
-        ThreadStaticPools.ReturnFactorScratchDictionary(scratch);
-        Array.Sort(pairs, 0, count, FactorPairValueComparer.Instance);
-
-        entry.Initialize(pairs, count, fromPool: true);
-        factorization = entry;
+        Array.Sort(primes, exponents);
+        factorization = new FactorCacheEntry(primes, exponents, count);
 
         if (cache is not null && cacheResult)
         {
-            cache[value] = entry;
-        }
-        else
-        {
-            releaseAfterUse = true;
+            cache[value] = factorization;
         }
 
         return true;
@@ -437,11 +398,10 @@ public class MersenneDivisorCycles
 
     private static bool TryFactorIntoCountsInternal(ulong value, Dictionary<ulong, int> counts)
     {
-        // The EvenPerfectBitScanner only requests factors for values >= 2 here. Leave the guard documented without executing it.
-        // if (value <= 1UL)
-        // {
-        //     return true;
-        // }
+        if (value <= 1UL)
+        {
+            return true;
+        }
 
         ulong remaining = value;
         uint[] smallPrimes = PrimesGenerator.SmallPrimes;
@@ -456,7 +416,6 @@ public class MersenneDivisorCycles
                 break;
             }
 
-            // TODO: Reuse the remainder-tracking update from the cycle heuristic so repeated divisions disappear from this loop.
             while ((remaining % prime) == 0UL)
             {
                 AddFactor(counts, prime);
@@ -469,43 +428,20 @@ public class MersenneDivisorCycles
             return true;
         }
 
-        // TODO: Benchmark PrimeTester.IsPrimeInternal against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites.
         if (PrimeTester.IsPrimeInternal(remaining, CancellationToken.None))
         {
             AddFactor(counts, remaining);
             return true;
         }
 
-        // TODO: Evaluate folding the remainder-tracking mechanism into PollardRho64 so the walk leverages previously computed residues.
         ulong factor = PollardRho64(remaining);
         if (factor == 0UL || factor == remaining)
         {
-            // Pollard-Rho can still stall on rare stubborn composites during test and benchmark runs, so keep the failure guard in place
-            // and let the caller fall back to the cached cycle calculator.
             return false;
         }
 
-        // TODO: Investigate computing this quotient via remainder tracking or reciprocal multiplication to avoid the division.
         ulong quotient = remaining / factor;
-        // Pollard-Rho never returns factors <= 1 on the EvenPerfectBitScanner path; keep the old guard documented
-        // for completeness without branching in the hot path.
-        // if (factor > 1UL && !TryFactorIntoCountsInternal(factor, counts))
-        // {
-        //     return false;
-        // }
-
-        if (!TryFactorIntoCountsInternal(factor, counts))
-        {
-            return false;
-        }
-
-        // Guard the recursive calls to terminate once Pollard-Rho reduces either branch to 1.
-        if (quotient > 1UL && !TryFactorIntoCountsInternal(quotient, counts))
-        {
-            return false;
-        }
-
-        return true;
+        return TryFactorIntoCountsInternal(factor, counts) && TryFactorIntoCountsInternal(quotient, counts);
     }
 
     private static void AddFactor(Dictionary<ulong, int> counts, ulong factor)
@@ -532,7 +468,7 @@ public class MersenneDivisorCycles
         }
     }
 
-    private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
+    private static ulong ReduceOrderCpu(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
     {
         if (factorCounts.Count == 0)
         {
@@ -591,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);
             }
 
@@ -615,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)
@@ -653,73 +593,19 @@ public class MersenneDivisorCycles
         }
     }
 
-    private static ulong NextRandomUInt64()
+    public readonly struct FactorCacheEntry
     {
-        Span<byte> buffer = stackalloc byte[8];
-        RandomNumberGenerator.Fill(buffer);
-        return BinaryPrimitives.ReadUInt64LittleEndian(buffer);
-    }
-
-    private sealed class FactorPairValueComparer : IComparer<KeyValuePair<ulong, int>>
-    {
-        public static readonly FactorPairValueComparer Instance = new FactorPairValueComparer();
-
-        private FactorPairValueComparer()
+        public FactorCacheEntry(ulong[] primes, byte[] exponents, int count)
         {
-        }
-
-        public int Compare(KeyValuePair<ulong, int> x, KeyValuePair<ulong, int> y)
-        {
-            if (x.Key < y.Key)
-            {
-                return -1;
-            }
-
-            if (x.Key > y.Key)
-            {
-                return 1;
-            }
-
-            return 0;
-        }
-    }
-
-    public sealed class FactorCacheEntry
-    {
-        internal FactorCacheEntry? Next;
-        internal KeyValuePair<ulong, int>[]? Pairs;
-        internal int Count;
-        private bool _fromPool;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal void Reset()
-        {
-            Pairs = null;
-            Count = 0;
-            _fromPool = false;
-            Next = null;
-        }
-
-        internal void Initialize(KeyValuePair<ulong, int>[]? pairs, int count, bool fromPool)
-        {
-            Pairs = pairs;
+            Primes = primes;
+            Exponents = exponents;
             Count = count;
-            _fromPool = fromPool;
-            Next = null;
         }
+        public ulong[] Primes { get; }
 
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal void ReturnToPool()
-        {
-            KeyValuePair<ulong, int>[]? pairs = Pairs;
-            if (_fromPool && pairs is not null)
-            {
-                ThreadStaticPools.FactorKeyValuePairPool.Return(pairs, clearArray: false);
-            }
+        public byte[] Exponents { get; }
 
-            Reset();
-            ThreadStaticPools.ReturnFactorCacheEntry(this);
-        }
+        public int Count { get; }
     }
 
     public static (long nextPosition, long completeCount) FindLast(string path)
@@ -1048,15 +934,15 @@ public class MersenneDivisorCycles
 
         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

@@ -1,262 +1,319 @@
+using System;
 using System.Buffers;
+using System.Collections.Generic;
 using System.Globalization;
+using System.Threading.Tasks;
+using Open.Numeric.Primes;
 using System.Runtime.InteropServices;
 
 namespace PerfectNumbers.Core;
 
 public static class MersenneNumberDivisorByDivisorTester
 {
-	public static void Run(
-			List<ulong> candidates,
-			IMersenneNumberDivisorByDivisorTester tester,
-			Dictionary<ulong, (bool DetailedCheck, bool PassedAllTests)>? previousResults,
-			ulong startPrime,
-			Action markComposite,
-			Action clearComposite,
-			Action<ulong, bool, bool, bool> printResult,
-			int threadCount)
-	{
-		if (candidates.Count == 0)
-		{
-			Console.WriteLine("No candidates were provided for --mersenne=bydivisor.");
-			return;
-		}
+    public static void Run(
+            List<ulong> candidates,
+            IMersenneNumberDivisorByDivisorTester tester,
+            Dictionary<ulong, (bool DetailedCheck, bool PassedAllTests)>? previousResults,
+            ulong startPrime,
+            Action markComposite,
+            Action clearComposite,
+            Action<ulong, bool, bool, bool> printResult,
+            int threadCount)
+    {
+        if (candidates.Count == 0)
+        {
+            Console.WriteLine("No candidates were provided for --mersenne=bydivisor.");
+            return;
+        }
 
-		candidates.Sort();
-		Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
-		int candidateCount = candidateSpan.Length;
+        candidates.Sort();
+        Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
+        int candidateCount = candidateSpan.Length;
 
-		bool applyStartPrime = startPrime > 0UL;
-		int skippedByPreviousResults = 0;
+        bool applyStartPrime = startPrime > 0UL;
+        int skippedByPreviousResults = 0;
 
-		ArrayPool<ulong> pool = ThreadStaticPools.UlongPool;
-		if (previousResults is not null && previousResults.Count > 0)
-		{
-			int recordedCount = previousResults.Count;
-			ulong[] recordedCandidates = pool.Rent(recordedCount);
+        if (previousResults is not null && previousResults.Count > 0)
+        {
+            int recordedCount = previousResults.Count;
+            bool preferSortedMerge = recordedCount <= 4096 || (recordedCount * 4L) <= candidateCount;
+            if (preferSortedMerge)
+            {
+                ulong[] recordedCandidates = ArrayPool<ulong>.Shared.Rent(recordedCount);
+                try
+                {
+                    Span<ulong> recordedSpan = recordedCandidates.AsSpan(0, recordedCount);
+                    int recordedIndex = 0;
+                    foreach (ulong value in previousResults.Keys)
+                    {
+                        recordedSpan[recordedIndex++] = value;
+                    }
 
-			Span<ulong> recordedSpan = recordedCandidates.AsSpan(0, recordedCount);
-			int recordedIndex = 0;
-			foreach (ulong value in previousResults.Keys)
-			{
-				recordedSpan[recordedIndex++] = value;
-			}
+                    recordedSpan.Sort();
 
-			recordedSpan.Sort();
+                    int writeIndex = 0;
+                    int skipIndex = 0;
 
-			int writeIndex = 0;
-			int skipIndex = 0;
+                    for (int readIndex = 0; readIndex < candidateCount; readIndex++)
+                    {
+                        ulong candidate = candidateSpan[readIndex];
 
-			for (int readIndex = 0; readIndex < candidateCount; readIndex++)
-			{
-				ulong candidate = candidateSpan[readIndex];
+                        while (skipIndex < recordedCount && recordedSpan[skipIndex] < candidate)
+                        {
+                            skipIndex++;
+                        }
 
-				while (skipIndex < recordedCount && recordedSpan[skipIndex] < candidate)
-				{
-					skipIndex++;
-				}
+                        if (skipIndex < recordedCount && recordedSpan[skipIndex] == candidate)
+                        {
+                            skippedByPreviousResults++;
+                            skipIndex++;
+                            continue;
+                        }
 
-				if (skipIndex < recordedCount && recordedSpan[skipIndex] == candidate)
-				{
-					skippedByPreviousResults++;
-					skipIndex++;
-					continue;
-				}
+                        candidateSpan[writeIndex++] = candidate;
+                    }
 
-				candidateSpan[writeIndex++] = candidate;
-			}
+                    if (writeIndex < candidateCount)
+                    {
+                        candidates.RemoveRange(writeIndex, candidateCount - writeIndex);
+                        candidateSpan = CollectionsMarshal.AsSpan(candidates); // candidateSpan now tracks the reduced list.
+                        candidateCount = candidateSpan.Length;
+                    }
+                }
+                finally
+                {
+                    ArrayPool<ulong>.Shared.Return(recordedCandidates, clearArray: true);
+                }
+            }
+            else
+            {
+                int writeIndex = 0;
+                for (int readIndex = 0; readIndex < candidateCount; readIndex++)
+                {
+                    ulong candidate = candidateSpan[readIndex];
 
-			if (writeIndex < candidateCount)
-			{
-				candidates.RemoveRange(writeIndex, candidateCount - writeIndex);
-				candidateSpan = CollectionsMarshal.AsSpan(candidates);
-				candidateCount = candidateSpan.Length;
-			}
+                    if (previousResults.ContainsKey(candidate))
+                    {
+                        skippedByPreviousResults++;
+                        continue;
+                    }
 
-			pool.Return(recordedCandidates, clearArray: true);
-		}
+                    candidateSpan[writeIndex++] = candidate;
+                }
 
-		if (skippedByPreviousResults > 0)
-		{
-			Console.WriteLine($"Skipped {skippedByPreviousResults.ToString(CultureInfo.InvariantCulture)} candidates excluded by previous results.");
-		}
+                if (writeIndex < candidateCount)
+                {
+                    candidates.RemoveRange(writeIndex, candidateCount - writeIndex);
+                    candidateSpan = CollectionsMarshal.AsSpan(candidates); // candidateSpan now tracks the reduced list.
+                    candidateCount = candidateSpan.Length;
+                }
+            }
+        }
 
-		int startIndex = 0;
-		if (applyStartPrime)
-		{
-			startIndex = candidates.BinarySearch(startPrime);
-			if (startIndex < 0)
-			{
-				startIndex = ~startIndex;
-			}
-			else
-			{
-				while (startIndex > 0 && candidates[startIndex - 1] >= startPrime)
-				{
-					startIndex--;
-				}
-			}
-		}
+        if (skippedByPreviousResults > 0)
+        {
+            Console.WriteLine($"Skipped {skippedByPreviousResults.ToString(CultureInfo.InvariantCulture)} candidates excluded by previous results.");
+        }
 
-		List<ulong> primesToTest = candidates;
-		Span<ulong> primesSpan = CollectionsMarshal.AsSpan(primesToTest);
-		ulong maxPrime = 0UL;
-		int primeWriteIndex = 0;
+        int startIndex = 0;
+        if (applyStartPrime)
+        {
+            startIndex = candidates.BinarySearch(startPrime);
+            if (startIndex < 0)
+            {
+                startIndex = ~startIndex;
+            }
+            else
+            {
+                while (startIndex > 0 && candidates[startIndex - 1] >= startPrime)
+                {
+                    startIndex--;
+                }
+            }
+        }
 
-		// 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.
+        List<ulong> primesToTest = candidates;
+        Span<ulong> primesSpan = CollectionsMarshal.AsSpan(primesToTest);
+        ulong maxPrime = 0UL;
+        int primeWriteIndex = 0;
 
-		// 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++)
-		{
-			ulong candidate = primesSpan[index];
+        for (int index = startIndex; index < candidateCount; index++)
+        {
+            ulong candidate = primesSpan[index];
 
-			if (candidate <= 1UL)
-			{
-				markComposite();
-				printResult(candidate, false, false, false);
-				continue;
-			}
+            if (candidate <= 1UL)
+            {
+                markComposite();
+                printResult(candidate, false, false, false);
+                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)
+            {
+                markComposite();
+                printResult(candidate, false, false, false);
+                continue;
+            }
 
-			// if (!hasPrime || currentPrime != candidate)
-			// if (!Number.IsPrime(candidate))
-			// {
-			//     markComposite();
-			//     printResult(candidate, false, false, false);
-			//     continue;
-			// }
+            primesSpan[primeWriteIndex++] = candidate;
 
-			primesSpan[primeWriteIndex++] = candidate;
+            if (candidate > maxPrime)
+            {
+                maxPrime = candidate;
+            }
+        }
 
-			if (candidate > maxPrime)
-			{
-				maxPrime = candidate;
-			}
-		}
+        if (primeWriteIndex < primesToTest.Count)
+        {
+            primesToTest.RemoveRange(primeWriteIndex, primesToTest.Count - primeWriteIndex);
+        }
 
-		if (primeWriteIndex < primesToTest.Count)
-		{
-			primesToTest.RemoveRange(primeWriteIndex, primesToTest.Count - primeWriteIndex);
-		}
+        if (primesToTest.Count == 0)
+        {
+            if (applyStartPrime)
+            {
+                Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor.");
+            }
+            else
+            {
+                Console.WriteLine("No prime candidates remain for --mersenne=bydivisor after filtering.");
+            }
 
-		if (primesToTest.Count == 0)
-		{
-			if (applyStartPrime)
-			{
-				Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor.");
-			}
-			else
-			{
-				Console.WriteLine("No prime candidates remain for --mersenne=bydivisor after filtering.");
-			}
+            return;
+        }
 
-			return;
-		}
+        if (maxPrime <= 1UL)
+        {
+            Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor.");
+            return;
+        }
 
-		if (maxPrime <= 1UL)
-		{
-			Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor.");
-			return;
-		}
+        tester.ConfigureFromMaxPrime(maxPrime);
 
-		tester.ConfigureFromMaxPrime(maxPrime);
+        int primeCount = primesToTest.Count;
+        ulong[] primeBatch = ArrayPool<ulong>.Shared.Rent(primeCount);
+        ulong[] allowedMaxBatch = ArrayPool<ulong>.Shared.Rent(primeCount);
+        List<ulong> filteredPrimes = new(primeCount);
 
-		int primeCount = primesToTest.Count;
-		ulong[] primeBatch = pool.Rent(primeCount);
-		ulong[] allowedMaxBatch = pool.Rent(primeCount);
-		List<ulong> filteredPrimes = new(primeCount);
+        try
+        {
+            Span<ulong> primeSpan = primeBatch.AsSpan(0, primeCount);
+            Span<ulong> allowedMaxSpan = allowedMaxBatch.AsSpan(0, primeCount);
 
-		Span<ulong> primeSpan = primeBatch.AsSpan(0, primeCount);
-		Span<ulong> allowedMaxSpan = allowedMaxBatch.AsSpan(0, primeCount);
+            for (int i = 0; i < primeCount; i++)
+            {
+                primeSpan[i] = primesToTest[i];
+            }
 
-		for (int i = 0; i < primeCount; i++)
-		{
-			primeSpan[i] = primesToTest[i];
-		}
+            tester.PrepareCandidates(primeSpan, allowedMaxSpan);
 
-		tester.PrepareCandidates(primeSpan, allowedMaxSpan);
+            for (int i = 0; i < primeCount; i++)
+            {
+                ulong prime = primeSpan[i];
+                ulong allowedMax = allowedMaxSpan[i];
 
-		for (int i = 0; i < primeCount; i++)
-		{
-			ulong prime = primeSpan[i];
-			ulong allowedMax = allowedMaxSpan[i];
+                if (allowedMax < 3UL)
+                {
+                    clearComposite();
+                    printResult(prime, true, true, true);
+                    continue;
+                }
 
-			if (allowedMax < 3UL)
-			{
-				clearComposite();
-				printResult(prime, true, true, true);
-				continue;
-			}
+                filteredPrimes.Add(prime);
+            }
+        }
+        finally
+        {
+            ArrayPool<ulong>.Shared.Return(allowedMaxBatch, clearArray: true);
+            ArrayPool<ulong>.Shared.Return(primeBatch, clearArray: true);
+        }
 
-			filteredPrimes.Add(prime);
-		}
+        if (filteredPrimes.Count == 0)
+        {
+            if (applyStartPrime)
+            {
+                Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor.");
+            }
 
-		pool.Return(allowedMaxBatch, clearArray: true);
-		pool.Return(primeBatch, clearArray: true);
+            return;
+        }
 
-		if (filteredPrimes.Count == 0)
-		{
-			if (applyStartPrime)
-			{
-				Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor.");
-			}
+        int workerCount = threadCount <= 0 ? Environment.ProcessorCount : threadCount;
+        if (workerCount < 1)
+        {
+            workerCount = 1;
+        }
 
-			return;
-		}
+        void ProcessPrime(ulong prime)
+        {
+            bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted);
 
-		int workerCount = threadCount <= 0 ? Environment.ProcessorCount : threadCount;
-		if (workerCount < 1)
-		{
-			workerCount = 1;
-		}
+            if (!isPrime)
+            {
+                markComposite();
+                printResult(prime, true, true, false);
+                return;
+            }
 
-		void ProcessPrime(ulong prime)
-		{
-			bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted);
+            clearComposite();
+            printResult(prime, true, divisorsExhausted, true);
+        }
 
-			if (!isPrime)
-			{
-				markComposite();
-				printResult(prime, true, true, false);
-				return;
-			}
+        if (workerCount == 1)
+        {
+            foreach (ulong prime in filteredPrimes)
+            {
+                ProcessPrime(prime);
+            }
+        }
+        else
+        {
+            ParallelOptions options = new()
+            {
+                MaxDegreeOfParallelism = workerCount
+            };
 
-			clearComposite();
-			printResult(prime, true, divisorsExhausted, true);
-		}
+            primeCount = filteredPrimes.Count; // Reusing primeCount to track the filtered count.
+            int primesPerWorker = (primeCount + workerCount - 1) / workerCount;
 
-		if (workerCount == 1)
-		{
-			foreach (ulong prime in filteredPrimes)
-			{
-				ProcessPrime(prime);
-			}
-		}
-		else
-		{
-			ParallelOptions options = new()
-			{
-				MaxDegreeOfParallelism = workerCount
-			};
+            Parallel.For(
+                    0,
+                    workerCount,
+                    options,
+                    workerIndex =>
+                    {
+                        int start = workerIndex * primesPerWorker;
+                        if (start >= primeCount)
+                        {
+                            return;
+                        }
 
-			Parallel.ForEach(filteredPrimes, options, ProcessPrime);
-		}
-	}
+                        int end = start + primesPerWorker;
+                        if (end > primeCount)
+                        {
+                            end = primeCount;
+                        }
+
+                        for (int i = start; i < end; i++)
+                        {
+                            ProcessPrime(filteredPrimes[i]);
+                        }
+                    });
+        }
+    }
 }

PerfectNumbers.Core/MontgomeryDivisorData.cs

@@ -1,3 +1,5 @@
+using System;
+using System.Numerics;
 using System.Runtime.CompilerServices;
 
 namespace PerfectNumbers.Core;

PerfectNumbers.Core/PartialFactorResult.cs

@@ -17,31 +17,27 @@ namespace PerfectNumbers.Core
         {
         }
 
-        public static PartialFactorResult Rent(in FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count)
-        {
-            PartialFactorResult? instance = s_poolHead;
-            if (instance is null)
-            {
-                instance = new PartialFactorResult();
-            }
-            else
-            {
-                s_poolHead = instance._next;
-            }
-
-            instance._next = null;
-            instance.Factors = factors;
-            instance.Cofactor = cofactor;
-            instance.FullyFactored = fullyFactored;
-            instance.Count = count;
-            return instance;
-        }
-
-        public static readonly PartialFactorResult Empty = new()
+		public static PartialFactorResult Rent(in FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count)
 		{
-            Cofactor = 1UL,
-            FullyFactored = true,
-        };
+			PartialFactorResult? instance = s_poolHead;
+			if (instance is null)
+			{
+				instance = new();
+			}
+			else
+			{
+				s_poolHead = instance._next;
+			}
+
+			instance._next = null;
+			instance.Factors = factors;
+			instance.Cofactor = cofactor;
+			instance.FullyFactored = fullyFactored;
+			instance.Count = count;
+			return instance;
+		}
+
+		public static readonly PartialFactorResult Empty = new() { Cofactor = 1UL, FullyFactored = true };
 
         public PartialFactorResult WithAdditionalPrime(ulong prime)
         {
@@ -62,19 +58,19 @@ namespace PerfectNumbers.Core
         }
 
         public void Dispose()
-        {
-            if (this != Empty)
-            {
-                FactorEntry[]? factors = Factors;
-                if (factors is not null)
-                {
-                    Factors = null;
-                    ThreadStaticPools.FactorEntryPool.Return(factors, clearArray: false);
-                }
+		{
+			if (this != Empty)
+			{
+				FactorEntry[]? factors = Factors;
+				if (factors is not null)
+				{
+					Factors = null;
+					ArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
+				}
 
-                _next = s_poolHead;
-                s_poolHead = this;
-            }
+				_next = s_poolHead;
+				s_poolHead = this;
+			}
         }
     }
-}
+}

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/PrimeOrderCalculator.Gpu.cs

@@ -1,45 +0,0 @@
-using PerfectNumbers.Core.Gpu;
-
-namespace PerfectNumbers.Core;
-
-internal static partial class PrimeOrderCalculator
-{
-        private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
-
-        private const int GpuSmallPrimeFactorSlots = 64;
-
-        private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary<ulong, int> counts, out ulong remaining)
-        {
-                ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots);
-                int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots);
-                Span<ulong> primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
-                Span<int> exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
-                primeBuffer.Clear();
-                exponentBuffer.Clear();
-
-                if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _))
-                {
-                        remaining = value;
-                        ThreadStaticPools.UlongPool.Return(primeBufferArray);
-                        ThreadStaticPools.IntPool.Return(exponentBufferArray);
-                        return false;
-                }
-
-                remaining = gpuRemaining;
-                for (int i = 0; i < factorCount; i++)
-                {
-                        ulong primeValue = primeBuffer[i];
-                        int exponent = exponentBuffer[i];
-                        if (primeValue == 0UL || exponent == 0)
-                        {
-                                continue;
-                        }
-
-                        counts[primeValue] = exponent;
-                }
-
-                ThreadStaticPools.UlongPool.Return(primeBufferArray);
-                ThreadStaticPools.IntPool.Return(exponentBufferArray);
-                return true;
-        }
-}

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -1,88 +1,58 @@
 using System.Buffers;
+using System.Collections.Generic;
 using System.Diagnostics;
 using System.Numerics;
+using System.Threading;
 using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
-        [ThreadStatic]
-        private static bool s_pow2ModeInitialized;
+    public static UInt128 Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
+        => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
 
-        [ThreadStatic]
-        private static bool s_allowGpuPow2;
-
-        [ThreadStatic]
-        private static PrimeOrderHeuristicDevice s_deviceMode;
-
-        [ThreadStatic]
-        private static bool s_debugLoggingEnabled;
-
-        private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu;
-
-        private readonly struct Pow2ModeScope
+    public static UInt128 Calculate(
+        UInt128 prime,
+        UInt128? previousOrder,
+        PrimeOrderSearchConfig config,
+        PrimeOrderHeuristicDevice device)
+    {
+        var scope = UsePow2Mode(device);
+        MontgomeryDivisorData divisorData;
+        UInt128 result;
+        if (prime <= ulong.MaxValue)
         {
-                private readonly bool _previousInitialized;
-                private readonly bool _previousAllow;
-                private readonly PrimeOrderHeuristicDevice _previousDevice;
-
-                public Pow2ModeScope(PrimeOrderHeuristicDevice device)
+            ulong? previous = null;
+            if (previousOrder.HasValue)
+            {
+                UInt128 previousValue = previousOrder.Value;
+                if (previousValue <= ulong.MaxValue)
                 {
-                        _previousInitialized = s_pow2ModeInitialized;
-                        _previousAllow = s_allowGpuPow2;
-                        _previousDevice = s_deviceMode;
-                        s_pow2ModeInitialized = true;
-                        s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu;
-                        s_deviceMode = device;
+                    previous = (ulong)previousValue;
                 }
-
-                public void Dispose()
+                else
                 {
-                        s_allowGpuPow2 = _previousAllow;
-                        s_pow2ModeInitialized = _previousInitialized;
-                        s_deviceMode = _previousDevice;
+                    previous = ulong.MaxValue;
                 }
+            }
+
+            ulong prime64 = (ulong)prime;
+            divisorData = MontgomeryDivisorData.FromModulus(prime64);
+            ulong order64 = Calculate(prime64, previous, divisorData, config, device);
+            result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
+        }
+        else
+        {
+            divisorData = default;
+            result = CalculateWideInternal(prime, previousOrder, divisorData, config);
         }
 
-        private readonly struct DebugLoggingScope
-        {
-                private readonly bool _previous;
+        scope.Dispose();
+        return result;
+    }
 
-                public DebugLoggingScope(bool enabled)
-                {
-                        _previous = s_debugLoggingEnabled;
-                        s_debugLoggingEnabled = enabled;
-                }
-
-                public void Dispose()
-                {
-                        s_debugLoggingEnabled = _previous;
-                }
-        }
-
-        private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device);
-
-        private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled);
-
-        [Conditional("DEBUG")]
-        private static void DebugLog(string message)
-        {
-                if (s_debugLoggingEnabled)
-                {
-                        Console.WriteLine(message);
-                }
-        }
-
-        [Conditional("DEBUG")]
-        private static void DebugLog(Func<string> messageFactory)
-        {
-                if (s_debugLoggingEnabled)
-                {
-                        Console.WriteLine(messageFactory());
-                }
-        }
-    private static UInt128 CalculateWideInternal(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
+    private static UInt128 CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
     {
         if (prime <= UInt128.One)
         {
@@ -112,12 +82,12 @@ internal static partial class PrimeOrderCalculator
     }
 
     private static UInt128 RunHeuristicPipelineWide(
-        in UInt128 prime,
-        in UInt128? previousOrder,
+        UInt128 prime,
+        UInt128? previousOrder,
         in MontgomeryDivisorData divisorData,
-        in PrimeOrderSearchConfig config,
-        in UInt128 phi,
-        in PartialFactorResult128 phiFactors)
+        PrimeOrderSearchConfig config,
+        UInt128 phi,
+        PartialFactorResult128 phiFactors)
     {
         if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
         {
@@ -145,12 +115,12 @@ internal static partial class PrimeOrderCalculator
         return candidateOrder;
     }
 
-    private static UInt128 FinishStrictlyWide(in UInt128 prime, in MontgomeryDivisorData divisorData)
+    private static UInt128 FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData)
     {
         return CalculateByFactorizationWide(prime, divisorData);
     }
 
-    private static bool TrySpecialMaxWide(in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
+    private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
     {
         ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
         int length = factors.Count;
@@ -167,7 +137,7 @@ internal static partial class PrimeOrderCalculator
         return true;
     }
 
-    private static UInt128 InitializeStartingOrderWide(in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData divisorData)
+    private static UInt128 InitializeStartingOrderWide(UInt128 prime, UInt128 phi, in MontgomeryDivisorData divisorData)
     {
         UInt128 order = phi;
         UInt128 mod8 = prime & (UInt128)7UL;
@@ -183,20 +153,19 @@ internal static partial class PrimeOrderCalculator
         return order;
     }
 
-    private static UInt128 ExponentLoweringWide(UInt128 order, in UInt128 prime, in MontgomeryDivisorData divisorData, in PartialFactorResult128 factors)
+    private static UInt128 ExponentLoweringWide(UInt128 order, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
     {
-		ArrayPool<FactorEntry128> pool = ThreadStaticPools.FactorEntry128Pool;
-		ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
-		int length = factors.Count;
-		// TODO: This should never trigger from production code - check
-		// if (length == 0)
-		// {
-		//     return order;
-		// }
-		FactorEntry128[]? tempArray = pool.Rent(length + 1);
+        FactorEntry128[]? tempArray = null;
         try
         {
+            ReadOnlySpan<FactorEntry128> factorSpan = factors.Factors;
+            int length = factors.Count;
+            if (length == 0)
+            {
+                return order;
+            }
 
+            tempArray = ArrayPool<FactorEntry128>.Shared.Rent(length + 1);
             Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
             factorSpan.CopyTo(buffer);
 
@@ -232,11 +201,14 @@ internal static partial class PrimeOrderCalculator
         }
         finally
         {
-			pool.Return(tempArray, clearArray: false);
+            if (tempArray is not null)
+            {
+                ArrayPool<FactorEntry128>.Shared.Return(tempArray, clearArray: false);
+            }
         }
     }
 
-    private static bool TryConfirmOrderWide(in UInt128 prime, in UInt128 order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
+    private static bool TryConfirmOrderWide(UInt128 prime, UInt128 order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
     {
         if (order == UInt128.Zero)
         {
@@ -294,11 +266,11 @@ internal static partial class PrimeOrderCalculator
     }
 
     private static bool TryHeuristicFinishWide(
-        in UInt128 prime,
-        in UInt128 order,
-        in UInt128? previousOrder,
+        UInt128 prime,
+        UInt128 order,
+        UInt128? previousOrder,
         in MontgomeryDivisorData divisorData,
-        in PrimeOrderSearchConfig config,
+        PrimeOrderSearchConfig config,
         out UInt128 result)
     {
         result = UInt128.Zero;
@@ -370,7 +342,7 @@ internal static partial class PrimeOrderCalculator
         return false;
     }
 
-    private static void SortCandidatesWide(UInt128 prime, in UInt128? previousOrder, List<UInt128> candidates)
+    private static void SortCandidatesWide(UInt128 prime, UInt128? previousOrder, List<UInt128> candidates)
     {
         UInt128 previous = previousOrder ?? UInt128.Zero;
         int previousGroup = previousOrder.HasValue ? GetGroupWide(previousOrder.Value, prime) : 1;
@@ -398,9 +370,9 @@ internal static partial class PrimeOrderCalculator
     }
 
     private static CandidateKey128 BuildKeyWide(
-        in UInt128 value,
-        in UInt128 prime,
-        in UInt128 previous,
+        UInt128 value,
+        UInt128 prime,
+        UInt128 previous,
         int previousGroup,
         bool hasPrevious)
     {
@@ -422,7 +394,7 @@ internal static partial class PrimeOrderCalculator
         return new CandidateKey128(primary, false, distance, false, value);
     }
 
-    private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY)
+    private static int CompareComponents(bool descendingX, UInt128 valueX, bool descendingY, UInt128 valueY)
     {
         if (descendingX == descendingY)
         {
@@ -447,7 +419,7 @@ internal static partial class PrimeOrderCalculator
         return valueX > valueY ? 1 : -1;
     }
 
-    private static int GetGroupWide(in UInt128 value, in UInt128 prime)
+    private static int GetGroupWide(UInt128 value, UInt128 prime)
     {
         UInt128 threshold1 = prime >> 3;
         if (value <= threshold1)
@@ -470,7 +442,7 @@ internal static partial class PrimeOrderCalculator
         return 0;
     }
 
-    private static void BuildCandidatesWide(in UInt128 order, FactorEntry128[] factors, int count, List<UInt128> candidates, int limit)
+    private static void BuildCandidatesWide(UInt128 order, FactorEntry128[] factors, int count, List<UInt128> candidates, int limit)
     {
         if (count == 0)
         {
@@ -483,10 +455,10 @@ internal static partial class PrimeOrderCalculator
     }
 
     private static void BuildCandidatesRecursiveWide(
-        in UInt128 order,
-        in ReadOnlySpan<FactorEntry128> factors,
+        UInt128 order,
+        ReadOnlySpan<FactorEntry128> factors,
         int index,
-        in UInt128 divisorProduct,
+        UInt128 divisorProduct,
         List<UInt128> candidates,
         int limit)
     {
@@ -542,7 +514,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool TryConfirmCandidateWide(in UInt128 prime, in UInt128 candidate, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData)
+    private static bool TryConfirmCandidateWide(UInt128 prime, UInt128 candidate, PrimeOrderSearchConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData divisorData)
     {
         PartialFactorResult128 factorization = PartialFactorWide(candidate, config);
         if (factorization.Factors is null)
@@ -595,7 +567,7 @@ internal static partial class PrimeOrderCalculator
         return true;
     }
 
-    private static PartialFactorResult128 PartialFactorWide(in UInt128 value, in PrimeOrderSearchConfig config)
+    private static PartialFactorResult128 PartialFactorWide(UInt128 value, PrimeOrderSearchConfig config)
     {
         if (value <= UInt128.One)
         {
@@ -636,12 +608,12 @@ internal static partial class PrimeOrderCalculator
 
         if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
         {
-            long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds);
+            Stopwatch stopwatch = Stopwatch.StartNew();
+            long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks;
             Stack<UInt128> stack = new();
             stack.Push(remaining);
             pending.Clear();
 
-            long timestamp = 0L; // reused for deadline checks.
             while (stack.Count > 0)
             {
                 UInt128 composite = stack.Pop();
@@ -656,14 +628,13 @@ internal static partial class PrimeOrderCalculator
                     continue;
                 }
 
-                timestamp = Stopwatch.GetTimestamp();
-                if (timestamp > deadlineTimestamp)
+                if (stopwatch.ElapsedTicks > budgetTicks)
                 {
                     pending.Add(composite);
                     continue;
                 }
 
-                if (!TryPollardRhoWide(composite, deadlineTimestamp, out UInt128 factor))
+                if (!TryPollardRhoWide(composite, stopwatch, budgetTicks, out UInt128 factor))
                 {
                     pending.Add(composite);
                     continue;
@@ -696,8 +667,7 @@ internal static partial class PrimeOrderCalculator
             return new PartialFactorResult128(null, value, false, 0);
         }
 
-		ArrayPool<FactorEntry128>? pool = null;
-		FactorEntry128[]? rented = counts.Count > 0 ?  (pool = ThreadStaticPools.FactorEntry128Pool).Rent(counts.Count) : null;
+        FactorEntry128[]? rented = counts.Count > 0 ? ArrayPool<FactorEntry128>.Shared.Rent(counts.Count) : null;
         int index = 0;
         if (rented is not null)
         {
@@ -709,20 +679,20 @@ internal static partial class PrimeOrderCalculator
         }
 
         FactorEntry128[]? resultArray = null;
-        if (pool is not null)
+        if (rented is not null)
         {
             Span<FactorEntry128> span = rented.AsSpan(0, index);
             span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
             resultArray = new FactorEntry128[index];
             span.CopyTo(resultArray);
-            pool.Return(rented!, clearArray: false);
+            ArrayPool<FactorEntry128>.Shared.Return(rented, clearArray: false);
         }
 
         bool fullyFactored = cofactor == UInt128.One;
         return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index);
     }
 
-    private static UInt128 PopulateSmallPrimeFactorsCpuWide(in UInt128 value, uint limit, Dictionary<UInt128, int> counts)
+    private static UInt128 PopulateSmallPrimeFactorsCpuWide(UInt128 value, uint limit, Dictionary<UInt128, int> counts)
     {
         UInt128 remaining = value;
         uint[] primes = PrimesGenerator.SmallPrimes;
@@ -762,7 +732,7 @@ internal static partial class PrimeOrderCalculator
 
         return remaining;
     }
-    private static bool TryPollardRhoWide(in UInt128 n, long deadlineTimestamp, out UInt128 factor)
+    private static bool TryPollardRhoWide(UInt128 n, Stopwatch stopwatch, long budgetTicks, out UInt128 factor)
     {
         factor = UInt128.Zero;
         if ((n & UInt128.One) == UInt128.Zero)
@@ -774,12 +744,10 @@ internal static partial class PrimeOrderCalculator
         UInt128 c = UInt128.One;
         UInt128 x = (UInt128)2UL;
         UInt128 y = x;
-        long timestamp = 0L; // reused for deadline checks.
 
         while (true)
         {
-            timestamp = Stopwatch.GetTimestamp();
-            if (timestamp > deadlineTimestamp)
+            if (stopwatch.ElapsedTicks > budgetTicks)
             {
                 return false;
             }
@@ -806,7 +774,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static UInt128 AdvancePolynomialWide(in UInt128 x, in UInt128 c, in UInt128 modulus)
+    private static UInt128 AdvancePolynomialWide(UInt128 x, UInt128 c, UInt128 modulus)
     {
         BigInteger value = (BigInteger)x;
         value = (value * value + (BigInteger)c) % (BigInteger)modulus;
@@ -844,7 +812,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static int TrailingZeroCountWide(in UInt128 value)
+    private static int TrailingZeroCountWide(UInt128 value)
     {
         if (value == UInt128.Zero)
         {
@@ -900,7 +868,7 @@ internal static partial class PrimeOrderCalculator
         return order;
     }
 
-    private static void FactorCompletelyWide(in UInt128 value, Dictionary<UInt128, int> counts)
+    private static void FactorCompletelyWide(UInt128 value, Dictionary<UInt128, int> counts)
     {
         if (value <= UInt128.One)
         {
@@ -918,7 +886,7 @@ internal static partial class PrimeOrderCalculator
         FactorCompletelyWide(value / factor, counts);
     }
 
-    private static UInt128 PollardRhoWide(in UInt128 n)
+    private static UInt128 PollardRhoWide(UInt128 n)
     {
         if ((n & UInt128.One) == UInt128.Zero)
         {
@@ -951,7 +919,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool IsPrimeWide(in UInt128 value)
+    private static bool IsPrimeWide(UInt128 value)
     {
         if (value <= ulong.MaxValue)
         {
@@ -1031,7 +999,7 @@ internal static partial class PrimeOrderCalculator
         return true;
     }
 
-    private static UInt128 Pow2ModWide(in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData)
+    private static UInt128 Pow2ModWide(UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData divisorData)
     {
         if (modulus == UInt128.One)
         {
@@ -1061,7 +1029,7 @@ internal static partial class PrimeOrderCalculator
         return exponent.Pow2MontgomeryModWindowed(modulus);
     }
 
-    private static void AddFactor(Dictionary<UInt128, int> counts, in UInt128 prime, int exponent)
+    private static void AddFactor(Dictionary<UInt128, int> counts, UInt128 prime, int exponent)
     {
         if (counts.TryGetValue(prime, out int existing))
         {
@@ -1073,30 +1041,60 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private readonly struct CandidateKey128(int primary, bool secondaryDescending, in UInt128 secondary, bool tertiaryDescending, in UInt128 tertiary)
-	{
-		public int Primary { get; } = primary;
+    private readonly struct CandidateKey128
+    {
+        public CandidateKey128(int primary, bool secondaryDescending, UInt128 secondary, bool tertiaryDescending, UInt128 tertiary)
+        {
+            Primary = primary;
+            SecondaryDescending = secondaryDescending;
+            Secondary = secondary;
+            TertiaryDescending = tertiaryDescending;
+            Tertiary = tertiary;
+        }
 
-		public bool SecondaryDescending { get; } = secondaryDescending;
+        public int Primary { get; }
 
-		public UInt128 Secondary { get; } = secondary;
+        public bool SecondaryDescending { get; }
 
-		public bool TertiaryDescending { get; } = tertiaryDescending;
+        public UInt128 Secondary { get; }
 
-		public UInt128 Tertiary { get; } = tertiary;
-	}
+        public bool TertiaryDescending { get; }
 
-    private readonly struct PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count)
-	{
-		public FactorEntry128[]? Factors { get; } = factors;
+        public UInt128 Tertiary { get; }
+    }
 
-		public UInt128 Cofactor { get; } = cofactor;
+    private readonly struct FactorEntry128
+    {
+        public FactorEntry128(UInt128 value, int exponent)
+        {
+            Value = value;
+            Exponent = exponent;
+        }
 
-		public bool FullyFactored { get; } = fullyFactored;
+        public UInt128 Value { get; }
 
-		public int Count { get; } = count;
+        public int Exponent { get; }
+    }
 
-		public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0);
+    private readonly struct PartialFactorResult128
+    {
+        public PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count)
+        {
+            Factors = factors;
+            Cofactor = cofactor;
+            FullyFactored = fullyFactored;
+            Count = count;
+        }
+
+        public FactorEntry128[]? Factors { get; }
+
+        public UInt128 Cofactor { get; }
+
+        public bool FullyFactored { get; }
+
+        public int Count { get; }
+
+        public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0);
 
         public PartialFactorResult128 WithAdditionalPrime(UInt128 prime)
         {

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/ThreadStaticPools.cs

@@ -1,320 +0,0 @@
-using System.Buffers;
-using System.Runtime.CompilerServices;
-using PerfectNumbers.Core.Cpu;
-
-namespace PerfectNumbers.Core
-{
-    public readonly struct ThreadStaticPools
-    {
-        [ThreadStatic]
-        private static ArrayPool<FactorEntry>? _factorEntryPool;
-
-        public static ArrayPool<FactorEntry> FactorEntryPool
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get
-            {
-                return _factorEntryPool ??= ArrayPool<FactorEntry>.Create();
-            }
-        }
-
-        [ThreadStatic]
-        private static ArrayPool<FactorEntry128>? _factorEntry128Pool;
-
-        public static ArrayPool<FactorEntry128> FactorEntry128Pool
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get
-            {
-                return _factorEntry128Pool ??= ArrayPool<FactorEntry128>.Create();
-            }
-        }
-
-
-        [ThreadStatic]
-        private static ArrayPool<KeyValuePair<ulong, int>>? _factorKeyValuePairPool;
-
-        public static ArrayPool<KeyValuePair<ulong, int>> FactorKeyValuePairPool
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get
-            {
-                return _factorKeyValuePairPool ??= ArrayPool<KeyValuePair<ulong, int>>.Create();
-            }
-        }
-        [ThreadStatic]
-        private static MersenneCpuDivisorScanSession? _mersenneCpuDivisorSession;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static MersenneCpuDivisorScanSession? RentMersenneCpuDivisorSession()
-        {
-            MersenneCpuDivisorScanSession? session = _mersenneCpuDivisorSession;
-            if (session is null)
-            {
-                return null;
-            }
-
-            _mersenneCpuDivisorSession = null;
-            session.Reset();
-            return session;
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static void ReturnMersenneCpuDivisorSession(MersenneCpuDivisorScanSession session)
-        {
-            _mersenneCpuDivisorSession = session;
-        }
-
-        [ThreadStatic]
-        private static ArrayPool<ulong>? _ulongPool;
-
-        public static ArrayPool<ulong> UlongPool
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get
-            {
-                return _ulongPool ??= ArrayPool<ulong>.Create();
-            }
-        }
-
-        [ThreadStatic]
-        private static ArrayPool<int>? _intPool;
-
-        public static ArrayPool<int> IntPool
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get
-            {
-                return _intPool ??= ArrayPool<int>.Create();
-            }
-        }
-
-
-        [ThreadStatic]
-        private static List<List<ulong>>? _ulongListPool;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static List<ulong> RentUlongList(int capacityHint)
-        {
-            List<List<ulong>>? pool = _ulongListPool;
-            if (pool is not null && pool.Count > 0)
-            {
-                int lastIndex = pool.Count - 1;
-                List<ulong> list = pool[lastIndex];
-                pool.RemoveAt(lastIndex);
-                if (list.Count > 0)
-                {
-                    list.Clear();
-                }
-
-                if (list.Capacity < capacityHint)
-                {
-                    list.Capacity = capacityHint;
-                }
-
-                return list;
-            }
-
-            return new List<ulong>(capacityHint);
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static void ReturnUlongList(List<ulong> list)
-        {
-            List<List<ulong>>? pool = _ulongListPool;
-            if (pool is null)
-            {
-                pool = new List<List<ulong>>(4);
-                _ulongListPool = pool;
-            }
-
-            pool.Add(list);
-        }
-
-        [ThreadStatic]
-        private static List<Stack<ulong>>? _ulongStackPool;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static Stack<ulong> RentUlongStack(int capacityHint)
-        {
-            List<Stack<ulong>>? pool = _ulongStackPool;
-            if (pool is not null && pool.Count > 0)
-            {
-                int lastIndex = pool.Count - 1;
-                Stack<ulong> stack = pool[lastIndex];
-                pool.RemoveAt(lastIndex);
-                if (stack.Count > 0)
-                {
-                    stack.Clear();
-                }
-
-                if (stack.Count < capacityHint)
-                {
-                    stack.EnsureCapacity(capacityHint);
-                }
-
-                return stack;
-            }
-
-            return new Stack<ulong>(capacityHint);
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static void ReturnUlongStack(Stack<ulong> stack)
-        {
-            List<Stack<ulong>>? pool = _ulongStackPool;
-            if (pool is null)
-            {
-                pool = new List<Stack<ulong>>(4);
-                _ulongStackPool = pool;
-            }
-
-            pool.Add(stack);
-        }
-
-        [ThreadStatic]
-        private static List<Dictionary<ulong, int>>? _ulongIntDictionaryPool;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static Dictionary<ulong, int> RentUlongIntDictionary(int capacityHint)
-        {
-            List<Dictionary<ulong, int>>? pool = _ulongIntDictionaryPool;
-            if (pool is not null && pool.Count > 0)
-            {
-                int lastIndex = pool.Count - 1;
-                Dictionary<ulong, int> dictionary = pool[lastIndex];
-                pool.RemoveAt(lastIndex);
-                if (dictionary.Count > 0)
-                {
-                    dictionary.Clear();
-                }
-
-                if (dictionary.Count < capacityHint)
-                {
-                    dictionary.EnsureCapacity(capacityHint);
-                }
-
-                return dictionary;
-            }
-
-            return new Dictionary<ulong, int>(capacityHint);
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static void ReturnUlongIntDictionary(Dictionary<ulong, int> dictionary)
-        {
-            List<Dictionary<ulong, int>>? pool = _ulongIntDictionaryPool;
-            if (pool is null)
-            {
-                pool = new List<Dictionary<ulong, int>>(4);
-                _ulongIntDictionaryPool = pool;
-            }
-
-            pool.Add(dictionary);
-        }
-
-        [ThreadStatic]
-        private static Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? _mersenneFactorCacheDictionary;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry> RentMersenneFactorCacheDictionary()
-        {
-            Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? cache = _mersenneFactorCacheDictionary;
-            if (cache is null)
-            {
-                return new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
-            }
-
-            _mersenneFactorCacheDictionary = null;
-            cache.Clear();
-            return cache;
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static void ReturnMersenneFactorCacheDictionary(Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry> cache)
-        {
-            for (Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>.Enumerator enumerator = cache.GetEnumerator(); enumerator.MoveNext();)
-            {
-                enumerator.Current.Value.ReturnToPool();
-            }
-
-            cache.Clear();
-            _mersenneFactorCacheDictionary = cache;
-        }
-
-        [ThreadStatic]
-        private static MersenneDivisorCycles.FactorCacheEntry? _factorCacheEntryPoolHead;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static MersenneDivisorCycles.FactorCacheEntry RentFactorCacheEntry()
-        {
-            MersenneDivisorCycles.FactorCacheEntry? entry = _factorCacheEntryPoolHead;
-            if (entry is null)
-            {
-                entry = new MersenneDivisorCycles.FactorCacheEntry();
-            }
-            else
-            {
-                _factorCacheEntryPoolHead = entry.Next;
-            }
-
-            entry.Reset();
-            return entry;
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static void ReturnFactorCacheEntry(MersenneDivisorCycles.FactorCacheEntry entry)
-        {
-            entry.Next = _factorCacheEntryPoolHead;
-            _factorCacheEntryPoolHead = entry;
-        }
-
-        [ThreadStatic]
-        private static Dictionary<ulong, int>? _factorCountDictionary;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static Dictionary<ulong, int> RentFactorCountDictionary()
-        {
-            Dictionary<ulong, int>? dictionary = _factorCountDictionary;
-            if (dictionary is null)
-            {
-                return new Dictionary<ulong, int>(8);
-            }
-
-            _factorCountDictionary = null;
-            dictionary.Clear();
-            return dictionary;
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static void ReturnFactorCountDictionary(Dictionary<ulong, int> dictionary)
-        {
-            dictionary.Clear();
-            _factorCountDictionary = dictionary;
-        }
-
-        [ThreadStatic]
-        private static Dictionary<ulong, int>? _factorScratchDictionary;
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static Dictionary<ulong, int> RentFactorScratchDictionary()
-        {
-            Dictionary<ulong, int>? dictionary = _factorScratchDictionary;
-            if (dictionary is null)
-            {
-                return new Dictionary<ulong, int>(8);
-            }
-
-            _factorScratchDictionary = null;
-            dictionary.Clear();
-            return dictionary;
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static void ReturnFactorScratchDictionary(Dictionary<ulong, int> dictionary)
-        {
-            dictionary.Clear();
-            _factorScratchDictionary = dictionary;
-        }
-    }
-}

PerfectNumbers.Core/UInt128Extensions.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using PerfectNumbers.Core.Gpu;

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`.
+

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.