Switch single-use GPU exponent checks to standard modular math

Clarify per-divisor Montgomery recomputation

EvenPerfectBitScanner.Benchmarks/EvenPerfectBitScanner.Benchmarks.csproj

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

EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/MutableUInt128AddBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/MutableUInt128Mul64Benchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/MutableUInt128MultiplyAddBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/MutableUInt128SubMutableBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs

@@ -22,320 +22,391 @@ namespace EvenPerfectBitScanner.Benchmarks;
 [SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
 public class Pow2MontgomeryModCycleComputationBenchmarks
 {
-    private const int SampleCount = 256;
-    private const ulong OneShiftedLeft63 = 1UL << 63;
-    private const long OneShiftedLeft60 = 1L << 60;
-    private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
-    private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
-    private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount];
-    private readonly UInt128[] _wideModuli = new UInt128[SampleCount];
-    private readonly ulong[] _smallExponents = new ulong[SampleCount];
-    private readonly ulong[] _largeExponents = new ulong[SampleCount];
-    private readonly ulong[] _veryLargeExponents = new ulong[SampleCount];
-    private readonly ulong[] _smallCycles = new ulong[SampleCount];
-    private readonly ulong[] _largeCycles = new ulong[SampleCount];
-    private readonly ulong[] _veryLargeCycles = new ulong[SampleCount];
-    private readonly UInt128[] _wideCycles = new UInt128[SampleCount];
+	private const int SampleCount = 256;
+	private const ulong OneShiftedLeft63 = 1UL << 63;
+	private const long OneShiftedLeft60 = 1L << 60;
+	private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount];
+	private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount];
+	private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount];
+	private readonly UInt128[] _wideModuli = new UInt128[SampleCount];
+	private readonly ulong[] _smallExponents = new ulong[SampleCount];
+	private readonly ulong[] _largeExponents = new ulong[SampleCount];
+	private readonly ulong[] _veryLargeExponents = new ulong[SampleCount];
+	private readonly ulong[] _smallCycles = new ulong[SampleCount];
+	private readonly ulong[] _largeCycles = new ulong[SampleCount];
+	private readonly ulong[] _veryLargeCycles = new ulong[SampleCount];
+	private readonly UInt128[] _wideCycles = new UInt128[SampleCount];
 
-    private readonly Random _random = new(113);
-    private ulong? _previousPrimeOrder;
-    private UInt128? _previousWidePrimeOrder;
+	private readonly Random _random = new(113);
+	private ulong? _previousPrimeOrder;
+	private UInt128? _previousWidePrimeOrder;
 
-    public enum InputScale
-    {
-        Small,
-        Large,
-        VeryLarge
-    }
+	public enum InputScale
+	{
+		Small,
+		Large,
+		VeryLarge
+	}
 
-    [ParamsAllValues]
-    public InputScale Scale { get; set; }
+	[ParamsAllValues]
+	public InputScale Scale { get; set; }
 
-    [GlobalSetup]
-    public void Setup()
-    {
-        for (int i = 0; i < SampleCount; i++)
-        {
+	[GlobalSetup]
+	public void Setup()
+	{
+		for (int i = 0; i < SampleCount; i++)
+		{
 #if DEBUG
-            Console.WriteLine($"Generating sample {i}");
+			Console.WriteLine($"Generating sample {i}");
 #endif
-            ulong smallModulus = NextSmallOddModulus();
-            _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
-            _smallExponents[i] = NextSmallExponent();
-            _smallCycles[i] = CalculateCycleLengthWithHeuristics(smallModulus);
+			ulong smallModulus = NextSmallOddModulus();
+			_smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus);
+			_smallExponents[i] = NextSmallExponent();
+			_smallCycles[i] = CalculateCycleLengthWithHeuristics(smallModulus);
 
-            (ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
-            _largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
-            _largeExponents[i] = NextLargeExponent();
-            _largeCycles[i] = CalculateCycleLengthWithHeuristics(largeModulus, largeCycle);
+			(ulong largeModulus, ulong largeCycle) = NextLargeModulusAndCycle();
+			_largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus);
+			_largeExponents[i] = NextLargeExponent();
+			_largeCycles[i] = CalculateCycleLengthWithHeuristics(largeModulus, largeCycle);
 
-            ulong veryLargeModulus = NextVeryLargeOddModulus();
-            _veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus);
-            _veryLargeExponents[i] = NextVeryLargeExponent();
+			ulong veryLargeModulus = NextVeryLargeOddModulus();
+			_veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus);
+			_veryLargeExponents[i] = NextVeryLargeExponent();
 #if DEBUG
-            Console.WriteLine($"Calculating cycle length {veryLargeModulus}");
+			Console.WriteLine($"Calculating cycle length {veryLargeModulus}");
 #endif
-            _veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus);
+			_veryLargeCycles[i] = CalculateCycleLengthWithHeuristics(veryLargeModulus);
 
-            Console.WriteLine($"Generating wide sample {i}");
-            UInt128 wideModulus = NextWideOddModulus();
-            _wideModuli[i] = wideModulus;
-            _wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus);
-        }
-    }
+			Console.WriteLine($"Generating wide sample {i}");
+			UInt128 wideModulus = NextWideOddModulus();
+			_wideModuli[i] = wideModulus;
+			_wideCycles[i] = CalculateCycleLengthWithHeuristics(wideModulus);
+		}
+	}
 
-    /// <summary>
-    /// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
-    /// large set.
-    /// </summary>
-    // [Benchmark(Baseline = true)]
-    public ulong MontgomeryWithoutCycle()
-    {
-        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
-        ulong checksum = 0UL;
+	/// <summary>
+	/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
+	/// large set.
+	/// </summary>
+	// [Benchmark(Baseline = true)]
+	public ulong MontgomeryWithoutCycleCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
 
-        for (int i = 0; i < SampleCount; i++)
-        {
-            checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false);
-        }
+		for (int i = 0; i < SampleCount; i++)
+		{
+			checksum ^= exponents[i].Pow2MontgomeryModWindowedCpu(divisors[i], keepMontgomery: false);
+		}
 
-        return checksum;
-    }
+		return checksum;
+	}
 
-    /// <summary>
-    /// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
-    /// on the large one (0.08×).
-    /// </summary>
-    [Benchmark]
-    public ulong MontgomeryWithPrecomputedCycle()
-    {
-        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
-        ulong checksum = 0UL;
+	/// <summary>
+	/// Baseline Montgomery reduction without precomputed cycle data; measured 35.41 μs on the small sample and 117.33 μs on the
+	/// large set.
+	/// </summary>
+	// [Benchmark(Baseline = true)]
+	public ulong MontgomeryWithoutCycleGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
 
-        for (int i = 0; i < SampleCount; i++)
-        {
-            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycles[i], divisors[i]);
-        }
+		for (int i = 0; i < SampleCount; i++)
+		{
+			checksum ^= exponents[i].Pow2MontgomeryModWindowedGpu(divisors[i], keepMontgomery: false);
+		}
 
-        return checksum;
-    }
+		return checksum;
+	}
 
-    /// <summary>
-    /// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
-    /// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
-    /// </summary>
-    // [Benchmark]
-    public ulong MontgomeryWithGpuCycleComputation()
-    {
-        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
-        ulong checksum = 0UL;
+	/// <summary>
+	/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
+	/// on the large one (0.08×).
+	/// </summary>
+	[Benchmark]
+	public ulong MontgomeryWithPrecomputedCycleCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
+		ulong checksum = 0UL;
 
-        for (int i = 0; i < SampleCount; i++)
-        {
-            ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
-            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
-        }
+		for (int i = 0; i < SampleCount; i++)
+		{
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycles[i], divisors[i]);
+		}
 
-        return checksum;
-    }
+		return checksum;
+	}
+
+	/// <summary>
+	/// Montgomery reduction with a known cycle length (cycle lookup only); 25.81 μs on the small set (0.73× baseline) and 9.12 μs
+	/// on the large one (0.08×).
+	/// </summary>
+	[Benchmark]
+	public ulong MontgomeryWithPrecomputedCycleGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycles[i], divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	/// <summary>
+	/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
+	/// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
+	/// </summary>
+	// [Benchmark]
+	public ulong MontgomeryWithGpuCycleComputationCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	/// <summary>
+	/// Montgomery reduction with the cycle computed on the fly using the GPU helper; costs 54.09 ms on the small benchmark due to
+	/// cycle discovery, but drops to 12.91 μs on the large set (0.11× baseline).
+	/// </summary>
+	// [Benchmark]
+	public ulong MontgomeryWithGpuCycleComputationGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	[Benchmark]
+	public ulong MontgomeryWithHeuristicCycleComputationCpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleCpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	[Benchmark]
+	public ulong MontgomeryWithHeuristicCycleComputationGpu()
+	{
+		GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
+		ulong checksum = 0UL;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
+			checksum ^= exponents[i].Pow2MontgomeryModWithCycleGpu(cycle, divisors[i]);
+		}
+
+		return checksum;
+	}
+
+	[Benchmark]
+	public UInt128 HeuristicWideCycleComputation()
+	{
+		UInt128 checksum = UInt128.Zero;
+
+		for (int i = 0; i < SampleCount; i++)
+		{
+			UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
+			_wideCycles[i] = cycle;
+			checksum ^= cycle;
+		}
+
+		return checksum;
+	}
+
+	private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
+	{
+		switch (Scale)
+		{
+			case InputScale.Small:
+				exponents = _smallExponents;
+				divisors = _smallDivisors;
+				cycles = _smallCycles;
+				break;
+			case InputScale.Large:
+				exponents = _largeExponents;
+				divisors = _largeDivisors;
+				cycles = _largeCycles;
+				break;
+			case InputScale.VeryLarge:
+				exponents = _veryLargeExponents;
+				divisors = _veryLargeDivisors;
+				cycles = _veryLargeCycles;
+				break;
+			default:
+				throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null);
+		}
+	}
+
+	private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
+
+	private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue);
+
+	private ulong NextVeryLargeExponent()
+	{
+		Span<byte> buffer = stackalloc byte[8];
+		_random.NextBytes(buffer);
+		ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
+		return value;
+	}
+
+	private ulong NextSmallOddModulus()
+	{
+		ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
+		return value | 1UL;
+	}
+
+	private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
+	{
+		int bitLength = _random.Next(48, 64);
+		ulong modulus = (1UL << bitLength) - 1UL;
+		return (modulus, (ulong)bitLength);
+	}
+
+	private ulong NextVeryLargeOddModulus()
+	{
+		Span<byte> buffer = stackalloc byte[8];
+		_random.NextBytes(buffer);
+		ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
+		return value;
+	}
+
+	private UInt128 NextWideOddModulus()
+	{
+		Span<byte> buffer = stackalloc byte[16];
+		_random.NextBytes(buffer);
+		ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
+		Span<byte> highSlice = buffer.Slice(8, 8);
+		ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
+		UInt128 value = ((UInt128)high << 64) | low;
+		return value | UInt128.One;
+	}
+
+	private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
+	{
+		if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+		{
+			return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL);
+		}
+
+		ulong nPrime = ComputeMontgomeryNPrime(modulus);
+		ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
+		ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
+		ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
+
+		return new MontgomeryDivisorData(
+			modulus,
+			nPrime,
+			montgomeryOne,
+			montgomeryTwo,
+			montgomeryTwoSquared);
+	}
 
 
-    [Benchmark]
-    public ulong MontgomeryWithHeuristicCycleComputation()
-    {
-        GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out _);
-        ulong checksum = 0UL;
+	private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
+	{
+		MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
+		if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+		{
+			return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
+		}
 
-        for (int i = 0; i < SampleCount; i++)
-        {
-            ulong cycle = CalculateCycleLengthWithHeuristics(divisors[i].Modulus);
-            checksum ^= exponents[i].Pow2MontgomeryModWithCycle(cycle, divisors[i]);
-        }
+		// When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
 
-        return checksum;
-    }
-
-    [Benchmark]
-    public UInt128 HeuristicWideCycleComputation()
-    {
-        UInt128 checksum = UInt128.Zero;
-
-        for (int i = 0; i < SampleCount; i++)
-        {
-            UInt128 cycle = CalculateCycleLengthWithHeuristics(_wideModuli[i]);
-            _wideCycles[i] = cycle;
-            checksum ^= cycle;
-        }
-
-        return checksum;
-    }
-
-    private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles)
-    {
-        switch (Scale)
-        {
-            case InputScale.Small:
-                exponents = _smallExponents;
-                divisors = _smallDivisors;
-                cycles = _smallCycles;
-                break;
-            case InputScale.Large:
-                exponents = _largeExponents;
-                divisors = _largeDivisors;
-                cycles = _largeCycles;
-                break;
-            case InputScale.VeryLarge:
-                exponents = _veryLargeExponents;
-                divisors = _veryLargeDivisors;
-                cycles = _veryLargeCycles;
-                break;
-            default:
-                throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null);
-        }
-    }
-
-    private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L);
-
-    private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue);
-
-    private ulong NextVeryLargeExponent()
-    {
-        Span<byte> buffer = stackalloc byte[8];
-        _random.NextBytes(buffer);
-        ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63);
-        return value;
-    }
-
-    private ulong NextSmallOddModulus()
-    {
-        ulong value = (ulong)_random.NextInt64(3L, 1_000_000L);
-        return value | 1UL;
-    }
-
-    private (ulong modulus, ulong cycleLength) NextLargeModulusAndCycle()
-    {
-        int bitLength = _random.Next(48, 64);
-        ulong modulus = (1UL << bitLength) - 1UL;
-        return (modulus, (ulong)bitLength);
-    }
-
-    private ulong NextVeryLargeOddModulus()
-    {
-        Span<byte> buffer = stackalloc byte[8];
-        _random.NextBytes(buffer);
-        ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL;
-        return value;
-    }
-
-    private UInt128 NextWideOddModulus()
-    {
-        Span<byte> buffer = stackalloc byte[16];
-        _random.NextBytes(buffer);
-        ulong low = BinaryPrimitives.ReadUInt64LittleEndian(buffer);
-        Span<byte> highSlice = buffer.Slice(8, 8);
-        ulong high = BinaryPrimitives.ReadUInt64LittleEndian(highSlice) | (1UL << 63);
-        UInt128 value = ((UInt128)high << 64) | low;
-        return value | UInt128.One;
-    }
-
-    private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
-    {
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL);
-        }
-
-        ulong nPrime = ComputeMontgomeryNPrime(modulus);
-        ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
-        ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus);
-        ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
-
-        return new MontgomeryDivisorData(
-            modulus,
-            nPrime,
-            montgomeryOne,
-            montgomeryTwo,
-            montgomeryTwoSquared);
-    }
-
-
-    private ulong CalculateCycleLengthWithHeuristics(ulong modulus, ulong fallbackOrder = 0UL)
-    {
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(modulus);
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
-        }
-
-        // When the below is uncommented, the benchmark never completes due to the primality test cost on very large random set.
-
-        // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(modulus);
-        // // bool isPrime = PrimeTester.IsPrimeInternal(modulus, CancellationToken.None);
-        // if (!isPrime)
-        // {
-        //     return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
-        // }
+		// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(modulus);
+		// // bool isPrime = PrimeTester.IsPrimeInternal(modulus, CancellationToken.None);
+		// if (!isPrime)
+		// {
+		//     return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus);
+		// }
 
 #if DEBUG
-        Console.WriteLine("Trying heuristic. Prime order calculation");
+		Console.WriteLine("Trying heuristic. Prime order calculation");
 #endif
-        PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
-            modulus,
-            _previousPrimeOrder,
-            divisorData,
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+		ulong order = PrimeOrderCalculator.Calculate(
+			modulus,
+			_previousPrimeOrder,
+			divisorData,
+			PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+			PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        if (orderResult.Order != 0UL)
-        {
-            _previousPrimeOrder = orderResult.Order;
-            return orderResult.Order;
-        }
+		if (order != 0UL)
+		{
+			_previousPrimeOrder = order;
+			return order;
+		}
 
-        _previousPrimeOrder = null;
+		_previousPrimeOrder = null;
 
-        Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation");
-        return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
-    }
+		Console.WriteLine($"Heuristic failed for {modulus}, falling back to full cycle calculation");
+		return fallbackOrder != 0UL ? fallbackOrder : MersenneDivisorCycles.CalculateCycleLength(modulus, divisorData);
+	}
 
-    private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
-    {
-        if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero)
-        {
-            return MersenneDivisorCycles.GetCycle(modulus);
-        }
+	private UInt128 CalculateCycleLengthWithHeuristics(UInt128 modulus)
+	{
+		if (modulus <= UInt128.One || (modulus & UInt128.One) == UInt128.Zero)
+		{
+			return MersenneDivisorCycles.GetCycle(modulus);
+		}
 
 #if DEBUG
-        Console.WriteLine("Trying heuristic. Wide prime order calculation");
+		Console.WriteLine("Trying heuristic. Wide prime order calculation");
 #endif
-        PrimeOrderCalculator.PrimeOrderResultWide orderResult = PrimeOrderCalculator.Calculate(
-            modulus,
-            _previousWidePrimeOrder,
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+        UInt128 order = PrimeOrderCalculator.Calculate(
+                modulus,
+                _previousWidePrimeOrder,
+                PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+                PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        if (orderResult.Order != UInt128.Zero)
-        {
-            _previousWidePrimeOrder = orderResult.Order;
-            return orderResult.Order;
-        }
+		if (order != UInt128.Zero)
+		{
+			_previousWidePrimeOrder = order;
+			return order;
+		}
 
-        _previousWidePrimeOrder = null;
+		_previousWidePrimeOrder = null;
 
-        Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation");
-        return MersenneDivisorCycles.GetCycle(modulus);
-    }
+		Console.WriteLine($"Wide heuristic failed for {modulus}, falling back to full cycle calculation");
+		return MersenneDivisorCycles.GetCycle(modulus);
+	}
 
-    private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
+	private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
 
-    private static ulong ComputeMontgomeryNPrime(ulong modulus)
-    {
-        ulong inv = modulus;
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        inv *= unchecked(2UL - modulus * inv);
-        return unchecked(0UL - inv);
-    }
+	private static ulong ComputeMontgomeryNPrime(ulong modulus)
+	{
+		ulong inv = modulus;
+		inv *= unchecked(2UL - modulus * inv);
+		inv *= unchecked(2UL - modulus * inv);
+		inv *= unchecked(2UL - modulus * inv);
+		inv *= unchecked(2UL - modulus * inv);
+		inv *= unchecked(2UL - modulus * inv);
+		inv *= unchecked(2UL - modulus * inv);
+		return unchecked(0UL - inv);
+	}
 }
 

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Program.cs

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

EvenPerfectBitScanner.ResultsParser/EvenPerfectBitScanner.ResultsParser.csproj

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

EvenPerfectBitScanner.Tests/EvenPerfectBitScanner.Tests.csproj

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

EvenPerfectBitScanner/EvenPerfectBitScanner.csproj

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

EvenPerfectBitScanner/Program.cs

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

PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs

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

PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs

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

PerfectNumbers.Core.Tests/Gpu/NttGpuMathTests.cs

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

PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs

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

PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs

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

PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs

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

PerfectNumbers.Core.Tests/MersenneNumberDivisorByDivisorTesterTests.cs

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

PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs

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

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

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

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

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

PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs

@@ -7,6 +7,9 @@ namespace PerfectNumbers.Core.Tests;
 
 public class PrimeOrderCalculatorTests
 {
+    // EvenPerfectBitScanner never exercises the GPU heuristic device or the small-prime guards that the legacy tests cover,
+    // so those scenarios are commented out because they are unreachable in the production pipeline.
+    /*
     public static IEnumerable<object[]> SmallPrimes()
     {
         yield return new object[] { 3UL };
@@ -45,14 +48,14 @@ public class PrimeOrderCalculatorTests
     {
         ulong expected = ComputeOrderByDoubling(prime);
 
-        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
+        ulong result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
-        result.Order.Should().Be(expected);
+        result.Should().Be(expected);
     }
 
     [Theory]
@@ -62,21 +65,21 @@ public class PrimeOrderCalculatorTests
     {
         ulong previousOrder = ComputeOrderByDoubling(previousPrime);
 
-        PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
+        ulong heuristic = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
+        ulong strict = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault);
+            PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault,
+			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
-        strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
-        heuristic.Order.Should().Be(strict.Order);
+        heuristic.Should().Be(strict);
     }
 
     [Theory]
@@ -85,14 +88,14 @@ public class PrimeOrderCalculatorTests
     [InlineData(3UL, 2UL)]
     public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder)
     {
-        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
+        ulong result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
-        result.Order.Should().Be(expectedOrder);
+        result.Should().Be(expectedOrder);
     }
 
     [Fact]
@@ -105,14 +108,14 @@ public class PrimeOrderCalculatorTests
             maxPowChecks: 8,
             mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
 
-        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
+        ulong result = PrimeOrderCalculator.Calculate(
             prime: 13UL,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(13UL),
-            config);
+            config,
+			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
-        result.Order.Should().Be(ComputeOrderByDoubling(13UL));
+        result.Should().Be(ComputeOrderByDoubling(13UL));
     }
 
     [Fact]
@@ -126,14 +129,14 @@ public class PrimeOrderCalculatorTests
             mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
 
         MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
-        PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
+        ulong heuristic = PrimeOrderCalculator.Calculate(
             prime: 239UL,
             previousOrder: null,
             divisorData,
-            heuristicConfig);
+            heuristicConfig,
+			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
-        heuristic.Order.Should().Be(ComputeOrderByDoubling(239UL));
+        heuristic.Should().Be(ComputeOrderByDoubling(239UL));
 
         var strictConfig = new PrimeOrderCalculator.PrimeOrderSearchConfig(
             smallFactorLimit: 2,
@@ -141,14 +144,14 @@ public class PrimeOrderCalculatorTests
             maxPowChecks: 1,
             mode: PrimeOrderCalculator.PrimeOrderMode.Strict);
 
-        PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
+        ulong strict = PrimeOrderCalculator.Calculate(
             prime: 239UL,
             previousOrder: null,
             divisorData,
-            strictConfig);
+            strictConfig,
+			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
-        strict.Order.Should().Be(heuristic.Order);
+        strict.Should().Be(heuristic);
     }
 
     [Fact]
@@ -157,13 +160,13 @@ public class PrimeOrderCalculatorTests
     {
         UInt128 prime = UInt128.Parse("18446744073709641691");
 
-        PrimeOrderCalculator.PrimeOrderResultWide result = PrimeOrderCalculator.Calculate(
-            prime,
-            previousOrder: null,
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+        UInt128 result = PrimeOrderCalculator.Calculate(
+                prime,
+                previousOrder: null,
+                PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+                PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
-        result.Order.Should().Be((UInt128)1229782938247309446UL);
+        result.Should().Be((UInt128)1229782938247309446UL);
     }
 
     [Fact]
@@ -197,4 +200,5 @@ public class PrimeOrderCalculatorTests
 
         return order;
     }
+*/
 }

PerfectNumbers.Core/ByDivisorDeltasDevice.cs

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

PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs

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

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -1,8 +1,10 @@
 using System;
+using System.Buffers;
 using System.Collections.Concurrent;
 using System.Collections.Generic;
 using System.Numerics;
 using PerfectNumbers.Core;
+using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core.Cpu;
 
@@ -15,17 +17,49 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
     private bool _isConfigured;
     private int _batchSize = 1_024;
 
+    private ulong _configuredMaxPrime;
+    private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu;
+    private ByDivisorMontgomeryDevice _montgomeryDevice = ByDivisorMontgomeryDevice.Cpu;
 
     public int BatchSize
     {
         get => _batchSize;
-        set => _batchSize = Math.Max(1, value);
+        set
+        {
+            int sanitized = Math.Max(1, value);
+            _batchSize = sanitized;
+        }
+    }
+
+    public ByDivisorDeltasDevice DeltasDevice
+    {
+        get => _deltasDevice;
+        set
+        {
+            lock (_sync)
+            {
+                _deltasDevice = value;
+            }
+        }
+    }
+
+    public ByDivisorMontgomeryDevice MontgomeryDevice
+    {
+        get => _montgomeryDevice;
+        set
+        {
+            lock (_sync)
+            {
+                _montgomeryDevice = value;
+            }
+        }
     }
 
     public void ConfigureFromMaxPrime(ulong maxPrime)
     {
         lock (_sync)
         {
+            _configuredMaxPrime = maxPrime;
             _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
             _lastStatusDivisor = 0UL;
             _isConfigured = true;
@@ -61,7 +95,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         }
     }
 
-    public bool IsPrime(ulong prime, out bool divisorsExhausted)
+    public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null)
     {
         ulong allowedMax;
         lock (_sync)
@@ -80,16 +114,17 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             return true;
         }
 
-        ulong processedCount = 0UL;
-        ulong lastProcessed = 0UL;
-        bool processedAll = false;
+        ulong processedCount;
+        ulong lastProcessed;
+        bool processedAll;
 
         bool composite = CheckDivisors(
             prime,
             allowedMax,
             out lastProcessed,
             out processedAll,
-            out processedCount);
+            out processedCount,
+            timeLimit);
 
         if (processedCount > 0UL)
         {
@@ -163,12 +198,38 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         ulong allowedMax,
         out ulong lastProcessed,
         out bool processedAll,
-        out ulong processedCount)
+        out ulong processedCount,
+        TimeSpan? timeLimit)
+    {
+        if (_deltasDevice == ByDivisorDeltasDevice.Cpu)
+        {
+            return CheckDivisorsCpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit);
+        }
+
+        return CheckDivisorsGpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit);
+    }
+
+    private bool CheckDivisorsCpu(
+        ulong prime,
+        ulong allowedMax,
+        out ulong lastProcessed,
+        out bool processedAll,
+        out ulong processedCount,
+        TimeSpan? timeLimit)
     {
         lastProcessed = 0UL;
         processedCount = 0UL;
         processedAll = false;
 
+        PrimeTestTimeLimit limitGuard;
+        if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard))
+        {
+            return false;
+        }
+
+        bool enforceLimit = limitGuard.IsActive;
+        ulong iterationCounter = 0UL;
+
         if (allowedMax < 3UL)
         {
             return false;
@@ -212,6 +273,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
         while (divisor <= limit)
         {
+            if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired())
+            {
+                processedAll = false;
+                return false;
+            }
+
+            iterationCounter++;
             ulong candidate = (ulong)divisor;
             processedCount++;
             lastProcessed = candidate;
@@ -222,7 +290,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
             if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0)
             {
-                MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(candidate);
+                // Each divisor appears only once in this monotonically increasing sequence, so
+                // rebuilding the Montgomery data per candidate avoids caching entries that we
+                // would never hit again.
+                MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
                 ulong divisorCycle;
 
                 if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
@@ -232,6 +303,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 else
                 {
                     factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
+                    // EvenPerfectBitScanner reaches this path only for divisors above the precomputed cycle snapshot.
+                    // Those candidates follow q = 2 * prime * k + 1 with prime >= 5, so the downstream order calculator
+                    // can rely on an odd divisor and a prime exponent when deriving the cycle length.
                     if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL)
                     {
                         divisorCycle = cycleCache.GetCycleLength(candidate);
@@ -266,9 +340,248 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         return false;
     }
 
+
+    private bool CheckDivisorsGpu(
+        ulong prime,
+        ulong allowedMax,
+        out ulong lastProcessed,
+        out bool processedAll,
+        out ulong processedCount,
+        TimeSpan? timeLimit)
+    {
+        lastProcessed = 0UL;
+        processedCount = 0UL;
+        processedAll = false;
+
+        PrimeTestTimeLimit limitGuard;
+        if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard))
+        {
+            return false;
+        }
+
+        bool enforceLimit = limitGuard.IsActive;
+        ulong iterationCounter = 0UL;
+
+        if (allowedMax < 3UL)
+        {
+            return false;
+        }
+
+        GpuUInt128 step128 = new GpuUInt128(prime);
+        step128.ShiftLeft(1);
+        if (step128.IsZero)
+        {
+            processedAll = true;
+            return false;
+        }
+
+        GpuUInt128 limit128 = new GpuUInt128(allowedMax);
+        GpuUInt128 divisor128 = step128;
+        divisor128.Add(1UL);
+        if (divisor128.CompareTo(limit128) > 0)
+        {
+            processedAll = true;
+            return false;
+        }
+
+        ulong step = step128.Low;
+        ulong limit64 = allowedMax;
+
+        int batchCapacity = Math.Max(1, _batchSize);
+        DivisorCycleCache cycleCache = DivisorCycleCache.Shared;
+        int preferredCycleBatch = cycleCache.PreferredBatchSize;
+        if (preferredCycleBatch > 0 && preferredCycleBatch < batchCapacity)
+        {
+            batchCapacity = preferredCycleBatch;
+        }
+
+        var candidateEvaluator = new MersenneNumberDivisorCandidateGpuEvaluator(batchCapacity);
+
+        ulong[] candidateBuffer = ArrayPool<ulong>.Shared.Rent(batchCapacity);
+        byte[] maskBuffer = ArrayPool<byte>.Shared.Rent(batchCapacity);
+
+        Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
+
+        const int RemainderTableSize = 6;
+        byte[] remainderBuffer = ArrayPool<byte>.Shared.Rent(RemainderTableSize);
+        byte[] stepBuffer = ArrayPool<byte>.Shared.Rent(RemainderTableSize);
+        Span<byte> remainderSpan = remainderBuffer.AsSpan(0, RemainderTableSize);
+        Span<byte> stepSpan = stepBuffer.AsSpan(0, RemainderTableSize);
+        // The remainder and step tables follow the fixed modulus order: 10, 8, 5, 3, 7, 11.
+
+        ulong startDivisor64 = divisor128.Low;
+        remainderSpan[0] = (byte)(startDivisor64 % 10UL);
+        remainderSpan[1] = (byte)(startDivisor64 % 8UL);
+        remainderSpan[2] = (byte)(startDivisor64 % 5UL);
+        remainderSpan[3] = (byte)(startDivisor64 % 3UL);
+        remainderSpan[4] = (byte)(startDivisor64 % 7UL);
+        remainderSpan[5] = (byte)(startDivisor64 % 11UL);
+
+        stepSpan[0] = (byte)(step % 10UL);
+        stepSpan[1] = (byte)(step % 8UL);
+        stepSpan[2] = (byte)(step % 5UL);
+        stepSpan[3] = (byte)(step % 3UL);
+        stepSpan[4] = (byte)(step % 7UL);
+        stepSpan[5] = (byte)(step % 11UL);
+
+        bool lastIsSeven = (prime & 3UL) == 3UL;
+
+        MersenneNumberDivisorRemainderGpuStepper remainderStepper = new(stepSpan);
+        MersenneNumberDivisorMontgomeryGpuBuilder? montgomeryBuilder = _montgomeryDevice == ByDivisorMontgomeryDevice.Gpu
+            ? new MersenneNumberDivisorMontgomeryGpuBuilder(batchCapacity)
+            : null;
+        MontgomeryDivisorData[]? montgomeryBuffer = montgomeryBuilder != null
+            ? ArrayPool<MontgomeryDivisorData>.Shared.Rent(batchCapacity)
+            : null;
+
+        ulong remainingCandidates = ComputeRemainingCandidateCount(startDivisor64, limit64, step);
+
+        try
+        {
+            while (remainingCandidates > 0UL && divisor128.CompareTo(limit128) <= 0)
+            {
+                if (divisor128.High != 0UL)
+                {
+                    break;
+                }
+
+                ulong currentDivisor = divisor128.Low;
+                int chunkCount = remainingCandidates > (ulong)batchCapacity
+                    ? batchCapacity
+                    : (int)remainingCandidates;
+                if (chunkCount <= 0)
+                {
+                    break;
+                }
+
+                Span<ulong> candidateSpan = candidateBuffer.AsSpan(0, chunkCount);
+                Span<byte> maskSpan = maskBuffer.AsSpan(0, chunkCount);
+
+                candidateEvaluator.EvaluateCandidates(
+                    currentDivisor,
+                    step,
+                    limit64,
+                    remainderSpan,
+                    stepSpan,
+                    lastIsSeven,
+                    candidateSpan,
+                    maskSpan);
+
+                Span<MontgomeryDivisorData> montgomerySpan = default;
+                if (montgomeryBuilder != null && montgomeryBuffer != null)
+                {
+                    montgomerySpan = montgomeryBuffer.AsSpan(0, chunkCount);
+                    // Every candidate divisor is unique across the scan, so build Montgomery data
+                    // for the current batch and drop it immediately instead of maintaining a cache.
+                    montgomeryBuilder.Build(candidateSpan, montgomerySpan);
+                }
+
+                for (int i = 0; i < chunkCount; i++)
+                {
+                    if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired())
+                    {
+                        processedAll = false;
+                        return false;
+                    }
+
+                    iterationCounter++;
+
+                    ulong candidate = candidateSpan[i];
+                    processedCount++;
+                    lastProcessed = candidate;
+
+                    if (maskSpan[i] == 0)
+                    {
+                        continue;
+                    }
+
+                    MontgomeryDivisorData divisorData;
+                    if (!montgomerySpan.IsEmpty)
+                    {
+                        divisorData = montgomerySpan[i];
+                        if (divisorData.Modulus != candidate)
+                        {
+                            // The GPU builder might skip entries when masks zero them out; recompute on the CPU
+                            // without caching because divisors never repeat within a session.
+                            divisorData = MontgomeryDivisorData.FromModulus(candidate);
+                        }
+                    }
+                    else
+                    {
+                        // Candidates remain unique even when GPU Montgomery support is disabled, so rebuild
+                        // the divisor data inline instead of persisting unused cache entries.
+                        divisorData = MontgomeryDivisorData.FromModulus(candidate);
+                    }
+
+                    ulong divisorCycle;
+                    if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
+                    {
+                        divisorCycle = cycleCache.GetCycleLength(candidate);
+                    }
+                    else
+                    {
+                        factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
+                        if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL)
+                        {
+                            divisorCycle = cycleCache.GetCycleLength(candidate);
+                        }
+                    }
+
+                    if (divisorCycle == prime)
+                    {
+                        processedAll = true;
+                        return true;
+                    }
+
+                    if (divisorCycle == 0UL)
+                    {
+                        Console.WriteLine($"Divisor cycle was not calculated for {prime}");
+                    }
+                }
+
+                GpuUInt128 increment = step128;
+                increment.Mul((ulong)chunkCount);
+                divisor128.Add(increment);
+
+                remainderStepper.Advance(chunkCount, remainderSpan);
+                remainingCandidates -= (ulong)chunkCount;
+            }
+        }
+        finally
+        {
+            montgomeryBuilder?.Dispose();
+            if (montgomeryBuffer != null)
+            {
+                ArrayPool<MontgomeryDivisorData>.Shared.Return(montgomeryBuffer, clearArray: false);
+            }
+
+            remainderStepper.Dispose();
+            ArrayPool<byte>.Shared.Return(remainderBuffer, clearArray: false);
+            ArrayPool<byte>.Shared.Return(stepBuffer, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(candidateBuffer, clearArray: false);
+            ArrayPool<byte>.Shared.Return(maskBuffer, clearArray: false);
+            candidateEvaluator.Dispose();
+        }
+
+        processedAll = remainingCandidates == 0UL && divisor128.CompareTo(limit128) > 0;
+        return false;
+    }
+
+    private static ulong ComputeRemainingCandidateCount(ulong currentDivisor, ulong limit, ulong step)
+    {
+        if (step == 0UL || currentDivisor > limit)
+        {
+            return 0UL;
+        }
+
+        ulong distance = limit - currentDivisor;
+        ulong stepsAvailable = distance / step;
+        return stepsAvailable + 1UL;
+    }
+
     private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
     {
-        ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData);
+        ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
         return residue == 1UL ? (byte)1 : (byte)0;
     }
 
@@ -298,16 +611,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         }
 
         ulong total = _lastStatusDivisor + processedCount;
-        if (total >= interval)
-        {
-            do
-            {
-                total -= interval;
-            }
-            while (total >= interval);
-        }
-
-        _lastStatusDivisor = total;
+        // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
+        // status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
+        _lastStatusDivisor = total % interval;
     }
 
     private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
@@ -344,6 +650,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
     {
         private readonly MersenneNumberDivisorByDivisorCpuTester _owner;
         private bool _disposed;
+        private ulong[]? _primeDeltas;
+        private ulong[]? _cycleRemainders;
+        private ulong[]? _residues;
+        private int _bufferCapacity;
+        private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator;
 
         internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner)
         {
@@ -368,10 +679,19 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 return;
             }
 
-            MontgomeryDivisorData cachedData = divisorData;
-            if (cachedData.Modulus != divisor)
+            MontgomeryDivisorData effectiveDivisorData = divisorData;
+            if (effectiveDivisorData.Modulus != divisor)
             {
-                cachedData = MontgomeryDivisorDataCache.Get(divisor);
+                // Each divisor surfaces only once across the monotonic scan, so rebuild the Montgomery data
+                // instead of retaining a cache entry that would never be reused within future checks.
+                effectiveDivisorData = MontgomeryDivisorData.FromModulus(divisor);
+            }
+
+            ulong modulus = effectiveDivisorData.Modulus;
+            if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+            {
+                hits.Clear();
+                return;
             }
 
             if (divisorCycle == 0UL)
@@ -384,9 +704,24 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 }
             }
 
+            if (_owner._deltasDevice == ByDivisorDeltasDevice.Cpu)
+            {
+                CheckDivisorCpu(effectiveDivisorData, divisorCycle, primes, hits);
+                return;
+            }
+
+            CheckDivisorGpu(effectiveDivisorData, divisorCycle, primes, hits, length);
+        }
+
+        private static void CheckDivisorCpu(
+            in MontgomeryDivisorData divisorData,
+            ulong divisorCycle,
+            ReadOnlySpan<ulong> primes,
+            Span<byte> hits)
+        {
             // 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);
+            var exponentStepper = new ExponentRemainderStepper(divisorData);
             if (!exponentStepper.IsValidModulus)
             {
                 hits.Clear();
@@ -400,6 +735,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0)
                 : (byte)0;
 
+            int length = primes.Length;
             for (int i = 1; i < length; i++)
             {
                 remainder = cycleStepper.ComputeNext(primes[i]);
@@ -413,6 +749,171 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             }
         }
 
+        private void CheckDivisorGpu(
+            in MontgomeryDivisorData divisorData,
+            ulong divisorCycle,
+            ReadOnlySpan<ulong> primes,
+            Span<byte> hits,
+            int length)
+        {
+            EnsureCapacity(length);
+
+            Span<ulong> deltaSpan = _primeDeltas!.AsSpan(0, length);
+            ComputePrimeDeltas(primes, deltaSpan);
+
+            Span<ulong> remainderSpan = _cycleRemainders!.AsSpan(0, length);
+            ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan);
+
+            EnsureGpuEvaluator();
+            Span<ulong> residueSpan = _residues!.AsSpan(0, length);
+            _gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan);
+
+            for (int i = 0; i < length; i++)
+            {
+                hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == 1UL ? (byte)1 : (byte)0;
+            }
+        }
+
+        private void EnsureGpuEvaluator()
+        {
+            if (_gpuResidueEvaluator is null)
+            {
+                _gpuResidueEvaluator = new MersenneNumberDivisorResidueGpuEvaluator(_owner._batchSize);
+            }
+        }
+
+        private void ReleaseGpuEvaluator()
+        {
+            if (_gpuResidueEvaluator is null)
+            {
+                return;
+            }
+
+            _gpuResidueEvaluator.Dispose();
+            _gpuResidueEvaluator = null;
+        }
+
+        private void EnsureCapacity(int requiredLength)
+        {
+            if (requiredLength <= 0)
+            {
+                return;
+            }
+
+            if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null)
+            {
+                return;
+            }
+
+            int newCapacity = _bufferCapacity;
+            if (newCapacity == 0)
+            {
+                newCapacity = 1;
+            }
+
+            while (newCapacity < requiredLength)
+            {
+                newCapacity <<= 1;
+            }
+
+            ReturnBuffers();
+
+            _primeDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
+            _cycleRemainders = ArrayPool<ulong>.Shared.Rent(newCapacity);
+            _residues = ArrayPool<ulong>.Shared.Rent(newCapacity);
+            _bufferCapacity = newCapacity;
+        }
+
+        private void ReturnBuffers()
+        {
+            if (_primeDeltas is not null)
+            {
+                ArrayPool<ulong>.Shared.Return(_primeDeltas, clearArray: false);
+                _primeDeltas = null;
+            }
+
+            if (_cycleRemainders is not null)
+            {
+                ArrayPool<ulong>.Shared.Return(_cycleRemainders, clearArray: false);
+                _cycleRemainders = null;
+            }
+
+            if (_residues is not null)
+            {
+                ArrayPool<ulong>.Shared.Return(_residues, clearArray: false);
+                _residues = null;
+            }
+
+            _bufferCapacity = 0;
+        }
+
+        private static void ComputePrimeDeltas(ReadOnlySpan<ulong> primes, Span<ulong> deltas)
+        {
+            if (primes.Length == 0)
+            {
+                return;
+            }
+
+            ulong previous = primes[0];
+            deltas[0] = previous;
+
+            for (int i = 1; i < primes.Length; i++)
+            {
+                ulong current = primes[i];
+                if (current <= previous)
+                {
+                    throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing.");
+                }
+
+                deltas[i] = current - previous;
+                previous = current;
+            }
+        }
+
+        private static void ComputeCycleRemainders(ulong cycleLength, ReadOnlySpan<ulong> primes, ReadOnlySpan<ulong> deltas, Span<ulong> remainders)
+        {
+            if (cycleLength == 0UL || primes.Length == 0)
+            {
+                remainders.Clear();
+                return;
+            }
+
+            ulong remainder = primes[0] % cycleLength;
+            remainders[0] = remainder;
+
+            for (int i = 1; i < primes.Length; i++)
+            {
+                remainder = AddMod(remainder, deltas[i], cycleLength);
+                remainders[i] = remainder;
+            }
+        }
+
+        private static ulong AddMod(ulong value, ulong delta, ulong modulus)
+        {
+            if (modulus == 0UL)
+            {
+                return 0UL;
+            }
+
+            if (ulong.MaxValue - value < delta)
+            {
+                UInt128 extended = (UInt128)value + delta;
+                return (ulong)(extended % modulus);
+            }
+
+            ulong sum = value + delta;
+            if (sum >= modulus)
+            {
+                sum -= modulus;
+                if (sum >= modulus)
+                {
+                    sum %= modulus;
+                }
+            }
+
+            return sum;
+        }
+
         public void Dispose()
         {
             if (_disposed)
@@ -421,6 +922,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             }
 
             _disposed = true;
+            ReleaseGpuEvaluator();
             _owner.ReturnSession(this);
         }
     }

PerfectNumbers.Core/DeterministicRandom.cs

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

PerfectNumbers.Core/DivisorCycleCache.cs

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

PerfectNumbers.Core/ExponentRemainderStepper.cs

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

PerfectNumbers.Core/FactorEntry.cs

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

PerfectNumbers.Core/Gpu/GpuContextPool.cs

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

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

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

PerfectNumbers.Core/Gpu/GpuUInt128.cs

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

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

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

PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs

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

PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs

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

PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs

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

PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs

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

PerfectNumbers.Core/Gpu/MontgomeryOddPowerGpu.cs

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

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -415,9 +415,9 @@ internal static partial class PrimeOrderGpuHeuristics
         ulong? previousOrder,
         PrimeOrderCalculator.PrimeOrderSearchConfig config,
         in MontgomeryDivisorData divisorData,
-        out PrimeOrderCalculator.PrimeOrderResult result)
+        out ulong order)
     {
-        result = default;
+        order = 0UL;
 
         var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
         var execution = lease.EnterExecutionScope();
@@ -497,18 +497,17 @@ internal static partial class PrimeOrderGpuHeuristics
             throw new InvalidOperationException("GPU Pollard Rho stack overflow; increase HeuristicStackCapacity.");
         }
 
-        ulong order = 0UL;
         resultBuffer.View.CopyToCPU(ref order, 1);
 
-        PrimeOrderCalculator.PrimeOrderStatus finalStatus = kernelStatus == PrimeOrderKernelStatus.Found
-            ? PrimeOrderCalculator.PrimeOrderStatus.Found
-            : PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved;
+        if (kernelStatus == PrimeOrderKernelStatus.FactoringFailure)
+        {
+            order = 0UL;
+        }
 
-        result = new PrimeOrderCalculator.PrimeOrderResult(finalStatus, order);
         DisposeResources();
         lease.Dispose();
 
-        return true;
+        return order != 0UL;
 
         void DisposeResources()
         {
@@ -812,7 +811,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             ulong reduced = phi / factor;
-            if (reduced.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL)
+            if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL)
             {
                 return false;
             }
@@ -1596,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
     {
-        return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL;
+        return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL;
     }
 
     private static ulong CalculateByDoublingKernel(ulong prime)
@@ -1716,10 +1715,10 @@ internal static partial class PrimeOrderGpuHeuristics
         GpuUInt128[]? rentedExponents = null;
         GpuUInt128[]? rentedResults = null;
         var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
+        var execution = lease.EnterExecutionScope();
 
         try
         {
-            var execution = lease.EnterExecutionScope();
             Accelerator accelerator = lease.Accelerator;
             AcceleratorStream stream = lease.Stream;
             var kernel = GetPow2ModWideKernel(accelerator);
@@ -1756,7 +1755,6 @@ internal static partial class PrimeOrderGpuHeuristics
 
             exponentBuffer.Dispose();
             remainderBuffer.Dispose();
-            execution.Dispose();
             return true;
         }
         catch (Exception)
@@ -1776,6 +1774,7 @@ internal static partial class PrimeOrderGpuHeuristics
                 ArrayPool<GpuUInt128>.Shared.Return(rentedResults, clearArray: false);
             }
 
+            execution.Dispose();
             lease.Dispose();
         }
     }
@@ -1796,7 +1795,7 @@ internal static partial class PrimeOrderGpuHeuristics
             return 0UL;
         }
 
-        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false);
+        return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
@@ -1812,7 +1811,57 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong exponent = exponents[index];
-        remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
+        ulong modulus = divisor.Modulus;
+		// ulong nPrime = divisor.NPrime;
+		// ulong result = divisor.MontgomeryOne;
+		// ulong baseValue = divisor.MontgomeryTwo;
+		// ulong remainingExponent = exponent;
+
+		// if (remainingExponent == 0UL)
+		// {
+		//     remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
+		//     return;
+		// }
+
+		// while (remainingExponent != 0UL)
+		// {
+		//     if ((remainingExponent & 1UL) != 0UL)
+		//     {
+		//         result = result.MontgomeryMultiply(baseValue, modulus, nPrime);
+		//     }
+
+		//     remainingExponent >>= 1;
+		//     if (remainingExponent == 0UL)
+		//     {
+		//         break;
+		//     }
+
+		//     // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term.
+		//     baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime);
+		// }
+
+		// remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
+
+		// Plain modular exponentiation without Montgomery reduction kept for reference.
+		ulong plainResult = 1UL;// % modulus;
+		ulong plainBase = 2UL;// % modulus;
+		while (exponent != 0UL)
+		{
+			if ((exponent & 1UL) != 0UL)
+			{
+				plainResult = (plainResult * plainBase) % modulus;
+			}
+
+			exponent >>= 1;
+			if (exponent == 0UL)
+			{
+				break;
+			}
+
+			plainBase = (plainBase * plainBase) % modulus;
+		}
+		
+        remainders[index] = plainResult;
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)

PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs

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

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Buffers;
 using System.Buffers.Binary;
 using System.Collections.Generic;
@@ -18,6 +19,12 @@ public class MersenneDivisorCycles
     // Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown
     private ulong[]? _smallCycles;
 
+    [ThreadStatic]
+    private static Dictionary<ulong, int>? s_factorCountsPool;
+
+    [ThreadStatic]
+    private static Dictionary<ulong, int>? s_factorScratchPool;
+
     public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles();
 
     public MersenneDivisorCycles()
@@ -102,7 +109,7 @@ public class MersenneDivisorCycles
             return cycleLength == exponent;
         }
 
-        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
     }
 
     public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
@@ -126,7 +133,7 @@ public class MersenneDivisorCycles
             }
         }
 
-        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
     }
 
     private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
@@ -234,14 +241,14 @@ public class MersenneDivisorCycles
             return (UInt128)cycle;
         }
 
-        PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
+        UInt128 wideOrder = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
             PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-        if (wideResult.Order != UInt128.Zero)
+        if (wideOrder != UInt128.Zero)
         {
-            return wideResult.Order;
+            return wideOrder;
         }
 
         // Otherwise, find order of 2 mod divisor
@@ -273,62 +280,85 @@ public class MersenneDivisorCycles
     {
         cycleLength = 0UL;
 
-        if ((divisor & (divisor - 1UL)) == 0UL)
-        {
-            cycleLength = 1UL;
-            return true;
-        }
+        // The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
+        // so this power-of-two fast-path never triggers on that call path.
+        // The historical fast-path is left here for reference but commented out because it is unreachable
+        // on the EvenPerfectBitScanner pipeline.
+        // if ((divisor & (divisor - 1UL)) == 0UL)
+        // {
+        //     cycleLength = 1UL;
+        //     return true;
+        // }
 
-        if (divisor <= 3UL || exponent <= 1UL)
-        {
-            cycleLength = CalculateCycleLength(divisor, divisorData);
-            return true;
-        }
+        // EvenPerfectBitScanner filters exponents <= 3 before invoking this helper and each
+        // divisor candidate is at least 2 * exponent + 1, so this guard never matches.
+        // The fallback calculation is therefore unreachable on the active call path.
+        // if (divisor <= 3UL || exponent <= 1UL)
+        // {
+        //     cycleLength = CalculateCycleLength(divisor, divisorData);
+        //     return true;
+        // }
 
-        PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
+        ulong computedOrder = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
             divisorData,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
             PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-        if (orderResult.Order != 0UL)
+        if (computedOrder != 0UL)
         {
-            cycleLength = orderResult.Order;
+            cycleLength = computedOrder;
             return true;
         }
 
         ulong phi = divisor - 1UL;
-        if (phi == 0UL)
-        {
-            return false;
-        }
+        // Divisors reach this helper only when they exceed the cached table, which means
+        // divisor >= 2 * exponent + 1 and phi cannot underflow to zero here.
+        // The defensive underflow guard is kept for posterity but disabled for the EvenPerfectBitScanner path.
+        // if (phi == 0UL)
+        // {
+        //     return false;
+        // }
 
         int twoCount = BitOperations.TrailingZeroCount(phi);
         ulong reducedPhi = phi >> twoCount;
-        if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
-        {
-            return false;
-        }
+        // Every divisor candidate generated by the by-divisor scan satisfies q ≡ 1 (mod 2 * exponent),
+        // ensuring reducedPhi keeps a factor of exponent and remains non-zero in this context.
+        // The divisibility check is therefore unreachable for the current caller and is commented out.
+        // if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
+        // {
+        //     return false;
+        // }
 
         ulong k = reducedPhi / exponent;
-        if (k == 0UL)
-        {
-            return false;
-        }
+        // reducedPhi always contains at least one factor of exponent, so dividing by exponent leaves k >= 1.
+        // if (k == 0UL)
+        // {
+        //     return false;
+        // }
 
-        Dictionary<ulong, int> factorCounts = new(8);
-        if (twoCount > 0)
+        Dictionary<ulong, int> factorCounts = AcquireFactorCountDictionary(8);
+        try
         {
+            // The divisor candidates are odd, so phi is even and this branch always records at least one factor of two.
+            // if (twoCount > 0)
+            // {
+            //     factorCounts[2UL] = twoCount;
+            // }
             factorCounts[2UL] = twoCount;
-        }
 
-        if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
+            if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
+            {
+                return false;
+            }
+
+            cycleLength = ReduceOrder(divisorData, phi, factorCounts, exponent);
+            return true;
+        }
+        finally
         {
-            return false;
+            ReleaseFactorCountDictionary(factorCounts);
         }
-
-        cycleLength = ReduceOrder(divisorData, phi, factorCounts);
-        return true;
     }
 
     private static bool AccumulateFactors(
@@ -342,61 +372,119 @@ public class MersenneDivisorCycles
             return true;
         }
 
-        if (!TryGetFactorization(value, cache, cacheResult, out FactorCacheEntry factorization))
+        FactorCacheEntry cachedFactorization;
+        PrimeFactorPartialResult? lease = null;
+        bool fromCache;
+
+        try
         {
-            return false;
+            if (!TryGetFactorization(value, cache, cacheResult, out cachedFactorization, out lease, out fromCache))
+            {
+                return false;
+            }
+
+            if (fromCache)
+            {
+                ulong[] primes = cachedFactorization.Primes;
+                byte[] exponents = cachedFactorization.Exponents;
+                int length = cachedFactorization.Count;
+
+                for (int i = 0; i < length; i++)
+                {
+                    AddFactor(counts, primes[i], exponents[i]);
+                }
+
+                return true;
+            }
+
+            if (lease is null || lease.Count == 0)
+            {
+                return true;
+            }
+
+            ReadOnlySpan<ulong> primesSpan = lease.Primes;
+            ReadOnlySpan<byte> exponentsSpan = lease.Exponents;
+            int spanLength = lease.Count;
+
+            for (int i = 0; i < spanLength; i++)
+            {
+                AddFactor(counts, primesSpan[i], exponentsSpan[i]);
+            }
+
+            return true;
         }
-
-        ulong[] primes = factorization.Primes;
-        byte[] exponents = factorization.Exponents;
-        int length = factorization.Count;
-
-        for (int i = 0; i < length; i++)
+        finally
         {
-            AddFactor(counts, primes[i], exponents[i]);
+            lease?.Dispose();
         }
-
-        return true;
     }
 
     private static bool TryGetFactorization(
         ulong value,
         Dictionary<ulong, FactorCacheEntry>? cache,
         bool cacheResult,
-        out FactorCacheEntry factorization)
+        out FactorCacheEntry factorization,
+        out PrimeFactorPartialResult? lease,
+        out bool fromCache)
     {
-        if (cache is not null && cache.TryGetValue(value, out factorization))
+        factorization = default;
+        lease = null;
+        fromCache = false;
+
+        if (value <= 1UL)
         {
             return true;
         }
 
-        Dictionary<ulong, int> scratch = new(8);
-        if (!TryFactorIntoCountsInternal(value, scratch))
+        if (cache is not null && cache.TryGetValue(value, out factorization))
         {
-            factorization = default!;
-            return false;
+            fromCache = true;
+            return true;
         }
 
-        int count = scratch.Count;
-        ulong[] primes = new ulong[count];
-        byte[] exponents = new byte[count];
-        int index = 0;
-        foreach (KeyValuePair<ulong, int> entry in scratch)
+        Dictionary<ulong, int> scratch = AcquireFactorScratchDictionary(8);
+        try
         {
-            primes[index] = entry.Key;
-            exponents[index] = checked((byte)entry.Value);
-            index++;
+            if (!TryFactorIntoCountsInternal(value, scratch))
+            {
+                factorization = default;
+                return false;
+            }
+
+            int count = scratch.Count;
+            if (count == 0)
+            {
+                return true;
+            }
+
+            PrimeFactorPartialResult result = PrimeFactorPartialResult.Rent(count);
+            lease = result;
+
+            Span<ulong> primes = result.PrimeWriteSpan;
+            Span<byte> exponents = result.ExponentWriteSpan;
+            int index = 0;
+            foreach (KeyValuePair<ulong, int> entry in scratch)
+            {
+                primes[index] = entry.Key;
+                exponents[index] = checked((byte)entry.Value);
+                index++;
+            }
+
+            result.CommitCount(index);
+            result.Sort();
+
+            if (cache is not null && cacheResult)
+            {
+                FactorCacheEntry cacheEntry = result.ToCacheEntry();
+                cache[value] = cacheEntry;
+            }
+
+            return true;
         }
-
-        Array.Sort(primes, exponents);
-        factorization = new FactorCacheEntry(primes, exponents, count);
-
-        if (cache is not null && cacheResult)
+        finally
         {
-            cache[value] = factorization;
+            ReleaseFactorScratchDictionary(scratch);
         }
-
-        return true;
     }
 
     private static bool TryFactorIntoCountsInternal(ulong value, Dictionary<ulong, int> counts)
@@ -471,15 +559,22 @@ public class MersenneDivisorCycles
         }
     }
 
-    private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
+    private static ulong ReduceOrder(
+        in MontgomeryDivisorData divisorData,
+        ulong initialOrder,
+        Dictionary<ulong, int> factorCounts,
+        ulong targetOrder)
     {
-        if (factorCounts.Count == 0)
+        if (factorCounts.Count == 0 || initialOrder == targetOrder)
         {
             return initialOrder;
         }
 
         int count = factorCounts.Count;
-        ulong[] primes = ArrayPool<ulong>.Shared.Rent(count);
+        bool rentFromPool = count >= PerfectNumberConstants.PooledArrayThreshold;
+        ulong[] primes = rentFromPool
+            ? ThreadLocalArrayPool<ulong>.Shared.Rent(count)
+            : new ulong[count];
         try
         {
             int index = 0;
@@ -505,9 +600,16 @@ public class MersenneDivisorCycles
                     }
 
                     ulong candidate = order / prime;
-                    if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+                    if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL)
                     {
                         order = candidate;
+                        if (order == targetOrder)
+                        {
+                            // Early-exit once the order matches the tested exponent;
+                            // further factor processing would not change the outcome for the by-divisor scan.
+                            return order;
+                        }
+
                         continue;
                     }
 
@@ -519,7 +621,10 @@ public class MersenneDivisorCycles
         }
         finally
         {
-            ArrayPool<ulong>.Shared.Return(primes, clearArray: false);
+            if (rentFromPool)
+            {
+                ThreadLocalArrayPool<ulong>.Shared.Return(primes, clearArray: false);
+            }
         }
     }
 
@@ -592,6 +697,53 @@ public class MersenneDivisorCycles
         }
     }
 
+    // Thread-local dictionary pools keep the factoring helpers reusable so the GC stays out of the by-divisor CPU loop.
+    private static Dictionary<ulong, int> AcquireFactorCountDictionary(int capacityHint)
+    {
+        Dictionary<ulong, int>? dictionary = s_factorCountsPool;
+        if (dictionary is null)
+        {
+            return new Dictionary<ulong, int>(capacityHint);
+        }
+
+        s_factorCountsPool = null;
+        dictionary.Clear();
+        dictionary.EnsureCapacity(capacityHint);
+        return dictionary;
+    }
+
+    private static void ReleaseFactorCountDictionary(Dictionary<ulong, int> dictionary)
+    {
+        dictionary.Clear();
+        if (s_factorCountsPool is null)
+        {
+            s_factorCountsPool = dictionary;
+        }
+    }
+
+    private static Dictionary<ulong, int> AcquireFactorScratchDictionary(int capacityHint)
+    {
+        Dictionary<ulong, int>? dictionary = s_factorScratchPool;
+        if (dictionary is null)
+        {
+            return new Dictionary<ulong, int>(capacityHint);
+        }
+
+        s_factorScratchPool = null;
+        dictionary.Clear();
+        dictionary.EnsureCapacity(capacityHint);
+        return dictionary;
+    }
+
+    private static void ReleaseFactorScratchDictionary(Dictionary<ulong, int> dictionary)
+    {
+        dictionary.Clear();
+        if (s_factorScratchPool is null)
+        {
+            s_factorScratchPool = dictionary;
+        }
+    }
+
     private static ulong NextRandomUInt64()
     {
         Span<byte> buffer = stackalloc byte[8];
@@ -658,7 +810,11 @@ public class MersenneDivisorCycles
             tasks[taskIndex] = Task.Run(() =>
             {
                 ulong rangeLength = localEnd - localStart + 1UL;
-                (ulong divisor, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent(checked((int)rangeLength));
+                int localLength = checked((int)rangeLength);
+                bool rentLocalCycles = localLength >= PerfectNumberConstants.PooledArrayThreshold;
+                (ulong divisor, ulong cycle)[] localCycles = rentLocalCycles
+                    ? ThreadLocalArrayPool<(ulong, ulong)>.Shared.Rent(localLength)
+                    : new (ulong divisor, ulong cycle)[localLength];
                 int localCycleIndex = 0;
 
                 try
@@ -706,7 +862,10 @@ public class MersenneDivisorCycles
                 }
                 finally
                 {
-                    ArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false);
+                    if (rentLocalCycles)
+                    {
+                        ThreadLocalArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false);
+                    }
                 }
             });
         }
@@ -721,7 +880,7 @@ public class MersenneDivisorCycles
         ulong start = 3UL;
         var locker = new object();
 
-        var pool = ArrayPool<ulong>.Shared;
+        var pool = ThreadLocalArrayPool<ulong>.Shared;
         ulong batchSizeUL = (ulong)batchSize, d, end;
         int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx;
         ulong[] divisors, outCycles, validDivisors;
@@ -926,29 +1085,34 @@ public class MersenneDivisorCycles
 
     internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength)
     {
-        if ((divisor & (divisor - 1UL)) == 0UL)
-        {
-            cycleLength = 1UL;
-            return true;
-        }
+        // The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
+        // so this power-of-two fast-path never triggers on that call path.
+        // The early exit is retained for reference but commented out because it is unreachable on the scanner path.
+        // if ((divisor & (divisor - 1UL)) == 0UL)
+        // {
+        //     cycleLength = 1UL;
+        //     return true;
+        // }
 
-        if (divisor <= 3UL)
-        {
-            cycleLength = CalculateCycleLengthFallback(divisor);
-            return true;
-        }
+        // By the divisor-by-divisor pipeline, candidates have the shape 2 * exponent * k + 1 with exponent >= 5,
+        // so the divisor never drops below 11 on this call path.
+        // if (divisor <= 3UL)
+        // {
+        //     cycleLength = CalculateCycleLengthFallback(divisor);
+        //     return true;
+        // }
 
         if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
         {
-            PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
+            ulong computedOrder = PrimeOrderCalculator.Calculate(
                     divisor,
                     previousOrder: null,
                     divisorData,
                     PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
                     PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-            if (orderResult.Order != 0UL)
+            if (computedOrder != 0UL)
             {
-                cycleLength = orderResult.Order;
+                cycleLength = computedOrder;
                 return true;
             }
         }

PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs

@@ -18,7 +18,8 @@ public static class MersenneNumberDivisorByDivisorTester
             Action markComposite,
             Action clearComposite,
             Action<ulong, bool, bool, bool> printResult,
-            int threadCount)
+            int threadCount,
+            TimeSpan? primeTestLimit)
     {
         if (candidates.Count == 0)
         {
@@ -112,9 +113,12 @@ public static class MersenneNumberDivisorByDivisorTester
         ulong maxPrime = 0UL;
         int primeWriteIndex = 0;
 
-        using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
-        bool hasPrime = primeEnumerator.MoveNext();
-        ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
+		// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
+		// We don't need to additionally check it.
+
+        // using IEnumerator<ulong> primeEnumerator = Prime.Numbers.GetEnumerator();
+        // bool hasPrime = primeEnumerator.MoveNext();
+        // ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL;
 
         for (int index = startIndex; index < candidateCount; index++)
         {
@@ -127,25 +131,29 @@ public static class MersenneNumberDivisorByDivisorTester
                 continue;
             }
 
-            while (hasPrime && currentPrime < candidate)
-            {
-                hasPrime = primeEnumerator.MoveNext();
-                if (!hasPrime)
-                {
-                    break;
-                }
+			// This implementation is terribly slow, while this method expect prime p given as --filter-p input already.
+			// We don't need to additionally check it.
+			
+			// while (hasPrime && currentPrime < candidate)
+			// {
+			//     hasPrime = primeEnumerator.MoveNext();
+			//     if (!hasPrime)
+			//     {
+			//         break;
+			//     }
 
-                currentPrime = primeEnumerator.Current;
-            }
+			//     currentPrime = primeEnumerator.Current;
+			// }
 
-            if (!hasPrime || currentPrime != candidate)
-            {
-                markComposite();
-                printResult(candidate, false, false, false);
-                continue;
-            }
+			// if (!hasPrime || currentPrime != candidate)
+			// if (!Number.IsPrime(candidate))
+			// {
+			//     markComposite();
+			//     printResult(candidate, false, false, false);
+			//     continue;
+			// }
 
-            primesSpan[primeWriteIndex++] = candidate;
+			primesSpan[primeWriteIndex++] = candidate;
 
             if (candidate > maxPrime)
             {
@@ -236,7 +244,7 @@ public static class MersenneNumberDivisorByDivisorTester
 
         void ProcessPrime(ulong prime)
         {
-            bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted);
+            bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted, primeTestLimit);
 
             if (!isPrime)
             {
@@ -249,21 +257,52 @@ public static class MersenneNumberDivisorByDivisorTester
             printResult(prime, true, divisorsExhausted, true);
         }
 
-        if (workerCount == 1)
+        if (workerCount > filteredPrimes.Count)
+        {
+            workerCount = filteredPrimes.Count;
+        }
+
+        if (workerCount <= 1)
         {
             foreach (ulong prime in filteredPrimes)
             {
                 ProcessPrime(prime);
             }
-        }
-        else
-        {
-            ParallelOptions options = new()
-            {
-                MaxDegreeOfParallelism = workerCount
-            };
 
-            Parallel.ForEach(filteredPrimes, options, ProcessPrime);
+            return;
         }
+
+        (int Start, int End)[] workerRanges = new (int, int)[workerCount];
+        int totalPrimeCount = filteredPrimes.Count;
+        int baseGroupSize = totalPrimeCount / workerCount;
+        int remaining = totalPrimeCount % workerCount;
+        int groupStartIndex = 0;
+
+        for (int workerIndex = 0; workerIndex < workerCount; workerIndex++)
+        {
+            int groupSize = baseGroupSize;
+            if (workerIndex < remaining)
+            {
+                groupSize++;
+            }
+
+            startIndex = groupStartIndex;
+            int endIndex = startIndex + groupSize;
+            workerRanges[workerIndex] = (startIndex, endIndex);
+            groupStartIndex = endIndex;
+        }
+
+        ParallelOptions parallelOptions = new()
+        {
+            MaxDegreeOfParallelism = workerCount
+        };
+
+        Parallel.ForEach(workerRanges, parallelOptions, range =>
+        {
+            for (int candidateIndex = range.Start; candidateIndex < range.End; candidateIndex++)
+            {
+                ProcessPrime(filteredPrimes[candidateIndex]);
+            }
+        });
     }
 }

PerfectNumbers.Core/MersennePrimeFactorTester.cs

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

PerfectNumbers.Core/MontgomeryDivisorData.cs

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

PerfectNumbers.Core/MutableUInt128.cs

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

PerfectNumbers.Core/PartialFactorResult.cs

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

PerfectNumbers.Core/PerfectNumberConstants.cs

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

PerfectNumbers.Core/PerfectNumbers.Core.csproj

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

PerfectNumbers.Core/PrimeFactorPartialResult.cs

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

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -9,30 +9,15 @@ namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
-    internal readonly struct PrimeOrderResultWide
-    {
-        public PrimeOrderResultWide(PrimeOrderStatus status, UInt128 order)
-        {
-            Status = status;
-            Order = order;
-        }
-
-        public PrimeOrderStatus Status { get; }
-
-        public UInt128 Order { get; }
-    }
-
-    public static PrimeOrderResultWide Calculate(UInt128 prime, UInt128? previousOrder, PrimeOrderSearchConfig config)
-        => Calculate(prime, previousOrder, config, PrimeOrderHeuristicDevice.Gpu);
-
-    public static PrimeOrderResultWide Calculate(
+    public static UInt128 Calculate(
         UInt128 prime,
         UInt128? previousOrder,
         PrimeOrderSearchConfig config,
         PrimeOrderHeuristicDevice device)
     {
-        using var scope = UsePow2Mode(device);
+        var scope = UsePow2Mode(device);
         MontgomeryDivisorData divisorData;
+        UInt128 result;
         if (prime <= ulong.MaxValue)
         {
             ulong? previous = null;
@@ -51,41 +36,49 @@ internal static partial class PrimeOrderCalculator
 
             ulong prime64 = (ulong)prime;
             divisorData = MontgomeryDivisorData.FromModulus(prime64);
-            PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device);
-            return new PrimeOrderResultWide(result.Status, result.Order);
+            ulong order64 = Calculate(prime64, previous, divisorData, config, device);
+            result = order64 == 0UL ? UInt128.Zero : (UInt128)order64;
         }
         else
         {
             divisorData = default;
+            result = CalculateWideInternal(prime, previousOrder, divisorData, config);
         }
 
-        return CalculateWideInternal(prime, previousOrder, divisorData, config);
+        scope.Dispose();
+        return result;
     }
 
-    private static PrimeOrderResultWide CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+    private static UInt128 CalculateWideInternal(UInt128 prime, UInt128? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
     {
         if (prime <= UInt128.One)
         {
-            return new PrimeOrderResultWide(PrimeOrderStatus.Found, UInt128.One);
+            return UInt128.One;
         }
 
         if (prime == (UInt128)3UL)
         {
-            return new PrimeOrderResultWide(PrimeOrderStatus.Found, (UInt128)2UL);
+            return (UInt128)2UL;
         }
 
         UInt128 phi = prime - UInt128.One;
-        using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
+        var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
         PartialFactorResult128 phiFactors = PartialFactorWide(phi, config);
+        UInt128 result;
         if (phiFactors.Factors is null)
         {
-            return FinishStrictlyWide(prime, divisorData, config.Mode);
+            result = FinishStrictlyWide(prime, divisorData);
+        }
+        else
+        {
+            result = RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
         }
 
-        return RunHeuristicPipelineWide(prime, previousOrder, divisorData, config, phi, phiFactors);
+        debugScope.Dispose();
+        return result;
     }
 
-    private static PrimeOrderResultWide RunHeuristicPipelineWide(
+    private static UInt128 RunHeuristicPipelineWide(
         UInt128 prime,
         UInt128? previousOrder,
         in MontgomeryDivisorData divisorData,
@@ -95,7 +88,7 @@ internal static partial class PrimeOrderCalculator
     {
         if (phiFactors.FullyFactored && TrySpecialMaxWide(phi, prime, divisorData, phiFactors))
         {
-            return new PrimeOrderResultWide(PrimeOrderStatus.Found, phi);
+            return phi;
         }
 
         UInt128 candidateOrder = InitializeStartingOrderWide(prime, phi, divisorData);
@@ -103,27 +96,25 @@ internal static partial class PrimeOrderCalculator
 
         if (TryConfirmOrderWide(prime, candidateOrder, divisorData, config))
         {
-            return new PrimeOrderResultWide(PrimeOrderStatus.Found, candidateOrder);
+            return candidateOrder;
         }
 
         if (config.Mode == PrimeOrderMode.Strict)
         {
-            return FinishStrictlyWide(prime, divisorData, PrimeOrderMode.Strict);
+            return FinishStrictlyWide(prime, divisorData);
         }
 
         if (TryHeuristicFinishWide(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order))
         {
-            return new PrimeOrderResultWide(PrimeOrderStatus.Found, order);
+            return order;
         }
 
-        return FinishStrictlyWide(prime, divisorData, config.Mode);
+        return candidateOrder;
     }
 
-    private static PrimeOrderResultWide FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData, PrimeOrderMode mode)
+    private static UInt128 FinishStrictlyWide(UInt128 prime, in MontgomeryDivisorData divisorData)
     {
-        UInt128 strictOrder = CalculateByFactorizationWide(prime, divisorData);
-        PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved;
-        return new PrimeOrderResultWide(status, strictOrder);
+        return CalculateByFactorizationWide(prime, divisorData);
     }
 
     private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)
@@ -171,7 +162,7 @@ internal static partial class PrimeOrderCalculator
                 return order;
             }
 
-            tempArray = ArrayPool<FactorEntry128>.Shared.Rent(length + 1);
+            tempArray = ThreadLocalArrayPool<FactorEntry128>.Shared.Rent(length + 1);
             Span<FactorEntry128> buffer = tempArray.AsSpan(0, length);
             factorSpan.CopyTo(buffer);
 
@@ -209,7 +200,7 @@ internal static partial class PrimeOrderCalculator
         {
             if (tempArray is not null)
             {
-                ArrayPool<FactorEntry128>.Shared.Return(tempArray, clearArray: false);
+                ThreadLocalArrayPool<FactorEntry128>.Shared.Return(tempArray, clearArray: false);
             }
         }
     }
@@ -673,7 +664,7 @@ internal static partial class PrimeOrderCalculator
             return new PartialFactorResult128(null, value, false, 0);
         }
 
-        FactorEntry128[]? rented = counts.Count > 0 ? ArrayPool<FactorEntry128>.Shared.Rent(counts.Count) : null;
+        FactorEntry128[]? rented = counts.Count > 0 ? ThreadLocalArrayPool<FactorEntry128>.Shared.Rent(counts.Count) : null;
         int index = 0;
         if (rented is not null)
         {
@@ -691,7 +682,7 @@ internal static partial class PrimeOrderCalculator
             span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
             resultArray = new FactorEntry128[index];
             span.CopyTo(resultArray);
-            ArrayPool<FactorEntry128>.Shared.Return(rented, clearArray: false);
+            ThreadLocalArrayPool<FactorEntry128>.Shared.Return(rented, clearArray: false);
         }
 
         bool fullyFactored = cofactor == UInt128.One;
@@ -1071,15 +1062,20 @@ internal static partial class PrimeOrderCalculator
 
     private readonly struct FactorEntry128
     {
-        public FactorEntry128(UInt128 value, int exponent)
-        {
-            Value = value;
-            Exponent = exponent;
-        }
+		public FactorEntry128(UInt128 value)
+		{
+			Value = value;
+		}
 
-        public UInt128 Value { get; }
+		public FactorEntry128(UInt128 value, int exponent)
+		{
+			Value = value;
+			Exponent = exponent;
+		}
 
-        public int Exponent { get; }
+        public readonly UInt128 Value;
+
+		public readonly int Exponent = 1;
     }
 
     private readonly struct PartialFactorResult128
@@ -1092,13 +1088,13 @@ internal static partial class PrimeOrderCalculator
             Count = count;
         }
 
-        public FactorEntry128[]? Factors { get; }
+		public readonly FactorEntry128[]? Factors;
 
-        public UInt128 Cofactor { get; }
+		public readonly UInt128 Cofactor;
 
-        public bool FullyFactored { get; }
+		public readonly bool FullyFactored;
 
-        public int Count { get; }
+		public readonly int Count;
 
         public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0);
 
@@ -1106,14 +1102,13 @@ internal static partial class PrimeOrderCalculator
         {
             if (Factors is null)
             {
-                FactorEntry128[] local = new FactorEntry128[1];
-                local[0] = new FactorEntry128(prime, 1);
-                return new PartialFactorResult128(local, UInt128.One, true, 1);
+                FactorEntry128[] local = [new FactorEntry128(prime)];
+				return new PartialFactorResult128(local, UInt128.One, true, 1);
             }
 
             FactorEntry128[] extended = new FactorEntry128[Count + 1];
             Array.Copy(Factors, extended, Count);
-            extended[Count] = new FactorEntry128(prime, 1);
+            extended[Count] = new FactorEntry128(prime);
             return new PartialFactorResult128(extended, UInt128.One, true, Count + 1);
         }
     }

PerfectNumbers.Core/PrimeTestTimeLimit.cs

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

PerfectNumbers.Core/ThreadLocalArrayPool.cs

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

PerfectNumbers.Core/UInt128Extensions.cs

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

PerfectNumbers.Core/ULongExtensions.cs

@@ -1,4 +1,5 @@
 using System;
+using System.Buffers;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using ILGPU.Algorithms;
@@ -76,7 +77,7 @@ public static class ULongExtensions
     public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift)
     {
         ulong high = value >> shift;
-        ulong mask = (1UL << shift) - 1UL;
+		ulong mask = (1UL << shift) - 1UL;
         ulong low = value & mask;
 
         UInt128 highContribution = (UInt128)high * multiplier;
@@ -309,6 +310,37 @@ public static class ULongExtensions
         return state.MulModWithNativeModulo(b, modulus);
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus)
+    {
+        if (modulus <= 1UL)
+        {
+            return 0UL;
+        }
+
+        ulong remainingExponent = exponent;
+        ulong result = 1UL % modulus;
+        ulong baseValue = 2UL % modulus;
+
+        while (remainingExponent != 0UL)
+        {
+            if ((remainingExponent & 1UL) != 0UL)
+            {
+                result = result.MulMod64GpuCompatible(baseValue, modulus);
+            }
+
+            remainingExponent >>= 1;
+            if (remainingExponent == 0UL)
+            {
+                break;
+            }
+
+            baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus);
+        }
+
+        return result;
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime)
     {
@@ -326,7 +358,7 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+    public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
     {
         ulong modulus = divisor.Modulus;
         if (exponent == 0UL)
@@ -345,7 +377,93 @@ public static class ULongExtensions
 
         ulong result = divisor.MontgomeryOne;
         ulong nPrime = divisor.NPrime;
-        ulong[] oddPowers = InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowerCount);
+
+        Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : default;
+        ulong[]? pooledOddPowers = null;
+
+		if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold)
+        {
+            oddPowerStorage = new ulong[oddPowerCount];
+        }
+        else
+        {
+            pooledOddPowers = ArrayPool<ulong>.Shared.Rent(oddPowerCount);
+            oddPowerStorage = pooledOddPowers;
+        }
+
+        Span<ulong> oddPowers = oddPowerStorage[..oddPowerCount];
+        InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowers);
+
+        int index = bitLength - 1;
+        while (index >= 0)
+        {
+            if (((exponent >> index) & 1UL) == 0UL)
+            {
+                result = result.MontgomeryMultiply(result, modulus, nPrime);
+                index--;
+                continue;
+            }
+
+            int windowStart = index - windowSize + 1;
+            if (windowStart < 0)
+            {
+                windowStart = 0;
+            }
+
+            while (((exponent >> windowStart) & 1UL) == 0UL)
+            {
+                windowStart++;
+            }
+
+            int windowLength = index - windowStart + 1;
+            for (int square = 0; square < windowLength; square++)
+            {
+                result = result.MontgomeryMultiply(result, modulus, nPrime);
+            }
+
+            ulong mask = (1UL << windowLength) - 1UL;
+            ulong windowValue = (exponent >> windowStart) & mask;
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+            ulong multiplier = oddPowers[tableIndex];
+            result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
+
+            index = windowStart - 1;
+        }
+
+        if (pooledOddPowers is not null)
+        {
+            ArrayPool<ulong>.Shared.Return(pooledOddPowers);
+        }
+
+        if (keepMontgomery)
+        {
+            return result;
+        }
+
+        return result.MontgomeryMultiply(1UL, modulus, nPrime);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+    {
+        ulong modulus = divisor.Modulus;
+        if (exponent == 0UL)
+        {
+            return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
+        }
+
+        if (exponent <= Pow2WindowFallbackThreshold)
+        {
+            return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
+        }
+
+        int bitLength = GetPortableBitLength(exponent);
+        int windowSize = GetWindowSize(bitLength);
+        int oddPowerCount = 1 << (windowSize - 1);
+
+        ulong result = divisor.MontgomeryOne;
+        ulong nPrime = divisor.NPrime;
+        ulong[] oddPowers = InitializeMontgomeryOddPowersGpu(divisor, modulus, nPrime, oddPowerCount);
 
         int index = bitLength - 1;
         while (index >= 0)
@@ -461,7 +579,38 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong[] InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
+    private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> destination)
+    {
+        int oddPowerCount = destination.Length;
+        if (oddPowerCount == 0)
+        {
+            return;
+        }
+
+        bool computedOnGpu = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
+            && MontgomeryOddPowerGpu.TryCompute(divisor, oddPowerCount, destination);
+        if (computedOnGpu)
+        {
+            return;
+        }
+
+        destination[0] = divisor.MontgomeryTwo;
+        if (oddPowerCount == 1)
+        {
+            return;
+        }
+
+        ulong square = divisor.MontgomeryTwoSquared;
+        ulong current = destination[0];
+        for (int i = 1; i < oddPowerCount; i++)
+        {
+            current = current.MontgomeryMultiply(square, modulus, nPrime);
+            destination[i] = current;
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong[] InitializeMontgomeryOddPowersGpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
     {
         ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount];
         oddPowers[0] = divisor.MontgomeryTwo;
@@ -481,17 +630,31 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
+    public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
     {
         ulong rotationCount = exponent % cycleLength;
-        return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false);
+        return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+    public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
     {
-        return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false);
+        return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false);
     }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
+    {
+        ulong rotationCount = exponent % cycleLength;
+        return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+    {
+        return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false);
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
     {

gpu-montgomery-usage.md

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

tools/HeuristicMultiplyShiftTiming/HeuristicMultiplyShiftTiming.csproj

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