Replace status modulo with subtract loop

EvenPerfectBitScanner.Benchmarks/GpuPow2ModBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs

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

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModUInt128Benchmarks.cs

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

PerfectNumbers.Core.Tests/Gpu/GpuContextPoolTests.cs

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

PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs

@@ -199,18 +199,11 @@ public class GpuKernelCompilationTests
     [MemberData(nameof(KernelLoaders))]
     public void Kernel_compiles_without_internal_compiler_errors(string name, Action<Accelerator> compile)
     {
-        GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
+        using var 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,9 +150,7 @@ public class NttGpuMathTests
         var modulus = new GpuUInt128(0UL, 17UL);
         var primitiveRoot = new GpuUInt128(0UL, 3UL);
 
-        GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
+        using var gpu = GpuContextPool.Rent();
-        try
-        {
         var accelerator = gpu.Accelerator;
 
         var method = typeof(NttGpuMath).GetMethod("GetSquareCache", BindingFlags.NonPublic | BindingFlags.Static);
@@ -177,11 +175,6 @@ public class NttGpuMathTests
             prod.Should().Be(r2);
         }
     }
-        finally
-        {
-            gpu.Dispose();
-        }
-    }
 
     [Fact]
     [Trait("Category", "Fast")]
@@ -345,24 +338,13 @@ public class NttGpuMathTests
         values[1] = new GpuUInt128(0UL, 2UL);
         var modulus = new GpuUInt128(0UL, 17UL);
         var primitiveRoot = new GpuUInt128(0UL, 3UL);
-        GpuContextPool.GpuContextLease gpu = GpuContextPool.Rent();
+        using (var 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,9 +33,7 @@ public class NttMontgomeryTests
         NttGpuMath.Square(cpu, modulus, primitiveRoot);
 
         // Device execution
-        GpuContextPool.GpuContextLease context = GpuContextPool.Rent();
+        using var context = GpuContextPool.Rent();
-        try
-        {
         var accelerator = context.Accelerator;
         using var buffer = accelerator.Allocate1D<GpuUInt128>(gpu.Length);
         buffer.View.CopyFromCPU(ref gpu[0], gpu.Length);
@@ -47,11 +45,6 @@ public class NttMontgomeryTests
             gpu[i].Should().Be(cpu[i]);
         }
     }
-        finally
-        {
-            context.Dispose();
-        }
-    }
 
     [Fact]
     [Trait("Category", "Fast")]

PerfectNumbers.Core.Tests/Gpu/PrimeOrderGpuHeuristicsTests.cs

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

PerfectNumbers.Core.Tests/GpuKernelPoolTests.cs

@@ -11,9 +11,7 @@ public class GpuKernelPoolTests
     [Trait("Category", "Fast")]
     public void EnsureSmallCyclesOnDevice_returns_same_view_on_subsequent_calls()
     {
-        GpuContextPool.GpuContextLease lease = GpuContextPool.RentPreferred(preferCpu: true);
+        using var lease = GpuContextPool.RentPreferred(preferCpu: true);
-        try
-        {
         var accelerator = lease.Accelerator;
 
         var first = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
@@ -22,11 +20,6 @@ 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);
 
-        ulong result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Should().Be(expected);
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        result.Order.Should().Be(expected);
     }
 
     [Theory]
@@ -62,21 +62,21 @@ public class PrimeOrderCalculatorTests
     {
         ulong previousOrder = ComputeOrderByDoubling(previousPrime);
 
-        ulong heuristic = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        ulong strict = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.StrictDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        heuristic.Should().Be(strict);
+        heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        heuristic.Order.Should().Be(strict.Order);
     }
 
     [Theory]
@@ -85,14 +85,14 @@ public class PrimeOrderCalculatorTests
     [InlineData(3UL, 2UL)]
     public void Calculate_handles_trivial_primes(ulong prime, ulong expectedOrder)
     {
-        ulong result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(prime),
-            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Should().Be(expectedOrder);
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        result.Order.Should().Be(expectedOrder);
     }
 
     [Fact]
@@ -105,14 +105,14 @@ public class PrimeOrderCalculatorTests
             maxPowChecks: 8,
             mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
 
-        ulong result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult result = PrimeOrderCalculator.Calculate(
             prime: 13UL,
             previousOrder: null,
             MontgomeryDivisorData.FromModulus(13UL),
-            config,
+            config);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Should().Be(ComputeOrderByDoubling(13UL));
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
+        result.Order.Should().Be(ComputeOrderByDoubling(13UL));
     }
 
     [Fact]
@@ -126,14 +126,14 @@ public class PrimeOrderCalculatorTests
             mode: PrimeOrderCalculator.PrimeOrderMode.Heuristic);
 
         MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(239UL);
-        ulong heuristic = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult heuristic = PrimeOrderCalculator.Calculate(
             prime: 239UL,
             previousOrder: null,
             divisorData,
-            heuristicConfig,
+            heuristicConfig);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        heuristic.Should().Be(ComputeOrderByDoubling(239UL));
+        heuristic.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.HeuristicUnresolved);
+        heuristic.Order.Should().Be(ComputeOrderByDoubling(239UL));
 
         var strictConfig = new PrimeOrderCalculator.PrimeOrderSearchConfig(
             smallFactorLimit: 2,
@@ -141,14 +141,14 @@ public class PrimeOrderCalculatorTests
             maxPowChecks: 1,
             mode: PrimeOrderCalculator.PrimeOrderMode.Strict);
 
-        ulong strict = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResult strict = PrimeOrderCalculator.Calculate(
             prime: 239UL,
             previousOrder: null,
             divisorData,
-            strictConfig,
+            strictConfig);
-			      PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        strict.Should().Be(heuristic);
+        strict.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        strict.Order.Should().Be(heuristic.Order);
     }
 
     [Fact]
@@ -157,13 +157,13 @@ public class PrimeOrderCalculatorTests
     {
         UInt128 prime = UInt128.Parse("18446744073709641691");
 
-        UInt128 result = PrimeOrderCalculator.Calculate(
+        PrimeOrderCalculator.PrimeOrderResultWide result = PrimeOrderCalculator.Calculate(
             prime,
             previousOrder: null,
-                PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
-                PrimeOrderCalculator.PrimeOrderHeuristicDevice.Gpu);
 
-        result.Should().Be((UInt128)1229782938247309446UL);
+        result.Status.Should().Be(PrimeOrderCalculator.PrimeOrderStatus.Found);
+        result.Order.Should().Be((UInt128)1229782938247309446UL);
     }
 
     [Fact]

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -80,9 +80,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             return true;
         }
 
-        ulong processedCount;
+        ulong processedCount = 0UL;
-        ulong lastProcessed;
+        ulong lastProcessed = 0UL;
-        bool processedAll;
+        bool processedAll = false;
 
         bool composite = CheckDivisors(
             prime,
@@ -268,7 +268,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
     private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
     {
-        ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData);
+        ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData);
         return residue == 1UL ? (byte)1 : (byte)0;
     }
 
@@ -298,9 +298,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         }
 
         ulong total = _lastStatusDivisor + processedCount;
-        // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
+        if (total >= interval)
-        // status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
+        {
-        _lastStatusDivisor = total % interval;
+            do
+            {
+                total -= interval;
+            }
+            while (total >= interval);
+        }
+
+        _lastStatusDivisor = total;
     }
 
     private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)

PerfectNumbers.Core/DeterministicRandom.cs

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

PerfectNumbers.Core/ExponentRemainderStepper.cs

@@ -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.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_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.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+        ulong multiplier = delta.Pow2MontgomeryModWindowed(_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.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_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.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+        ulong multiplier = delta.Pow2MontgomeryModWindowed(_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.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
             _previousExponent = exponent;
             _hasState = true;
             isUnity = _currentMontgomery == _montgomeryOne;

PerfectNumbers.Core/FactorEntry.cs

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

PerfectNumbers.Core/Gpu/GpuContextPool.cs

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

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

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

PerfectNumbers.Core/Gpu/GpuUInt128.cs

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

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

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

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

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

PerfectNumbers.Core/MersenneDivisorCycles.cs

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

PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs

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

PerfectNumbers.Core/MersennePrimeFactorTester.cs

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

PerfectNumbers.Core/PartialFactorResult.cs

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

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

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

PerfectNumbers.Core/PrimeOrderCalculator.cs

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

PerfectNumbers.Core/ULongExtensions.cs

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