TODO added

Remove default prime order overloads

EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs

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

PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs

@@ -15,9 +15,15 @@ public class GpuContextPoolTests
     {
         GpuContextPool.DisposeAll();
 
-        using var lease = GpuContextPool.RentPreferred(preferCpu: true);
+        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,11 +199,18 @@ 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);
+        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,7 +150,9 @@ public class NttGpuMathTests
         var modulus = new GpuUInt128(0UL, 17UL);
         var primitiveRoot = new GpuUInt128(0UL, 3UL);
 
-        using var gpu = GpuContextPool.Rent();
+        GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
+        try
+        {
             var accelerator = gpu.Accelerator;
 
             var method = typeof(NttGpuMath).GetMethod("GetSquareCache", BindingFlags.NonPublic | BindingFlags.Static);
@@ -175,6 +177,11 @@ public class NttGpuMathTests
                 prod.Should().Be(r2);
             }
         }
+        finally
+        {
+            gpu.Dispose();
+        }
+    }
 
     [Fact]
     [Trait("Category", "Fast")]
@@ -338,13 +345,24 @@ 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);
+            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.
         var field = typeof(NttGpuMath).GetField("SquareCache", BindingFlags.NonPublic | BindingFlags.Static)!;

PerfectNumbers.Core.Tests/Gpu/NttMontgomeryTests.cs

@@ -33,7 +33,9 @@ public class NttMontgomeryTests
         NttGpuMath.Square(cpu, modulus, primitiveRoot);
 
         // Device execution
-        using var context = GpuContextPool.Rent();
+        GpuContextPool.GpuContextLease context = GpuContextPool.Rent();
+        try
+        {
             var accelerator = context.Accelerator;
             using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
             buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
@@ -45,6 +47,11 @@ public class NttMontgomeryTests
                 gpu[i].Should().Be(cpu[i]);
             }
         }
+        finally
+        {
+            context.Dispose();
+        }
+    }
 
     [Fact]
     [Trait("Category", "Fast")]

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.Should().Be(3UL);
-            result.Order.Should().Be(3UL);
         }
         finally
         {

PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs

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

PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs

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

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -238,20 +238,19 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                     }
                 }
 
-				if (divisorCycle != 0 && divisorCycle == prime)
+                if (divisorCycle == prime)
                 {
+                    // A cycle equal to the tested exponent (which is prime in this path)
+                    // guarantees that the candidate divides the corresponding Mersenne
+                    // number because the order of 2 modulo the divisor is exactly p.
                     processedAll = true;
                     return true;
                 }
-				else if (divisorCycle == 0)
+
+                if (divisorCycle == 0UL)
                 {
-					Console.WriteLine($"Divisor cycle wasn't calculated for ${prime}");
+                    Console.WriteLine($"Divisor cycle was not calculated for {prime}");
                 }
-				// if (divisorCycle != 0UL && CheckDivisor(prime, divisorCycle, divisorData) != 0)
-				// {
-				// 	processedAll = true;
-				// 	return true;
-				// }
             }
 
             divisor += step;
@@ -269,7 +268,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
     private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
     {
-        ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData);
+        ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
         return residue == 1UL ? (byte)1 : (byte)0;
     }
 

PerfectNumbers.Core/DeterministicRandom.cs

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

PerfectNumbers.Core/ExponentRemainderStepper.cs

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

PerfectNumbers.Core/FactorEntry.cs

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

PerfectNumbers.Core/Gpu/GpuContextPool.cs

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

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

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

PerfectNumbers.Core/Gpu/GpuUInt128.cs

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

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -551,7 +551,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulong exponent = exponents[index];
-        hits[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -770,7 +770,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         ulong exponent = exponents[index];
 
-        results[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: true);
+        results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
     }
 
     public ulong DivisorLimit

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

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

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -102,7 +102,7 @@ public class MersenneDivisorCycles
             return cycleLength == exponent;
         }
 
-        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
     }
 
     public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
@@ -126,7 +126,7 @@ public class MersenneDivisorCycles
             }
         }
 
-        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
     }
 
     private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
@@ -234,14 +234,14 @@ public class MersenneDivisorCycles
             return (UInt128)cycle;
         }
 
-        PrimeOrderCalculator.PrimeOrderResultWide wideResult = PrimeOrderCalculator.Calculate(
+        UInt128 wideOrder = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
             PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-        if (wideResult.Order != UInt128.Zero)
+        if (wideOrder != UInt128.Zero)
         {
-            return wideResult.Order;
+            return wideOrder;
         }
 
         // Otherwise, find order of 2 mod divisor
@@ -285,15 +285,15 @@ public class MersenneDivisorCycles
             return true;
         }
 
-        PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
+        ulong computedOrder = PrimeOrderCalculator.Calculate(
             divisor,
             previousOrder: null,
             divisorData,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
             PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-        if (orderResult.Order != 0UL)
+        if (computedOrder != 0UL)
         {
-            cycleLength = orderResult.Order;
+            cycleLength = computedOrder;
             return true;
         }
 
@@ -505,7 +505,7 @@ public class MersenneDivisorCycles
                     }
 
                     ulong candidate = order / prime;
-                    if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+                    if (candidate.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL)
                     {
                         order = candidate;
                         continue;
@@ -940,15 +940,15 @@ public class MersenneDivisorCycles
 
         if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
         {
-            PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
+            ulong computedOrder = PrimeOrderCalculator.Calculate(
                     divisor,
                     previousOrder: null,
                     divisorData,
                     PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
                     PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
-            if (orderResult.Order != 0UL)
+            if (computedOrder != 0UL)
             {
-                cycleLength = orderResult.Order;
+                cycleLength = computedOrder;
                 return true;
             }
         }

PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs

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

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);
+            c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
-            x = (ulong)rng.Next(2, int.MaxValue);
+            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);
+            c = (DeterministicRandom.NextUInt128() % (n - 1)) + 1;
-            x = (UInt128)rng.NextInt64(2, long.MaxValue);
+            x = (DeterministicRandom.NextUInt128() % (n - 2)) + 2;
             y = x;
             d = 1;
 

PerfectNumbers.Core/PartialFactorResult.cs

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

PerfectNumbers.Core/PerfectNumberConstants.cs

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

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -9,30 +9,15 @@ namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
-    internal readonly struct PrimeOrderResultWide
+    public static UInt128 Calculate(
-    {
-        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(
         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);
+            ulong order64 = Calculate(prime64, previous, divisorData, config, device);
-            return new PrimeOrderResultWide(result.Status, result.Order);
+            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);
+        return CalculateByFactorizationWide(prime, divisorData);
-        PrimeOrderStatus status = mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved;
-        return new PrimeOrderResultWide(status, strictOrder);
     }
 
     private static bool TrySpecialMaxWide(UInt128 phi, UInt128 prime, in MontgomeryDivisorData divisorData, PartialFactorResult128 factors)

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -1,65 +1,30 @@
-using System;
 using System.Buffers;
 using System.Buffers.Binary;
-using System.Collections.Generic;
 using System.Diagnostics;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
-using System.Security.Cryptography;
 using PerfectNumbers.Core.Gpu;
-using System.Threading;
 
 namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
-    internal enum PrimeOrderStatus
-    {
-        Found,
-        HeuristicUnresolved,
-    }
-
-    internal readonly struct PrimeOrderResult
-    {
-        public PrimeOrderResult(PrimeOrderStatus status, ulong order)
-        {
-            Status = status;
-            Order = order;
-        }
-
-        public PrimeOrderStatus Status { get; }
-
-        public ulong Order { get; }
-    }
-
     internal enum PrimeOrderMode
     {
         Heuristic,
         Strict,
     }
 
-    internal readonly struct PrimeOrderSearchConfig
+    internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode)
     {
-        public PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode)
+        public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 10_240, maxPowChecks: 24, PrimeOrderMode.Heuristic);
-        {
+        public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict);
-            SmallFactorLimit = smallFactorLimit;
-            PollardRhoMilliseconds = pollardRhoMilliseconds;
-            MaxPowChecks = maxPowChecks;
-            Mode = mode;
-        }
 
-        public uint SmallFactorLimit { get; }
+        public readonly uint SmallFactorLimit = smallFactorLimit;
-
+        public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds;
-        public int PollardRhoMilliseconds { get; }
+        public readonly int MaxPowChecks = maxPowChecks;
-
+        public readonly PrimeOrderMode Mode = mode;
-        public int MaxPowChecks { get; }
-
-        public PrimeOrderMode Mode { get; }
-
-        public static PrimeOrderSearchConfig HeuristicDefault => new(100_000, 128, 24, PrimeOrderMode.Heuristic);
-
-        public static PrimeOrderSearchConfig StrictDefault => new(1_000_000, 0, 0, PrimeOrderMode.Strict);
     }
 
     internal enum PrimeOrderHeuristicDevice
@@ -76,15 +41,24 @@ internal static partial class PrimeOrderCalculator
 
     private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
 
-
     [ThreadStatic]
     private static PrimeOrderHeuristicDevice s_deviceMode;
 
     [ThreadStatic]
     private static bool s_debugLoggingEnabled;
 
+    [ThreadStatic]
+    private static List<ulong>? s_partialFactorPending;
+
+    [ThreadStatic]
+    private static Stack<ulong>? s_partialFactorCompositeStack;
+
+    [ThreadStatic]
+    private static Dictionary<ulong, int>? s_partialFactorCounts;
+
     private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu;
-    private readonly struct Pow2ModeScope : IDisposable
+
+    private readonly struct Pow2ModeScope
     {
         private readonly bool _previousInitialized;
         private readonly bool _previousAllow;
@@ -108,7 +82,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private readonly struct DebugLoggingScope : IDisposable
+    private readonly struct DebugLoggingScope
     {
         private readonly bool _previous;
 
@@ -146,93 +120,89 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+    public static ulong Calculate(
-        => Calculate(prime, previousOrder, divisorData, config, PrimeOrderHeuristicDevice.Gpu);
-
-    public static PrimeOrderResult Calculate(
         ulong prime,
         ulong? previousOrder,
         in MontgomeryDivisorData divisorData,
-        PrimeOrderSearchConfig config,
+        in PrimeOrderSearchConfig config,
         PrimeOrderHeuristicDevice device)
     {
-        using var scope = UsePow2Mode(device);
+        var scope = UsePow2Mode(device);
-        return CalculateInternal(prime, previousOrder, divisorData, config);
+        ulong order = CalculateInternal(prime, previousOrder, divisorData, config);
+        scope.Dispose();
+        return order;
     }
 
-    private static PrimeOrderResult CalculateInternal(ulong prime, ulong? previousOrder, MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+    private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
 	{
+		// TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31.
         if (prime <= 3UL)
         {
-            return new PrimeOrderResult(PrimeOrderStatus.Found, prime == 3UL ? 2UL : 1UL);
+            return prime == 3UL ? 2UL : 1UL;
         }
 
         ulong phi = prime - 1UL;
 
-        if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out PrimeOrderResult gpuResult))
+        if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder))
         {
-            return gpuResult;
+            return gpuOrder;
         }
 
-        using var debugScope = UseDebugLogging(!IsGpuHeuristicDevice);
-        DebugLog("Partial factoring φ(p)");
-
         PartialFactorResult phiFactors = PartialFactor(phi, config);
 
+        ulong result;
         if (phiFactors.Factors is null)
         {
-            DebugLog("No factors found");
+            result = FinishStrictlyCpu(prime, divisorData);
-            return FinishStrictly(prime, config.Mode);
+            phiFactors.Dispose();
+            return result;
         }
 
-        return RunHeuristicPipeline(prime, previousOrder, config, divisorData, phi, phiFactors);
+        result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors);
+        phiFactors.Dispose();
+        return result;
     }
 
-    private static PrimeOrderResult RunHeuristicPipeline(
+    private static ulong RunHeuristicPipelineCpu(
         ulong prime,
         ulong? previousOrder,
-        PrimeOrderSearchConfig config,
+        in PrimeOrderSearchConfig config,
         in MontgomeryDivisorData divisorData,
         ulong phi,
         PartialFactorResult phiFactors)
     {
-        DebugLog("Trying special max check");
+        if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData))
-        if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData))
         {
-            return new PrimeOrderResult(PrimeOrderStatus.Found, phi);
+            return phi;
         }
 
-        DebugLog("Initializing starting order");
+        ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData);
-        ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData);
+        candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData);
-        candidateOrder = ExponentLowering(candidateOrder, prime, phiFactors, divisorData);
 
-        DebugLog("Trying to confirm order");
+        if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config))
-        if (TryConfirmOrder(prime, candidateOrder, divisorData, config))
         {
-            return new PrimeOrderResult(PrimeOrderStatus.Found, candidateOrder);
+            return candidateOrder;
         }
 
         if (config.Mode == PrimeOrderMode.Strict)
         {
-            return FinishStrictly(prime, PrimeOrderMode.Strict);
+            return FinishStrictlyCpu(prime, divisorData);
         }
 
-        if (TryHeuristicFinish(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order))
+        if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order))
         {
-            return new PrimeOrderResult(PrimeOrderStatus.Found, order);
+            return order;
         }
 
-        DebugLog("Heuristic unresolved, finishing strictly");
+        return candidateOrder;
-        return FinishStrictly(prime, config.Mode);
     }
 
-    private static PrimeOrderResult FinishStrictly(ulong prime, PrimeOrderMode mode)
+    private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData)
     {
-        ulong strictOrder = CalculateByDoubling(prime);
+        return CalculateByFactorizationCpu(prime, divisorData);
-        return new PrimeOrderResult(mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved, strictOrder);
     }
 
-    private static bool TrySpecialMax(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
+    private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
     {
         ReadOnlySpan<FactorEntry> factorSpan = factors.Factors;
         int length = factors.Count;
@@ -240,7 +210,7 @@ internal static partial class PrimeOrderCalculator
         {
             ulong factor = factorSpan[i].Value;
             ulong reduced = phi / factor;
-            if (Pow2EqualsOne(reduced, prime, divisorData))
+            if (Pow2EqualsOneCpu(reduced, prime, divisorData))
             {
                 return false;
             }
@@ -249,13 +219,13 @@ internal static partial class PrimeOrderCalculator
         return true;
     }
 
-    private static ulong InitializeStartingOrder(ulong prime, ulong phi, in MontgomeryDivisorData divisorData)
+    private static ulong InitializeStartingOrderCpu(ulong prime, ulong phi, in MontgomeryDivisorData divisorData)
     {
         ulong order = phi;
         if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL)
         {
             ulong half = phi >> 1;
-            if (Pow2EqualsOne(half, prime, divisorData))
+            if (Pow2EqualsOneCpu(half, prime, divisorData))
             {
                 order = half;
             }
@@ -264,7 +234,7 @@ internal static partial class PrimeOrderCalculator
         return order;
     }
 
-    private static ulong ExponentLowering(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
+    private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
     {
         FactorEntry[]? tempArray = null;
         try
@@ -297,7 +267,7 @@ internal static partial class PrimeOrderCalculator
                     }
 
                     ulong reduced = order / primeFactor;
-                    if (Pow2EqualsOne(reduced, prime, divisorData))
+                    if (Pow2EqualsOneCpu(reduced, prime, divisorData))
                     {
                         order = reduced;
                         continue;
@@ -318,7 +288,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool TryConfirmOrder(ulong prime, ulong order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config)
+    private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
     {
         if (order == 0UL)
         {
@@ -326,16 +296,18 @@ internal static partial class PrimeOrderCalculator
         }
 
         // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`.
-        DebugLog("Verifying a^order ≡ 1 (mod p)");
+        // DebugLog("Verifying a^order ≡ 1 (mod p)");
-        if (!Pow2EqualsOne(order, prime, divisorData))
+        if (!Pow2EqualsOneCpu(order, prime, divisorData))
         {
             return false;
         }
 
-        DebugLog("Partial factoring order");
+        // DebugLog("Partial factoring order");
 
         // TODO: Do we do partial factoring of order multiple times?
         PartialFactorResult factorization = PartialFactor(order, config);
+        try
+        {
             if (factorization.Factors is null)
             {
                 return false;
@@ -345,12 +317,10 @@ internal static partial class PrimeOrderCalculator
             {
                 if (factorization.Cofactor <= 1UL)
                 {
-            DebugLog("Cofactor <= 1. No factors found");
                     return false;
                 }
 
-            // TODO: Use Open.Numerics.Primality for this final check once it's available.
+                // DebugLog("Cofactor > 1. Testing primality of cofactor");
-            DebugLog("Cofactor > 1. Testing primality of cofactor");
                 bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
                 // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
                 if (!isPrime)
@@ -358,12 +328,14 @@ internal static partial class PrimeOrderCalculator
                     return false;
                 }
 
-            DebugLog("Adding cofactor as prime factor");
+                // DebugLog("Adding cofactor as prime factor");
-            factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
+                PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
+                factorization.Dispose();
+                factorization = extended;
             }
 
             ReadOnlySpan<FactorEntry> span = factorization.Factors;
-        DebugLog("Verifying prime-power reductions");
+            // DebugLog("Verifying prime-power reductions");
             int length = factorization.Count;
             for (int i = 0; i < length; i++)
             {
@@ -377,7 +349,7 @@ internal static partial class PrimeOrderCalculator
                     }
 
                     reduced /= primeFactor;
-                if (Pow2EqualsOne(reduced, prime, divisorData))
+                    if (Pow2EqualsOneCpu(reduced, prime, divisorData))
                     {
                         return false;
                     }
@@ -386,13 +358,18 @@ internal static partial class PrimeOrderCalculator
 
             return true;
         }
+        finally
+        {
+            factorization.Dispose();
+        }
+    }
 
-    private static bool TryHeuristicFinish(
+    private static bool TryHeuristicFinishCpu(
         ulong prime,
         ulong order,
         ulong? previousOrder,
         in MontgomeryDivisorData divisorData,
-        PrimeOrderSearchConfig config,
+        in PrimeOrderSearchConfig config,
         PartialFactorResult phiFactors,
         out ulong result)
     {
@@ -403,8 +380,10 @@ internal static partial class PrimeOrderCalculator
         }
 
         // TODO: Do we do partial factoring of order multiple times?
-        DebugLog("Trying heuristic. Partial factoring order");
+        // DebugLog("Trying heuristic. Partial factoring order");
         PartialFactorResult orderFactors = PartialFactor(order, config);
+        try
+        {
             if (orderFactors.Factors is null)
             {
                 return false;
@@ -424,20 +403,22 @@ internal static partial class PrimeOrderCalculator
                     return false;
                 }
 
-            orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
+                PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
+                orderFactors.Dispose();
+                orderFactors = extended;
             }
 
             int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4;
             List<ulong> candidates = new(capacity);
             FactorEntry[] factorArray = orderFactors.Factors!;
-        DebugLog("Building candidates list");
+            // DebugLog("Building candidates list");
             BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity);
             if (candidates.Count == 0)
             {
                 return false;
             }
 
-        DebugLog("Sorting candidates");
+            // DebugLog("Sorting candidates");
             SortCandidates(prime, previousOrder, candidates);
 
             int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
@@ -446,10 +427,10 @@ internal static partial class PrimeOrderCalculator
             bool allowGpuBatch = true;
             Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
 
-        DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
+            // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
             int index = 0;
-        const int MaxGpuBatchSize = 1024;
+            const int MaxGpuBatchSize = 256;
-        const int StackGpuBatchSize = 1024;
+            const int StackGpuBatchSize = 64;
             Span<ulong> stackGpuRemainders = stackalloc ulong[StackGpuBatchSize];
             while (index < candidateCount && powUsed < powBudget)
             {
@@ -516,14 +497,14 @@ internal static partial class PrimeOrderCalculator
                     }
                     else
                     {
-                    equalsOne = Pow2EqualsOne(candidate, prime, divisorData);
+                        equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData);
                     }
                     if (!equalsOne)
                     {
                         continue;
                     }
 
-                if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget))
+                    if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget))
                     {
                         continue;
                     }
@@ -545,9 +526,14 @@ internal static partial class PrimeOrderCalculator
                 index += batchSize;
             }
 
-        DebugLog("No candidate confirmed");
+            // DebugLog("No candidate confirmed");
             return false;
         }
+        finally
+        {
+            orderFactors.Dispose();
+        }
+    }
 
     private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates)
     {
@@ -716,9 +702,11 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static bool TryConfirmCandidate(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, ref int powUsed, int powBudget)
+    private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget)
     {
         PartialFactorResult factorization = PartialFactor(candidate, config);
+        try
+        {
             if (factorization.Factors is null)
             {
                 return false;
@@ -738,7 +726,9 @@ internal static partial class PrimeOrderCalculator
                     return false;
                 }
 
-            factorization = factorization.WithAdditionalPrime(factorization.Cofactor);
+                PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
+                factorization.Dispose();
+                factorization = extended;
             }
 
             ReadOnlySpan<FactorEntry> span = factorization.Factors;
@@ -761,7 +751,7 @@ internal static partial class PrimeOrderCalculator
                     }
 
                     powUsed++;
-                if (Pow2EqualsOne(reduced, prime, divisorData))
+                    if (Pow2EqualsOneCpu(reduced, prime, divisorData))
                     {
                         return false;
                     }
@@ -770,8 +760,13 @@ internal static partial class PrimeOrderCalculator
 
             return true;
         }
+        finally
+        {
+            factorization.Dispose();
+        }
+    }
 
-    private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
+    private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
     {
         if (IsGpuPow2Allowed)
         {
@@ -782,10 +777,59 @@ internal static partial class PrimeOrderCalculator
             }
         }
 
-        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
     }
 
-    private static PartialFactorResult PartialFactor(ulong value, PrimeOrderSearchConfig config)
+    private static List<ulong> AcquirePendingCompositeList(int capacityHint)
+    {
+        List<ulong>? list = s_partialFactorPending;
+        if (list is null)
+        {
+            list = new List<ulong>(capacityHint);
+            s_partialFactorPending = list;
+            return list;
+        }
+
+        list.Clear();
+        if (list.Capacity < capacityHint)
+        {
+            list.Capacity = capacityHint;
+        }
+
+        return list;
+    }
+
+    private static Stack<ulong> AcquireCompositeStack(int capacityHint)
+    {
+        Stack<ulong>? stack = s_partialFactorCompositeStack;
+        if (stack is null)
+        {
+            stack = new Stack<ulong>(capacityHint);
+            s_partialFactorCompositeStack = stack;
+            return stack;
+        }
+
+        stack.Clear();
+        stack.EnsureCapacity(capacityHint);
+        return stack;
+    }
+
+    private static Dictionary<ulong, int> AcquireFactorCountDictionary(int capacityHint)
+    {
+        Dictionary<ulong, int>? dictionary = s_partialFactorCounts;
+        if (dictionary is null)
+        {
+            dictionary = new Dictionary<ulong, int>(capacityHint);
+            s_partialFactorCounts = dictionary;
+            return dictionary;
+        }
+
+        dictionary.Clear();
+        dictionary.EnsureCapacity(capacityHint);
+        return dictionary;
+    }
+
+    private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config)
     {
         if (value <= 1UL)
         {
@@ -803,18 +847,19 @@ internal static partial class PrimeOrderCalculator
         bool useDictionary = false;
 
         uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit;
-        ulong remaining;
+        ulong remaining = value;
 
         bool gpuFactored = false;
         if (IsGpuHeuristicDevice)
         {
             gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _);
         }
-        else
-        {
-            remaining = value;
-        }
 
+        List<ulong> pending = AcquirePendingCompositeList(2);
+        Stack<ulong>? compositeStack = null;
+
+        try
+        {
             if (!gpuFactored)
             {
                 primeSlots.Clear();
@@ -822,12 +867,11 @@ internal static partial class PrimeOrderCalculator
                 if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining))
                 {
                     useDictionary = true;
-                counts = new Dictionary<ulong, int>(capacity: 8);
+                    counts = AcquireFactorCountDictionary(Math.Max(factorCount, 8));
                     remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts);
                 }
             }
 
-        List<ulong> pending = new();
             if (remaining > 1UL)
             {
                 pending.Add(remaining);
@@ -835,16 +879,15 @@ internal static partial class PrimeOrderCalculator
 
             if (config.PollardRhoMilliseconds > 0 && pending.Count > 0)
             {
-            DebugLog("Processing pending composites with Pollard's Rho");
                 Stopwatch stopwatch = Stopwatch.StartNew();
                 long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks;
-            Stack<ulong> stack = new();
+                compositeStack = AcquireCompositeStack(Math.Max(pending.Count * 2, 4));
-            stack.Push(remaining);
+                compositeStack.Push(remaining);
                 pending.Clear();
 
-            while (stack.Count > 0)
+                while (compositeStack.Count > 0)
                 {
-                ulong composite = stack.Pop();
+                    ulong composite = compositeStack.Pop();
                     if (composite <= 1UL)
                     {
                         continue;
@@ -870,8 +913,8 @@ internal static partial class PrimeOrderCalculator
                     }
 
                     ulong quotient = composite / factor;
-                stack.Push(factor);
+                    compositeStack.Push(factor);
-                stack.Push(quotient);
+                    compositeStack.Push(quotient);
                 }
             }
 
@@ -895,18 +938,16 @@ internal static partial class PrimeOrderCalculator
             {
                 if ((counts is null || counts.Count == 0) && cofactor == value)
                 {
-                DebugLog("cofactor is the same as value, no factors found");
+                    return PartialFactorResult.Rent(null, value, false, 0);
-                return new PartialFactorResult(null, value, false, 0);
                 }
 
                 if (counts is null || counts.Count == 0)
                 {
-                return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
+                    return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
                 }
 
                 FactorEntry[] factors = ArrayPool<FactorEntry>.Shared.Rent(counts.Count);
                 int index = 0;
-            DebugLog(() => $"Collecting {counts.Count} prime factors");
                 foreach (KeyValuePair<ulong, int> entry in counts)
                 {
                     factors[index] = new FactorEntry(entry.Key, entry.Value);
@@ -916,12 +957,8 @@ internal static partial class PrimeOrderCalculator
                 Span<FactorEntry> span = factors.AsSpan(0, index);
                 span.Sort(static (a, b) => a.Value.CompareTo(b.Value));
 
-            FactorEntry[] resultArray = new FactorEntry[index];
-            span.CopyTo(resultArray);
-            ArrayPool<FactorEntry>.Shared.Return(factors, clearArray: false);
-
                 bool fullyFactored = cofactor == 1UL;
-            return new PartialFactorResult(resultArray, cofactor, fullyFactored, index);
+                return PartialFactorResult.Rent(factors, cofactor, fullyFactored, index);
             }
 
             int actualCount = 0;
@@ -937,16 +974,14 @@ internal static partial class PrimeOrderCalculator
             {
                 if (cofactor == value)
                 {
-                DebugLog("cofactor is the same as value, no factors found");
+                    return PartialFactorResult.Rent(null, value, false, 0);
-                return new PartialFactorResult(null, value, false, 0);
                 }
 
-            return new PartialFactorResult(null, cofactor, cofactor == 1UL, 0);
+                return PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0);
             }
 
-        FactorEntry[] array = new FactorEntry[actualCount];
+            FactorEntry[] array = ArrayPool<FactorEntry>.Shared.Rent(actualCount);
             int arrayIndex = 0;
-        DebugLog(() => $"Collecting {actualCount} prime factors");
             for (int i = 0; i < factorCount; i++)
             {
                 ulong primeValue = primeSlots[i];
@@ -960,9 +995,17 @@ internal static partial class PrimeOrderCalculator
                 arrayIndex++;
             }
 
-        Array.Sort(array, static (a, b) => a.Value.CompareTo(b.Value));
+            Span<FactorEntry> arraySpan = array.AsSpan(0, arrayIndex);
+            arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value));
             bool fullyFactoredArray = cofactor == 1UL;
-        return new PartialFactorResult(array, cofactor, fullyFactoredArray, array.Length);
+            return PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex);
+        }
+        finally
+        {
+            pending.Clear();
+            counts?.Clear();
+            compositeStack?.Clear();
+        }
     }
 
     private const int GpuSmallPrimeFactorSlots = 64;
@@ -1114,8 +1157,6 @@ internal static partial class PrimeOrderCalculator
             return true;
         }
 
-        Span<byte> buffer = stackalloc byte[8];
-
         while (true)
         {
             if (stopwatch.ElapsedTicks > budgetTicks)
@@ -1123,10 +1164,8 @@ internal static partial class PrimeOrderCalculator
                 return false;
             }
 
-            RandomNumberGenerator.Fill(buffer);
+            ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
-            ulong c = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 1UL)) + 1UL;
+            ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
-            RandomNumberGenerator.Fill(buffer);
-            ulong x = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 2UL)) + 2UL;
             ulong y = x;
             ulong d = 1UL;
 
@@ -1254,7 +1293,7 @@ internal static partial class PrimeOrderCalculator
             return;
         }
 
-        counts ??= new Dictionary<ulong, int>(Math.Max(count, 8));
+        counts ??= AcquireFactorCountDictionary(Math.Max(count, 8));
         CopyFactorsToDictionary(primes, exponents, count, counts);
         count = 0;
         useDictionary = true;
@@ -1280,59 +1319,98 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static ulong CalculateByDoubling(ulong prime)
+    private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData)
     {
-        ulong order = 1UL;
+        ulong phi = prime - 1UL;
-        ulong pow = 2UL;
+        Dictionary<ulong, int> counts = new(capacity: 8);
-
+        FactorCompletely(phi, counts);
-        while (pow != 1UL)
+        if (counts.Count == 0)
         {
-            pow <<= 1;
+            return phi;
-            if (pow >= prime)
-            {
-                pow -= prime;
         }
 
-            order++;
+        List<KeyValuePair<ulong, int>> entries = new(counts);
+        entries.Sort(static (a, b) => a.Key.CompareTo(b.Key));
+
+        ulong order = phi;
+        int entryCount = entries.Count;
+        for (int i = 0; i < entryCount; i++)
+        {
+            ulong primeFactor = entries[i].Key;
+            int exponent = entries[i].Value;
+            for (int iteration = 0; iteration < exponent; iteration++)
+            {
+                if ((order % primeFactor) != 0UL)
+                {
+                    break;
+                }
+
+                ulong candidate = order / primeFactor;
+                if (Pow2EqualsOneCpu(candidate, prime, divisorData))
+                {
+                    order = candidate;
+                    continue;
+                }
+
+                break;
+            }
         }
 
         return order;
     }
 
+    private static void FactorCompletely(ulong value, Dictionary<ulong, int> counts)
+    {
+        if (value <= 1UL)
+        {
+            return;
+        }
+
+        if (Open.Numeric.Primes.Prime.Numbers.IsPrime(value))
+        {
+            AddFactor(counts, value, 1);
+            return;
+        }
+
+        ulong factor = PollardRhoStrict(value);
+        FactorCompletely(factor, counts);
+        FactorCompletely(value / factor, counts);
+    }
+
+    private static ulong PollardRhoStrict(ulong n)
+    {
+        if ((n & 1UL) == 0UL)
+        {
+            return 2UL;
+        }
+
+        while (true)
+        {
+            ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL;
+            ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL;
+            ulong y = x;
+            ulong d = 1UL;
+
+            while (d == 1UL)
+            {
+                x = AdvancePolynomial(x, c, n);
+                y = AdvancePolynomial(y, c, n);
+                y = AdvancePolynomial(y, c, n);
+                ulong diff = x > y ? x - y : y - x;
+                d = BinaryGcd(diff, n);
+            }
+
+            if (d != n)
+            {
+                return d;
+            }
+        }
+    }
+
     private readonly struct CandidateKey(int primary, long secondary, long tertiary)
     {
         public readonly int Primary = primary;
         public readonly long Secondary = secondary;
         public readonly long Tertiary = tertiary;
     }
-
-    private readonly struct FactorEntry(ulong value, int exponent)
-    {
-        public readonly ulong Value = value;
-        public readonly int Exponent = exponent;
-    }
-
-    private readonly struct PartialFactorResult(FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count)
-    {
-        public readonly FactorEntry[]? Factors = factors;
-        public readonly ulong Cofactor = cofactor;
-        public readonly bool FullyFactored = fullyFactored;
-        public readonly int Count = count;
-        public readonly static PartialFactorResult Empty = new(null, 1UL, true, 0);
-
-        public PartialFactorResult WithAdditionalPrime(ulong prime)
-        {
-            if (Factors is null)
-            {
-                FactorEntry[] local = [new FactorEntry(prime, 1)];
-                return new PartialFactorResult(local, 1UL, true, 1);
-            }
-
-            FactorEntry[] extended = new FactorEntry[Count + 1];
-            Array.Copy(Factors, 0, extended, 0, Count);
-            extended[Count] = new FactorEntry(prime, 1);
-            Array.Sort(extended, static (a, b) => a.Value.CompareTo(b.Value));
-            return new PartialFactorResult(extended, 1UL, true, Count + 1);
-        }
-    }
 }

PerfectNumbers.Core/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;
+            exponent.ShiftRight(windowStart);
-            shifted.ShiftRight(windowStart);
+            return exponent.Low & mask;
-            ulong mask = (1UL << windowBitCount) - 1UL;
-            return shifted.Low & mask;
         }
 
-        ulong directMask = (1UL << windowBitCount) - 1UL;
+        return exponent.Low & mask;
-        return exponent.Low & directMask;
     }
 
     [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;
@@ -326,7 +327,7 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+    public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
     {
         ulong modulus = divisor.Modulus;
         if (exponent == 0UL)
@@ -345,7 +346,78 @@ public static class ULongExtensions
 
         ulong result = divisor.MontgomeryOne;
         ulong nPrime = divisor.NPrime;
-        ulong[] oddPowers = InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowerCount);
+        ulong[] oddPowers = InitializeMontgomeryOddPowersCpu(divisor, modulus, nPrime, oddPowerCount);
+
+        int index = bitLength - 1;
+		while (index >= 0)
+		{
+			if (((exponent >> index) & 1UL) == 0UL)
+			{
+				result = result.MontgomeryMultiply(result, modulus, nPrime);
+				index--;
+				continue;
+			}
+
+			int windowStart = index - windowSize + 1;
+			if (windowStart < 0)
+			{
+				windowStart = 0;
+			}
+
+			while (((exponent >> windowStart) & 1UL) == 0UL)
+			{
+				windowStart++;
+			}
+
+			int windowLength = index - windowStart + 1;
+			for (int square = 0; square < windowLength; square++)
+			{
+				result = result.MontgomeryMultiply(result, modulus, nPrime);
+			}
+
+			ulong mask = (1UL << windowLength) - 1UL;
+			ulong windowValue = (exponent >> windowStart) & mask;
+			int tableIndex = (int)((windowValue - 1UL) >> 1);
+			ulong multiplier = oddPowers[tableIndex];
+			result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
+
+			index = windowStart - 1;
+		}
+		
+		if (oddPowerCount >= PerfectNumberConstants.PooledArrayThreshold)
+		{
+			ArrayPool<ulong>.Shared.Return(oddPowers);
+		}
+
+        if (keepMontgomery)
+        {
+            return result;
+        }
+
+        return result.MontgomeryMultiply(1UL, modulus, nPrime);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+    {
+        ulong modulus = divisor.Modulus;
+        if (exponent == 0UL)
+        {
+            return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
+        }
+
+        if (exponent <= Pow2WindowFallbackThreshold)
+        {
+            return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
+        }
+
+        int bitLength = GetPortableBitLength(exponent);
+        int windowSize = GetWindowSize(bitLength);
+        int oddPowerCount = 1 << (windowSize - 1);
+
+        ulong result = divisor.MontgomeryOne;
+        ulong nPrime = divisor.NPrime;
+        ulong[] oddPowers = InitializeMontgomeryOddPowersGpu(divisor, modulus, nPrime, oddPowerCount);
 
         int index = bitLength - 1;
         while (index >= 0)
@@ -461,7 +533,30 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong[] InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
+    private static ulong[] InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
+    {
+		ulong[] oddPowers = oddPowerCount < PerfectNumberConstants.PooledArrayThreshold ? new ulong[oddPowerCount] : ArrayPool<ulong>.Shared.Rent(oddPowerCount);
+		
+        oddPowers[0] = divisor.MontgomeryTwo;
+        if (oddPowerCount == 1)
+        {
+            return oddPowers;
+        }
+
+        ulong previous,
+			  square = divisor.MontgomeryTwoSquared;
+
+        for (int i = 1; i < oddPowerCount; i++)
+        {
+            previous = oddPowers[i - 1];
+            oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime);
+        }
+
+        return oddPowers;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong[] InitializeMontgomeryOddPowersGpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, int oddPowerCount)
     {
         ulong[] oddPowers = new ulong[PerfectNumberConstants.MaxOddPowersCount];
         oddPowers[0] = divisor.MontgomeryTwo;
@@ -481,17 +576,31 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
+    public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
     {
         ulong rotationCount = exponent % cycleLength;
-        return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false);
+        return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+    public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
     {
-        return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false);
+        return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false);
     }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
+    {
+        ulong rotationCount = exponent % cycleLength;
+        return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor)
+    {
+        return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false);
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
     {