Merge pull request #249 from mlemanczyk/codex/optimize-heuristiccombinedprimetester-calculations

Use divisor square tables in heuristic combined prime tester

EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs

@@ -29,8 +29,7 @@ public class PrimeTesterBenchmarks
             return;
         }
 
-        CancellationToken cancellationToken = CancellationToken.None;
-        _ = _heuristicTester.IsPrimeCpu(_candidates[0], cancellationToken);
+        _ = _heuristicTester.IsPrimeCpu(_candidates[0]);
         _ = _heuristicTester.IsPrimeGpu(_candidates[0]);
         _ = Prime.Numbers.IsPrime(_candidates[0]);
     }
@@ -40,12 +39,11 @@ public class PrimeTesterBenchmarks
     {
         HeuristicPrimeTester tester = _heuristicTester;
         ulong[] values = _candidates;
-        CancellationToken cancellationToken = CancellationToken.None;
         int primeCount = 0;
 
         for (int i = 0; i < values.Length; i++)
         {
-            if (tester.IsPrimeCpu(values[i], cancellationToken))
+            if (tester.IsPrimeCpu(values[i]))
             {
                 primeCount++;
             }
@@ -75,7 +73,6 @@ public class PrimeTesterBenchmarks
     [Benchmark]
     public int HeuristicCombinedGpu()
     {
-        HeuristicCombinedPrimeTester tester = _heuristicCombinedTester;
         ulong[] values = _candidates;
         int primeCount = 0;
 
@@ -93,7 +90,6 @@ public class PrimeTesterBenchmarks
     [Benchmark]
     public int HeuristicCombinedCpu()
     {
-        HeuristicCombinedPrimeTester tester = _heuristicCombinedTester;
         ulong[] values = _candidates;
         int primeCount = 0;
 
@@ -115,12 +111,11 @@ public class PrimeTesterBenchmarks
 	{
 		PrimeTester tester = _nonHeuristicTester;
         ulong[] values = _candidates;
-        CancellationToken cancellationToken = CancellationToken.None;
         int primeCount = 0;
 
         for (int i = 0; i < values.Length; i++)
         {
-            if (tester.IsPrimeGpu(values[i], cancellationToken))
+            if (tester.IsPrimeGpu(values[i]))
             {
                 primeCount++;
             }
@@ -133,12 +128,11 @@ public class PrimeTesterBenchmarks
     public int NonHeuristicCpu()
     {
         ulong[] values = _candidates;
-        CancellationToken cancellationToken = CancellationToken.None;
         int primeCount = 0;
 
         for (int i = 0; i < values.Length; i++)
         {
-            if (PrimeTester.IsPrimeCpu(values[i], cancellationToken))
+            if (PrimeTester.IsPrime(values[i]))
             {
                 primeCount++;
             }

EvenPerfectBitScanner.Tests/ProgramTests.cs

@@ -55,7 +55,7 @@ public class ProgramTests
     [InlineData(81UL, false)]
     public void IsPrime_identifies_primes_correctly(ulong n, bool expected)
     {
-        PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected);
+        PrimeTester.IsPrime(n).Should().Be(expected);
     }
 
     [Fact]

EvenPerfectBitScanner/Program.cs

@@ -622,8 +622,8 @@ internal static class Program
 		if (!_cliArguments.UseByDivisor)
 		{
 			if (!(GpuContextPool.ForceCpu
-											? PrimeTester.IsPrime(p, CancellationToken.None)
-											: PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None)))
+					? PrimeTester.IsPrime(p)
+					: PrimeTesters.Value!.IsPrimeGpu(p)))
 			{
 				_lastCompositeP = true;
 				return false;

PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs

@@ -17,7 +17,7 @@ public class PrimeTesterGpuTests
             ulong[] primes = [31UL, 61UL, 89UL, 107UL, 127UL, 521UL];
             foreach (ulong prime in primes)
             {
-                tester.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue();
+                tester.IsPrimeGpu(prime).Should().BeTrue();
             }
         }
         finally
@@ -37,7 +37,7 @@ public class PrimeTesterGpuTests
             ulong[] composites = [33UL, 39UL, 51UL, 77UL, 91UL];
             foreach (ulong composite in composites)
             {
-                tester.IsPrimeGpu(composite, CancellationToken.None).Should().BeFalse();
+                tester.IsPrimeGpu(composite).Should().BeFalse();
             }
         }
         finally

PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs

@@ -19,7 +19,7 @@ public class HeuristicCombinedPrimeTesterTests
         Span<HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState> groupBState =
             stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4];
 
-        var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500, lastDigit, groupBState);
+        var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500UL * 500UL, lastDigit, groupBState);
         var results = new List<HeuristicCombinedPrimeTester.HeuristicDivisorCandidate>(32);
 
         while (results.Count < 20 && enumerator.TryGetNext(out var candidate))
@@ -230,14 +230,10 @@ public class HeuristicCombinedPrimeTesterTests
     [InlineData(18_446_744_073_709_551_615UL)]
     [InlineData(18_446_744_073_709_551_614UL)]
     [InlineData(18_446_744_073_709_551_613UL)]
-    public void ComputeHeuristicSqrt_returns_floor_square_root(ulong value)
+    public void ComputeHeuristicDivisorSquareLimit_returns_input_value(ulong value)
     {
-        ulong sqrt = HeuristicCombinedPrimeTester.ComputeHeuristicSqrt(value);
+        ulong limit = HeuristicCombinedPrimeTester.ComputeHeuristicDivisorSquareLimit(value);
 
-        bool squareBelowOrEqual = (UInt128)sqrt * sqrt <= (UInt128)value;
-        bool nextSquareAbove = (UInt128)(sqrt + 1UL) * (UInt128)(sqrt + 1UL) > (UInt128)value;
-
-        squareBelowOrEqual.Should().BeTrue();
-        nextSquareAbove.Should().BeTrue();
+        limit.Should().Be(value);
     }
 }

PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs

@@ -26,7 +26,7 @@ public class PrimeTesterIsPrimeTests
     {
         var tester = new PrimeTester();
 
-		PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected);
+		PrimeTester.IsPrime(n).Should().Be(expected);
     }
 
     [Theory]
@@ -40,7 +40,7 @@ public class PrimeTesterIsPrimeTests
         GpuContextPool.ForceCpu = false;
         var tester = new PrimeTester();
 
-        tester.IsPrimeGpu(n, CancellationToken.None).Should().Be(expected);
+        tester.IsPrimeGpu(n).Should().Be(expected);
     }
 
     [Fact]
@@ -50,9 +50,9 @@ public class PrimeTesterIsPrimeTests
         var tester = new PrimeTester();
         GpuContextPool.ForceCpu = true;
 
-        tester.IsPrimeGpu(137UL, CancellationToken.None).Should().BeTrue();
+        tester.IsPrimeGpu(137UL).Should().BeTrue();
 
-        tester.IsPrimeGpu(341UL, CancellationToken.None).Should().BeFalse();
+        tester.IsPrimeGpu(341UL).Should().BeFalse();
 
         GpuContextPool.ForceCpu = false;
     }
@@ -71,7 +71,7 @@ public class PrimeTesterIsPrimeTests
         var testers = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
         Parallel.For(0, primes.Length, i =>
         {
-            results[i] = testers.Value!.IsPrimeGpu(primes[i], CancellationToken.None);
+            results[i] = testers.Value!.IsPrimeGpu(primes[i]);
         });
 
         results.Should().AllBeEquivalentTo(true);
@@ -80,7 +80,7 @@ public class PrimeTesterIsPrimeTests
         {
             foreach (ulong prime in primes)
             {
-                testers.Value!.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue();
+                testers.Value!.IsPrimeGpu(prime).Should().BeTrue();
             }
         });
     }

PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs

@@ -48,7 +48,7 @@ internal static class PrimeTesterKernels
                 break;
         }
 
-        int length = XMath.Min((int)primes.Length, (int)primeSquares.Length);
+        int length = (int)primes.Length;
         for (int i = 0; i < length; i++)
         {
             ulong prime = primes[i];
@@ -68,16 +68,10 @@ internal static class PrimeTesterKernels
         results[index] = result;
     }
 
-    public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView<ulong> divisors, ulong n, ArrayView<byte> results)
+    public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView<ulong> divisors, ulong n, ArrayView<int> resultFlag)
     {
         ulong divisor = divisors[index];
-        if (divisor <= 1UL)
-        {
-            results[index] = 0;
-            return;
-        }
-
-        results[index] = n % divisor == 0UL ? (byte)1 : (byte)0;
+        Atomic.Or(ref resultFlag[0], divisor <= 1UL ? 0 : (n % divisor == 0UL ? 1 : 0));
     }
 
     public static void SharesFactorKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<byte> results)

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -292,6 +292,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                         coveredRange = true;
                         return true;
                     }
+
+                    divisorSpan[admissibleCount] = divisorValue;
+                    divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue);
+                    offsetSpan[admissibleCount] = admissibleCount;
+                    countSpan[admissibleCount] = 1;
+                    cycleSpan[admissibleCount] = divisorCycle;
+                    exponentSpan[admissibleCount] = prime;
+                    admissibleCount++;
                 }
 
                 nextDivisor128 += twoP128;
@@ -444,13 +452,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         // {
         //     return 0UL;
         // }
+        ulong cappedLimit = divisorLimit;
+        const ulong TestPrimeDivisorCap = 170_000_000UL;
+        if (prime < 10_000_000UL && cappedLimit > TestPrimeDivisorCap)
+        {
+            cappedLimit = TestPrimeDivisorCap;
+        }
+
         if (prime - 1UL >= 64UL)
         {
-            return divisorLimit;
+            return cappedLimit;
         }
 
         ulong computedLimit = (1UL << (int)(prime - 1UL)) - 1UL;
-        return computedLimit < divisorLimit ? computedLimit : divisorLimit;
+        if (prime < 10_000_000UL && computedLimit > cappedLimit)
+        {
+            computedLimit = cappedLimit;
+        }
+
+        return computedLimit < cappedLimit ? computedLimit : cappedLimit;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession

PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs

@@ -17,29 +17,33 @@ public sealed class HeuristicCombinedPrimeTester
 	{
 		HeuristicPrimeSieves.EnsureInitialized();
 
-		CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB);
-		CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB);
-		CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB);
-		CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB);
+		CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding1TwoAOneBSquares);
+		CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding3TwoAOneBSquares);
+		CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding7TwoAOneBSquares);
+		CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding9TwoAOneBSquares);
 
-		CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB);
-		CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB);
-		CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB);
-		CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB);
+		CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding1OneAOneBSquares);
+		CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding3OneAOneBSquares);
+		CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding7OneAOneBSquares);
+		CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding9OneAOneBSquares);
 
-		CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB);
-		CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB);
-		CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB);
-		CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB);
+		CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding1ThreeAOneBSquares);
+		CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding3ThreeAOneBSquares);
+		CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding7ThreeAOneBSquares);
+		CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding9ThreeAOneBSquares);
 
-		CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB);
-		CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB);
-		CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB);
-		CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB);
+		CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding1ThreeATwoBSquares);
+		CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding3ThreeATwoBSquares);
+		CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding7ThreeATwoBSquares);
+		CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding9ThreeATwoBSquares);
 		CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB;
 		CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB;
 		CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB;
 		CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB;
+		CombinedDivisorsEnding1Squares = CombinedDivisorsEnding1OneAOneBSquares;
+		CombinedDivisorsEnding3Squares = CombinedDivisorsEnding3OneAOneBSquares;
+		CombinedDivisorsEnding7Squares = CombinedDivisorsEnding7OneAOneBSquares;
+		CombinedDivisorsEnding9Squares = CombinedDivisorsEnding9OneAOneBSquares;
 	}
 
 	private const ulong Wheel210 = 210UL;
@@ -71,6 +75,31 @@ public sealed class HeuristicCombinedPrimeTester
 	private static readonly uint[] CombinedDivisorsEnding7;
 	private static readonly uint[] CombinedDivisorsEnding9;
 
+	private static readonly ulong[] CombinedDivisorsEnding1TwoAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding3TwoAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding7TwoAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding9TwoAOneBSquares;
+
+	private static readonly ulong[] CombinedDivisorsEnding1OneAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding3OneAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding7OneAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding9OneAOneBSquares;
+
+	private static readonly ulong[] CombinedDivisorsEnding1ThreeAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding3ThreeAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding7ThreeAOneBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding9ThreeAOneBSquares;
+
+	private static readonly ulong[] CombinedDivisorsEnding1ThreeATwoBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding3ThreeATwoBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding7ThreeATwoBSquares;
+	private static readonly ulong[] CombinedDivisorsEnding9ThreeATwoBSquares;
+
+	private static readonly ulong[] CombinedDivisorsEnding1Squares;
+	private static readonly ulong[] CombinedDivisorsEnding3Squares;
+	private static readonly ulong[] CombinedDivisorsEnding7Squares;
+	private static readonly ulong[] CombinedDivisorsEnding9Squares;
+
 	private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
 
 	internal enum CombinedDivisorPattern : byte
@@ -123,53 +152,51 @@ public sealed class HeuristicCombinedPrimeTester
 	public static bool IsPrimeCpu(ulong n)
 	{
 		byte nMod10 = (byte)n.Mod10();
-		ulong sqrtLimit = ComputeHeuristicSqrt(n);
+		ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
 
-		if (sqrtLimit < 3UL)
+		if (maxDivisorSquare < 9UL)
 		{
 			return EvaluateWithOpenNumericFallback(n);
 		}
 
-		return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10);
+		return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10);
 	}
 
 	public static bool IsPrimeGpu(ulong n)
 	{
 		byte nMod10 = (byte)n.Mod10();
-		ulong sqrtLimit = ComputeHeuristicSqrt(n);
-		return IsPrimeGpu(n, sqrtLimit, nMod10);
+		ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
+		return IsPrimeGpu(n, maxDivisorSquare, nMod10);
 	}
 
-	public static bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10)
+	public static bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10)
 	{
 		// TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
-		if (sqrtLimit < 3UL)
+		if (maxDivisorSquare < 9UL)
 		{
 			return EvaluateWithOpenNumericFallback(n);
 		}
 
-		bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10);
+		bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10);
 		return !compositeDetected;
 	}
 
-	private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10)
+	private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10)
 	{
 		ReadOnlySpan<uint> combinedDivisors = GetCombinedDivisors(nMod10);
-		if (combinedDivisors.IsEmpty)
-		{
-			return EvaluateWithOpenNumericFallback(n);
-		}
+		ReadOnlySpan<ulong> combinedDivisorSquares = GetCombinedDivisorSquares(nMod10);
 
 		int length = combinedDivisors.Length;
+
 		for (int i = 0; i < length; i++)
 		{
-			uint entry = combinedDivisors[i];
-			if (entry > sqrtLimit)
+			ulong divisorSquare = combinedDivisorSquares[i];
+			if (divisorSquare > maxDivisorSquare)
 			{
 				break;
 			}
 
-			if (n % entry == 0UL)
+			if (n % combinedDivisors[i] == 0UL)
 			{
 				return false;
 			}
@@ -178,127 +205,87 @@ public sealed class HeuristicCombinedPrimeTester
 		return EvaluateWithOpenNumericFallback(n);
 	}
 
-	private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10)
+	private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
 	{
-		int batchCapacity = HeuristicGpuDivisorBatchSize;
-		var divisorPool = ThreadStaticPools.UlongPool;
-		var hitPool = ThreadStaticPools.BytePool;
+	        int batchCapacity = HeuristicGpuDivisorBatchSize;
+	        var divisorPool = ThreadStaticPools.UlongPool;
 
-		ulong[] divisorArray = divisorPool.Rent(batchCapacity);
-		byte[] hitFlags = hitPool.Rent(batchCapacity);
+	        ulong[] divisorArray = divisorPool.Rent(batchCapacity);
 
-		var limiter = GpuPrimeWorkLimiter.Acquire();
-		var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent();
-		var accelerator = gpu.Accelerator;
-		var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator);
-		var scratch = state.RentScratch(batchCapacity, accelerator);
+	        var limiter = GpuPrimeWorkLimiter.Acquire();
+                var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity);
+                var accelerator = gpu.Accelerator;
+                var kernel = gpu.HeuristicTrialDivisionKernel;
 
-		bool compositeDetected = false;
-		int count = 0;
+                bool compositeDetected = false;
+                int count = 0;
 
-		var input = scratch.Input;
-		var output = scratch.Output;
+                var inputView = gpu.Input.View;
+                var flagView = gpu.HeuristicFlag.View;
 
-		bool ProcessBatch(int length)
-		{
-			input.View.CopyFromCPU(ref divisorArray[0], length);
-			state.HeuristicTrialDivisionKernel(length, input.View, n, output.View);
-			accelerator.Synchronize();
-			output.View.CopyToCPU(ref hitFlags[0], length);
+                bool ProcessBatch(int length)
+                {
+                        inputView.CopyFromCPU(ref divisorArray[0], length);
+                        int compositeFlag = 0;
+                        flagView.CopyFromCPU(ref compositeFlag, 1);
+                        kernel(length, inputView, n, flagView);
+                        accelerator.Synchronize();
+                        flagView.CopyToCPU(ref compositeFlag, 1);
 
-			for (int i = 0; i < length; i++)
-			{
-				if (hitFlags[i] == 0)
-				{
-					continue;
-				}
+                        return compositeFlag != 0;
+                }
 
-				return true;
-				// ulong divisor = divisorArray[i];
-				// if (divisor > 1UL && n % divisor == 0UL)
-				// {
-				// 	return true;
-				// }
-			}
+                ReadOnlySpan<uint> combinedDivisors = GetCombinedDivisors(nMod10);
+                ReadOnlySpan<ulong> combinedDivisorSquares = GetCombinedDivisorSquares(nMod10);
 
-			return false;
-		}
+                int length = combinedDivisors.Length;
 
-		ReadOnlySpan<uint> combinedDivisors = GetCombinedDivisors(nMod10);
-		if (combinedDivisors.IsEmpty)
-		{
-			goto Cleanup;
-		}
+                for (int i = 0; i < length; i++)
+                {
+                        ulong divisorSquare = combinedDivisorSquares[i];
+                        if (divisorSquare > maxDivisorSquare)
+                        {
+                                break;
+                        }
 
-		int length = combinedDivisors.Length;
-		for (int i = 0; i < length; i++)
-		{
-			uint entry = combinedDivisors[i];
-			if (entry > sqrtLimit)
-			{
-				break;
-			}
+                        divisorArray[count] = combinedDivisors[i];
+                        count++;
 
-			divisorArray[count] = entry;
-			count++;
+                        if (count == batchCapacity)
+                        {
+                                if (ProcessBatch(count))
+                                {
+                                        compositeDetected = true;
+                                        goto Cleanup;
+                                }
 
-			if (count == batchCapacity)
-			{
-				if (ProcessBatch(count))
-				{
-					compositeDetected = true;
-					goto Cleanup;
-				}
+                                count = 0;
+                        }
+                }
 
-				count = 0;
-			}
-		}
+                if (count > 0 && ProcessBatch(count))
+                {
+                        compositeDetected = true;
+                        goto Cleanup;
+                }
 
-		if (!compositeDetected && count > 0)
-		{
-			if (ProcessBatch(count))
-			{
-				compositeDetected = true;
-				goto Cleanup;
-			}
-		}
 
-	Cleanup:
-		state.ReturnScratch(scratch);
-		divisorPool.Return(divisorArray);
-		hitPool.Return(hitFlags);
-		gpu.Dispose();
-		limiter.Dispose();
-		return compositeDetected;
+Cleanup:
+	        divisorPool.Return(divisorArray);
+	        gpu.Dispose();
+	        limiter.Dispose();
+	        return compositeDetected;
 	}
 
+
 	private static bool EvaluateWithOpenNumericFallback(ulong n)
 	{
 		return Prime.Numbers.IsPrime(n);
 	}
-
-	internal static ulong ComputeHeuristicSqrt(ulong n)
-	{
-		ulong sqrt = (ulong)Math.Sqrt(n);
-		UInt128 square = (UInt128)sqrt * sqrt;
-
-		while (square > n)
-		{
-			sqrt--;
-			square = (UInt128)sqrt * sqrt;
-		}
-
-		ulong next = sqrt + 1UL;
-		UInt128 nextSquare = (UInt128)next * next;
-		while (nextSquare <= n)
-		{
-			sqrt = next;
-			next++;
-			nextSquare = (UInt128)next * next;
-		}
-
-		return sqrt;
-	}
+        internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n)
+        {
+                return n;
+        }
 
 	private static ReadOnlySpan<uint> GetCombinedDivisors(byte nMod10)
 	{
@@ -308,7 +295,19 @@ public sealed class HeuristicCombinedPrimeTester
 			3 => CombinedDivisorsEnding3,
 			7 => CombinedDivisorsEnding7,
 			9 => CombinedDivisorsEnding9,
-			_ => ReadOnlySpan<uint>.Empty,
+			_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
+		};
+	}
+
+	private static ReadOnlySpan<ulong> GetCombinedDivisorSquares(byte nMod10)
+	{
+		return nMod10 switch
+		{
+			1 => CombinedDivisorsEnding1Squares,
+			3 => CombinedDivisorsEnding3Squares,
+			7 => CombinedDivisorsEnding7Squares,
+			9 => CombinedDivisorsEnding9Squares,
+			_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
 		};
 	}
 
@@ -322,7 +321,7 @@ public sealed class HeuristicCombinedPrimeTester
 				3 => CombinedDivisorsEnding3TwoAOneB,
 				7 => CombinedDivisorsEnding7TwoAOneB,
 				9 => CombinedDivisorsEnding9TwoAOneB,
-				_ => ReadOnlySpan<uint>.Empty,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
 			},
 			CombinedDivisorPattern.OneAOneB => nMod10 switch
 			{
@@ -330,7 +329,7 @@ public sealed class HeuristicCombinedPrimeTester
 				3 => CombinedDivisorsEnding3OneAOneB,
 				7 => CombinedDivisorsEnding7OneAOneB,
 				9 => CombinedDivisorsEnding9OneAOneB,
-				_ => ReadOnlySpan<uint>.Empty,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
 			},
 			CombinedDivisorPattern.ThreeAOneB => nMod10 switch
 			{
@@ -338,7 +337,7 @@ public sealed class HeuristicCombinedPrimeTester
 				3 => CombinedDivisorsEnding3ThreeAOneB,
 				7 => CombinedDivisorsEnding7ThreeAOneB,
 				9 => CombinedDivisorsEnding9ThreeAOneB,
-				_ => ReadOnlySpan<uint>.Empty,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
 			},
 			CombinedDivisorPattern.ThreeATwoB => nMod10 switch
 			{
@@ -346,13 +345,53 @@ public sealed class HeuristicCombinedPrimeTester
 				3 => CombinedDivisorsEnding3ThreeATwoB,
 				7 => CombinedDivisorsEnding7ThreeATwoB,
 				9 => CombinedDivisorsEnding9ThreeATwoB,
-				_ => ReadOnlySpan<uint>.Empty,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."),
 			},
-			_ => ReadOnlySpan<uint>.Empty,
+			_ => throw new InvalidOperationException($"Unsupported combined divisor pattern: {pattern}."),
 		};
 	}
 
-	private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern)
+	internal static ReadOnlySpan<ulong> GetCombinedDivisorSquares(byte nMod10, CombinedDivisorPattern pattern)
+	{
+		return pattern switch
+		{
+			CombinedDivisorPattern.TwoAOneB => nMod10 switch
+			{
+				1 => CombinedDivisorsEnding1TwoAOneBSquares,
+				3 => CombinedDivisorsEnding3TwoAOneBSquares,
+				7 => CombinedDivisorsEnding7TwoAOneBSquares,
+				9 => CombinedDivisorsEnding9TwoAOneBSquares,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
+			},
+			CombinedDivisorPattern.OneAOneB => nMod10 switch
+			{
+				1 => CombinedDivisorsEnding1OneAOneBSquares,
+				3 => CombinedDivisorsEnding3OneAOneBSquares,
+				7 => CombinedDivisorsEnding7OneAOneBSquares,
+				9 => CombinedDivisorsEnding9OneAOneBSquares,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
+			},
+			CombinedDivisorPattern.ThreeAOneB => nMod10 switch
+			{
+				1 => CombinedDivisorsEnding1ThreeAOneBSquares,
+				3 => CombinedDivisorsEnding3ThreeAOneBSquares,
+				7 => CombinedDivisorsEnding7ThreeAOneBSquares,
+				9 => CombinedDivisorsEnding9ThreeAOneBSquares,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
+			},
+			CombinedDivisorPattern.ThreeATwoB => nMod10 switch
+			{
+				1 => CombinedDivisorsEnding1ThreeATwoBSquares,
+				3 => CombinedDivisorsEnding3ThreeATwoBSquares,
+				7 => CombinedDivisorsEnding7ThreeATwoBSquares,
+				9 => CombinedDivisorsEnding9ThreeATwoBSquares,
+				_ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."),
+			},
+			_ => throw new InvalidOperationException($"Unsupported combined divisor pattern for squares: {pattern}."),
+		};
+	}
+
+	private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern, out ulong[] squares)
 	{
 		ReadOnlySpan<int> groupA = HeuristicPrimeSieves.GroupADivisors;
 		ReadOnlySpan<uint> groupB = GetGroupBDivisors(nMod10);
@@ -412,7 +451,16 @@ public sealed class HeuristicCombinedPrimeTester
 			}
 		}
 
-		return combined.ToArray();
+		uint[] result = combined.ToArray();
+		ulong[] resultSquares = new ulong[result.Length];
+		for (int i = 0; i < result.Length; i++)
+		{
+			ulong divisor = result[i];
+			resultSquares[i] = divisor * divisor;
+		}
+
+		squares = resultSquares;
+		return result;
 	}
 
 
@@ -447,9 +495,9 @@ public sealed class HeuristicCombinedPrimeTester
 	};
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
+	internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
 	{
-		return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
+		return new HeuristicDivisorEnumerator(maxDivisorSquare, nMod10, groupBBuffer);
 	}
 
 	internal struct HeuristicGroupBSequenceState
@@ -530,44 +578,50 @@ public sealed class HeuristicCombinedPrimeTester
 	}
 
 	internal ref struct HeuristicDivisorEnumerator
-	{
-		private readonly ulong sqrtLimit;
-		private readonly ReadOnlySpan<uint> combinedDivisors;
-		private int index;
+        {
+                private readonly ulong maxDivisorSquare;
+                private readonly ReadOnlySpan<uint> combinedDivisors;
+                private readonly ReadOnlySpan<ulong> combinedDivisorSquares;
+                private int index;
 
-		public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
-		{
-			this.sqrtLimit = sqrtLimit;
-			_ = groupBBuffer;
-			combinedDivisors = GetCombinedDivisors(nMod10);
-			index = 0;
-		}
+                public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
+                {
+                        this.maxDivisorSquare = maxDivisorSquare;
+                        _ = groupBBuffer;
+                        combinedDivisors = GetCombinedDivisors(nMod10);
+                        combinedDivisorSquares = GetCombinedDivisorSquares(nMod10);
+                        index = 0;
+                }
 
-		public bool TryGetNext(out HeuristicDivisorCandidate candidate)
-		{
-			while (index < combinedDivisors.Length)
-			{
-				uint entry = combinedDivisors[index++];
-				if (entry > sqrtLimit)
-				{
-					index = combinedDivisors.Length;
-					break;
-				}
+                public bool TryGetNext(out HeuristicDivisorCandidate candidate)
+                {
+                        while (index < combinedDivisors.Length)
+                        {
+                                int currentIndex = index;
+                                ulong divisorSquare = combinedDivisorSquares[currentIndex];
+                                if (divisorSquare > maxDivisorSquare)
+                                {
+                                        index = combinedDivisors.Length;
+                                        break;
+                                }
 
-				HeuristicDivisorGroup group = ResolveGroup(entry);
-				ulong value = entry;
-				byte ending = (byte)(value % 10UL);
-				ushort residue = (ushort)(value % Wheel210);
-				candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue);
-				return true;
-			}
+                                uint entry = combinedDivisors[currentIndex];
+                                index = currentIndex + 1;
 
-			candidate = default;
-			return false;
-		}
-	}
+                                HeuristicDivisorGroup group = ResolveGroup(entry);
+                                ulong value = entry;
+                                byte ending = (byte)(value % 10UL);
+                                ushort residue = (ushort)(value % Wheel210);
+                                candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue);
+                                return true;
+                        }
 
-	internal ref struct MersenneHeuristicDivisorEnumerator
+                        candidate = default;
+                        return false;
+                }
+        }
+
+        internal ref struct MersenneHeuristicDivisorEnumerator
 	{
 		private readonly GpuUInt128 step;
 		private readonly GpuUInt128 limit;

PerfectNumbers.Core/HeuristicPrimeTester.cs

@@ -21,8 +21,6 @@ public sealed class HeuristicPrimeTester
     }
 
     // TODO: Don't use these tiny arrays. Hard-code the values / checks instead
-    private static readonly ulong[] GroupAConstantDivisors = [3UL, 7UL, 11UL, 13UL];
-    private static readonly byte[] GroupAIncrementPattern = [20, 10];
     private static readonly byte[] GroupBEndingOrderMod1 = [9, 1];
     private static readonly byte[] GroupBEndingOrderMod3 = [9, 7];
     private static readonly byte[] GroupBEndingOrderMod7 = [7, 1];
@@ -30,10 +28,8 @@ public sealed class HeuristicPrimeTester
 
     // TODO: Would it speed up things if we implement something like ExponentRemainderStepper,
     // MersenneDivisorResidueStepper, or CycleRemainderStepper work with these?
-    private static readonly ushort[] Wheel210ResiduesEnding1 = [1, 11, 31, 41, 61, 71, 101, 121, 131, 151, 181, 191];
-    private static readonly ushort[] Wheel210ResiduesEnding7 = [17, 37, 47, 67, 97, 107, 127, 137, 157, 167, 187, 197];
-    private static readonly ushort[] Wheel210ResiduesEnding9 = [19, 29, 59, 79, 89, 109, 139, 149, 169, 179, 199, 209];
     private const ulong Wheel210 = 210UL;
+    private const int GroupAConstantCount = 4;
     internal const int MaxGroupBSequences = 36;
 
     private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B.
@@ -81,76 +77,73 @@ public sealed class HeuristicPrimeTester
 
     public static HeuristicPrimeTester Exclusive => _tester ??= new();
 
-    public bool IsPrimeCpu(ulong n, CancellationToken ct)
+    public bool IsPrimeCpu(ulong n)
     {
-        ct.ThrowIfCancellationRequested();
-
         byte nMod10 = (byte)n.Mod10();
-        ulong sqrtLimit = ComputeHeuristicSqrt(n);
+        ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
 
-        if (sqrtLimit < 3UL)
+        if (maxDivisorSquare < 9UL)
         {
             return EvaluateWithOpenNumericFallback(n);
         }
 
         bool includeGroupB = UseHeuristicGroupBTrialDivision;
-        return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB);
-    }
-
-    public bool IsPrimeCpu(ulong n)
-    {
-        return IsPrimeCpu(n, CancellationToken.None);
+        return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB);
     }
 
     public bool IsPrimeGpu(ulong n)
     {
         byte nMod10 = (byte)n.Mod10();
-        ulong sqrtLimit = ComputeHeuristicSqrt(n);
-        return IsPrimeGpu(n, sqrtLimit, nMod10);
+        ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n);
+        return IsPrimeGpu(n, maxDivisorSquare, nMod10);
     }
 
-    public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10)
+    public bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10)
     {
-        return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10);
+        return HeuristicIsPrimeGpuCore(n, maxDivisorSquare, nMod10);
     }
 
-    private bool HeuristicIsPrimeGpuCore(ulong n, ulong sqrtLimit, byte nMod10)
+    private bool HeuristicIsPrimeGpuCore(ulong n, ulong maxDivisorSquare, byte nMod10)
     {
         // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
-        if (sqrtLimit < 3UL)
+        if (maxDivisorSquare < 9UL)
         {
             return EvaluateWithOpenNumericFallback(n);
         }
 
         // if (!UseHeuristicGroupBTrialDivision)
         // {
-        //     return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false);
+        //     return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB: false);
         // }
 
-        bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10);
+        bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10);
         return !compositeDetected;
     }
 
-    private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB)
+
+    private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10, bool includeGroupB)
     {
         ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
+        ReadOnlySpan<ulong> groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares;
         int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize);
         int groupAIndex = 0;
 
         if (!includeGroupB)
         {
-            while (groupAIndex < groupADivisors.Length)
+            int groupALength = groupADivisors.Length;
+            while (groupAIndex < groupALength)
             {
                 int processed = 0;
-                while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
+                while (processed < interleaveBatchSize && groupAIndex < groupALength)
                 {
-                    ulong divisor = (ulong)groupADivisors[groupAIndex];
-                    if (divisor > sqrtLimit)
+                    ulong divisorSquare = groupADivisorSquares[groupAIndex];
+                    if (divisorSquare > maxDivisorSquare)
                     {
-                        groupAIndex = groupADivisors.Length;
+                        groupAIndex = groupALength;
                         break;
                     }
 
+                    ulong divisor = (ulong)groupADivisors[groupAIndex];
                     if (n % divisor == 0UL)
                     {
                         return false;
@@ -160,14 +153,14 @@ public sealed class HeuristicPrimeTester
                     processed++;
                 }
 
-                if (groupAIndex >= groupADivisors.Length)
+                if (groupAIndex >= groupALength)
                 {
                     break;
                 }
 
-                if ((ulong)groupADivisors[groupAIndex] > sqrtLimit)
+                if (groupADivisorSquares[groupAIndex] > maxDivisorSquare)
                 {
-                    groupAIndex = groupADivisors.Length;
+                    groupAIndex = groupALength;
                     break;
                 }
             }
@@ -185,8 +178,8 @@ public sealed class HeuristicPrimeTester
             ? stackalloc int[endingOrder.Length]
             : new int[endingOrder.Length];
 
-        bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
-        bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit);
+        bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare;
+        bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
 
         while (groupAHasMore || groupBHasMore)
         {
@@ -195,14 +188,15 @@ public sealed class HeuristicPrimeTester
                 int processed = 0;
                 while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
                 {
-                    ulong divisor = (ulong)groupADivisors[groupAIndex];
-                    if (divisor > sqrtLimit)
+                    ulong divisorSquare = groupADivisorSquares[groupAIndex];
+                    if (divisorSquare > maxDivisorSquare)
                     {
                         groupAIndex = groupADivisors.Length;
                         groupAHasMore = false;
                         break;
                     }
 
+                    ulong divisor = (ulong)groupADivisors[groupAIndex];
                     if (n % divisor == 0UL)
                     {
                         return false;
@@ -218,7 +212,7 @@ public sealed class HeuristicPrimeTester
                 }
                 else
                 {
-                    groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
+                    groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare;
                     if (!groupAHasMore)
                     {
                         groupAIndex = groupADivisors.Length;
@@ -231,7 +225,7 @@ public sealed class HeuristicPrimeTester
                 int processed = 0;
                 while (processed < interleaveBatchSize)
                 {
-                    if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor))
+                    if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor))
                     {
                         groupBHasMore = false;
                         break;
@@ -247,7 +241,7 @@ public sealed class HeuristicPrimeTester
 
                 if (groupBHasMore)
                 {
-                    groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit);
+                    groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
                 }
             }
         }
@@ -258,198 +252,162 @@ public sealed class HeuristicPrimeTester
 
 
 
-    private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10)
+
+
+
+    private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
     {
         int batchCapacity = HeuristicGpuDivisorBatchSize;
         int interleaveBatchSize = HeuristicDivisorInterleaveBatchSize;
         var divisorPool = ThreadStaticPools.UlongPool;
-        var hitPool = ThreadStaticPools.BytePool;
 
         ulong[] divisorArray = divisorPool.Rent(batchCapacity);
-        byte[] hitFlags = hitPool.Rent(batchCapacity);
 
         var limiter = GpuPrimeWorkLimiter.Acquire();
-        var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent();
+        var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity);
         var accelerator = gpu.Accelerator;
-        var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator);
-        var scratch = state.RentScratch(batchCapacity, accelerator);
+        var kernel = gpu.HeuristicTrialDivisionKernel;
 
         bool compositeDetected = false;
         int count = 0;
 
-        // lock (gpu.ExecutionLock)
+        var inputView = gpu.Input.View;
+        var flagView = gpu.HeuristicFlag.View;
+
+        bool ProcessBatch(int length)
         {
-            var input = scratch.Input;
-            var output = scratch.Output;
-
-            bool ProcessBatch(int length)
+            if (length <= 0)
             {
-                if (length <= 0)
-                {
-                    return false;
-                }
-
-                input.View.CopyFromCPU(ref divisorArray[0], length);
-                state.HeuristicTrialDivisionKernel(length, input.View, n, output.View);
-                accelerator.Synchronize();
-                output.View.CopyToCPU(ref hitFlags[0], length);
-
-                for (int i = 0; i < length; i++)
-                {
-                    if (hitFlags[i] == 0)
-                    {
-                        continue;
-                    }
-
-                    ulong divisor = divisorArray[i];
-                    if (divisor > 1UL && n % divisor == 0UL)
-                    {
-                        return true;
-                    }
-                }
-
                 return false;
             }
 
-            ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
-            int groupAIndex = 0;
-            ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
-            Span<int> indices = endingOrder.Length <= 8
-                ? stackalloc int[endingOrder.Length]
-                : new int[endingOrder.Length];
+            inputView.CopyFromCPU(ref divisorArray[0], length);
+            int compositeFlag = 0;
+            flagView.CopyFromCPU(ref compositeFlag, 1);
+            kernel(length, inputView, n, flagView);
+            accelerator.Synchronize();
+            flagView.CopyToCPU(ref compositeFlag, 1);
+            return compositeFlag != 0;
+        }
 
-            bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
-            bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, sqrtLimit);
+        ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
+        ReadOnlySpan<ulong> groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares;
+        int groupAIndex = 0;
+        ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
+        Span<int> indices = endingOrder.Length <= 8
+            ? stackalloc int[endingOrder.Length]
+            : new int[endingOrder.Length];
 
-            while (!compositeDetected && (groupAHasMore || groupBHasMore))
+        bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare;
+        bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
+
+        while (!compositeDetected && (groupAHasMore || groupBHasMore))
+        {
+            if (groupAHasMore)
             {
-                if (groupAHasMore)
+                int processed = 0;
+                while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
                 {
-                    int processed = 0;
-                    while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
+                    ulong divisorSquare = groupADivisorSquares[groupAIndex];
+                    if (divisorSquare > maxDivisorSquare)
                     {
-                        ulong divisor = (ulong)groupADivisors[groupAIndex];
-                        if (divisor > sqrtLimit)
+                        groupAIndex = groupADivisors.Length;
+                        groupAHasMore = false;
+                        break;
+                    }
+
+                    divisorArray[count] = (ulong)groupADivisors[groupAIndex];
+                    count++;
+
+                    if (count == batchCapacity)
+                    {
+                        if (ProcessBatch(count))
                         {
-                            groupAIndex = groupADivisors.Length;
-                            groupAHasMore = false;
+                            compositeDetected = true;
                             break;
                         }
 
-                        divisorArray[count] = divisor;
-                        count++;
+                        count = 0;
+                    }
 
-                        if (count == batchCapacity)
+                    groupAIndex++;
+                    processed++;
+                }
+
+                if (compositeDetected)
+                {
+                    break;
+                }
+
+                if (groupAIndex >= groupADivisors.Length)
+                {
+                    groupAHasMore = false;
+                }
+                else
+                {
+                    groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare;
+                    if (!groupAHasMore)
+                    {
+                        groupAIndex = groupADivisors.Length;
+                    }
+                }
+            }
+
+            if (!compositeDetected && groupBHasMore)
+            {
+                int processed = 0;
+                while (processed < interleaveBatchSize)
+                {
+                    if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor))
+                    {
+                        groupBHasMore = false;
+                        break;
+                    }
+
+                    divisorArray[count] = divisor;
+                    count++;
+
+                    if (count == batchCapacity)
+                    {
+                        if (ProcessBatch(count))
                         {
-                            if (ProcessBatch(count))
-                            {
-                                compositeDetected = true;
-                                goto Cleanup;
-                            }
-
-                            count = 0;
+                            compositeDetected = true;
+                            break;
                         }
 
-                        groupAIndex++;
-                        processed++;
+                        count = 0;
                     }
 
-                    if (groupAIndex >= groupADivisors.Length)
-                    {
-                        groupAHasMore = false;
-                    }
-                    else
-                    {
-                        groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
-                        if (!groupAHasMore)
-                        {
-                            groupAIndex = groupADivisors.Length;
-                        }
-                    }
+                    processed++;
                 }
 
                 if (!compositeDetected && groupBHasMore)
                 {
-                    int processed = 0;
-                    while (processed < interleaveBatchSize)
-                    {
-                        if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor))
-                        {
-                            groupBHasMore = false;
-                            break;
-                        }
-
-                        divisorArray[count] = divisor;
-                        count++;
-
-                        if (count == batchCapacity)
-                        {
-                            if (ProcessBatch(count))
-                            {
-                                compositeDetected = true;
-                                goto Cleanup;
-                            }
-
-                            count = 0;
-                        }
-
-                        processed++;
-                    }
-
-                    if (groupBHasMore)
-                    {
-                        groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit);
-                    }
+                    groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare);
                 }
             }
-
-            if (!compositeDetected && count > 0)
-            {
-                if (ProcessBatch(count))
-                {
-                    compositeDetected = true;
-                    goto Cleanup;
-                }
-
-                count = 0;
-            }
         }
 
-Cleanup:
-        state.ReturnScratch(scratch);
+        if (!compositeDetected && count > 0 && ProcessBatch(count))
+        {
+            compositeDetected = true;
+        }
+
         divisorPool.Return(divisorArray);
-        hitPool.Return(hitFlags);
         gpu.Dispose();
         limiter.Dispose();
         return compositeDetected;
     }
 
+
+
     private static bool EvaluateWithOpenNumericFallback(ulong n)
     {
         return Prime.Numbers.IsPrime(n);
     }
-
-    internal static ulong ComputeHeuristicSqrt(ulong n)
+    internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n)
     {
-        ulong sqrt = (ulong)Math.Sqrt(n);
-        UInt128 square = (UInt128)sqrt * sqrt;
-
-        while (square > n)
-        {
-            sqrt--;
-            square = (UInt128)sqrt * sqrt;
-        }
-
-        ulong next = sqrt + 1UL;
-        UInt128 nextSquare = (UInt128)next * next;
-        while (nextSquare <= n)
-        {
-            sqrt = next;
-            next++;
-            nextSquare = (UInt128)next * next;
-        }
-
-        return sqrt;
+        return n;
     }
 
     private static ReadOnlySpan<byte> GetGroupBEndingOrder(byte nMod10) => nMod10 switch
@@ -462,8 +420,9 @@ Cleanup:
     };
 
 
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static bool TrySelectNextGroupBDivisor(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong sqrtLimit, out ulong divisor)
+    private static bool TrySelectNextGroupBDivisor(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong maxDivisorSquare, out ulong divisor)
     {
         ulong bestCandidate = ulong.MaxValue;
         int bestEndingIndex = -1;
@@ -471,19 +430,21 @@ Cleanup:
         for (int i = 0; i < endingOrder.Length; i++)
         {
             ReadOnlySpan<uint> divisors = GetGroupBDivisors(endingOrder[i]);
+            ReadOnlySpan<ulong> squares = GetGroupBDivisorSquares(endingOrder[i]);
             int index = indices[i];
             if ((uint)index >= (uint)divisors.Length)
             {
                 continue;
             }
 
-            ulong candidate = (ulong)divisors[index];
-            if (candidate > sqrtLimit)
+            ulong candidateSquare = squares[index];
+            if (candidateSquare > maxDivisorSquare)
             {
                 indices[i] = divisors.Length;
                 continue;
             }
 
+            ulong candidate = (ulong)divisors[index];
             if (bestEndingIndex == -1 || candidate < bestCandidate)
             {
                 bestCandidate = candidate;
@@ -497,25 +458,29 @@ Cleanup:
             return false;
         }
 
-        divisor = bestCandidate;
+        ReadOnlySpan<uint> selectedDivisors = GetGroupBDivisors(endingOrder[bestEndingIndex]);
+        divisor = (ulong)selectedDivisors[indices[bestEndingIndex]];
         indices[bestEndingIndex]++;
         return true;
     }
 
+
+
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static bool HasGroupBCandidates(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong sqrtLimit)
+    private static bool HasGroupBCandidates(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong maxDivisorSquare)
     {
         for (int i = 0; i < endingOrder.Length; i++)
         {
             ReadOnlySpan<uint> divisors = GetGroupBDivisors(endingOrder[i]);
+            ReadOnlySpan<ulong> squares = GetGroupBDivisorSquares(endingOrder[i]);
             int index = indices[i];
             if ((uint)index >= (uint)divisors.Length)
             {
                 continue;
             }
 
-            ulong candidate = (ulong)divisors[index];
-            if (candidate > sqrtLimit)
+            if (squares[index] > maxDivisorSquare)
             {
                 indices[i] = divisors.Length;
                 continue;
@@ -527,6 +492,8 @@ Cleanup:
         return false;
     }
 
+
+
     private static ReadOnlySpan<uint> GetGroupBDivisors(byte ending) => ending switch
     {
         1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree,
@@ -536,35 +503,32 @@ Cleanup:
         _ => ReadOnlySpan<uint>.Empty,
     };
 
-    private static ReadOnlySpan<ushort> GetWheel210ResiduesForEnding(byte ending) => ending switch
+    private static ReadOnlySpan<ulong> GetGroupBDivisorSquares(byte ending) => ending switch
     {
-        1 => Wheel210ResiduesEnding1,
-        7 => Wheel210ResiduesEnding7,
-        9 => Wheel210ResiduesEnding9,
-        _ => ReadOnlySpan<ushort>.Empty,
+        1 => DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree,
+        3 => ReadOnlySpan<ulong>.Empty,
+        7 => DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree,
+        9 => DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree,
+        _ => ReadOnlySpan<ulong>.Empty,
     };
 
     private static int InitializeGroupBStates(byte nMod10, Span<HeuristicGroupBSequenceState> buffer)
     {
         ReadOnlySpan<byte> endings = GetGroupBEndingOrder(nMod10);
-        int count = 0;
+        int count = endings.Length;
 
-        for (int priorityIndex = 0; priorityIndex < endings.Length; priorityIndex++)
+        for (int i = 0; i < count; i++)
         {
-            ReadOnlySpan<ushort> residues = GetWheel210ResiduesForEnding(endings[priorityIndex]);
-            for (int residueIndex = 0; residueIndex < residues.Length; residueIndex++)
-            {
-                buffer[count++] = new HeuristicGroupBSequenceState(residues[residueIndex], endings[priorityIndex], (byte)priorityIndex);
-            }
+            buffer[i] = new HeuristicGroupBSequenceState(endings[i], (byte)i);
         }
 
         return count;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
+    internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
     {
-        return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
+        return new HeuristicDivisorEnumerator(maxDivisorSquare, nMod10, groupBBuffer);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -640,158 +604,143 @@ Cleanup:
         return false;
     }
 
-    internal ref struct HeuristicDivisorEnumerator
+
+internal ref struct HeuristicDivisorEnumerator
+{
+    private readonly ulong maxDivisorSquare;
+    private readonly ReadOnlySpan<int> groupADivisors;
+    private readonly ReadOnlySpan<ulong> groupADivisorSquares;
+    private int groupAIndex;
+    private Span<HeuristicGroupBSequenceState> groupBStates;
+
+    public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
     {
-        private readonly ulong sqrtLimit;
-        private int groupAConstantIndex;
-        private ulong groupAWheelCandidate;
-        private int groupAIncrementIndex;
-        private bool groupAWheelActive;
-        private Span<HeuristicGroupBSequenceState> groupBStates;
+        this.maxDivisorSquare = maxDivisorSquare;
+        groupADivisors = HeuristicPrimeSieves.GroupADivisors;
+        groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares;
+        groupAIndex = 0;
 
-        public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
+        int count = InitializeGroupBStates(nMod10, groupBBuffer);
+        groupBStates = count == 0 ? Span<HeuristicGroupBSequenceState>.Empty : groupBBuffer[..count];
+    }
+
+    public bool TryGetNext(out HeuristicDivisorCandidate candidate)
+    {
+        while (groupAIndex < groupADivisors.Length)
         {
-            this.sqrtLimit = sqrtLimit;
-            groupAConstantIndex = 0;
-            groupAWheelCandidate = 23UL;
-            groupAIncrementIndex = 0;
-            groupAWheelActive = sqrtLimit >= 23UL;
-
-            int count = InitializeGroupBStates(nMod10, groupBBuffer);
-            groupBStates = count == 0 ? Span<HeuristicGroupBSequenceState>.Empty : groupBBuffer[..count];
-
-            for (int i = 0; i < groupBStates.Length; i++)
+            ulong divisorSquare = groupADivisorSquares[groupAIndex];
+            if (divisorSquare > maxDivisorSquare)
             {
-                groupBStates[i].TryNormalize();
-            }
-        }
-
-        public bool TryGetNext(out HeuristicDivisorCandidate candidate)
-		{
-			// TODO: Hard-code constant A divisor checks here instead of using an array
-            while (groupAConstantIndex < GroupAConstantDivisors.Length)
-            {
-                ulong value = GroupAConstantDivisors[groupAConstantIndex++];
-				if (value > sqrtLimit)
-				{
-					groupAConstantIndex = GroupAConstantDivisors.Length;
-					break;
-				}
-
-				// TODO: We're constantly dividing constants mod 10 which is constant too.
-                candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAConstant, (byte)(value % 10UL), 0, (ushort)value);
-                return true;
+                groupAIndex = groupADivisors.Length;
+                break;
             }
 
-            if (groupAWheelActive)
-            {
-                ulong value = groupAWheelCandidate;
-                if (value <= sqrtLimit)
-                {
-                    candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAWheel, 3, 0, (ushort)(value % Wheel210));
-                    byte increment = GroupAIncrementPattern[groupAIncrementIndex];
-                    groupAIncrementIndex ^= 1;
+            int currentIndex = groupAIndex;
+            ulong divisor = (ulong)groupADivisors[currentIndex];
+            groupAIndex++;
 
-                    ulong next = value + increment;
-                    if (next <= value)
-                    {
-                        groupAWheelActive = false;
-                    }
-                    else
-                    {
-                        groupAWheelCandidate = next;
-                    }
+            HeuristicDivisorGroup group = currentIndex < GroupAConstantCount
+                ? HeuristicDivisorGroup.GroupAConstant
+                : HeuristicDivisorGroup.GroupAWheel;
 
-                    return true;
-                }
-
-                groupAWheelActive = false;
-            }
-
-            ulong bestCandidate = ulong.MaxValue;
-            int bestIndex = -1;
-
-            for (int i = 0; i < groupBStates.Length; i++)
-            {
-                ref HeuristicGroupBSequenceState state = ref groupBStates[i];
-                ulong value = state.Candidate;
-                if (value == ulong.MaxValue)
-                {
-                    continue;
-                }
-
-                if (value > sqrtLimit)
-                {
-                    state.Candidate = ulong.MaxValue;
-                    continue;
-                }
-
-                if (bestIndex == -1 || value < bestCandidate ||
-                    (value == bestCandidate && state.PriorityIndex < groupBStates[bestIndex].PriorityIndex))
-                {
-                    bestCandidate = value;
-                    bestIndex = i;
-                }
-            }
-
-            if (bestIndex == -1)
-            {
-                candidate = default;
-                return false;
-            }
-
-            ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex];
-            var result = new HeuristicDivisorCandidate(bestState.Candidate, HeuristicDivisorGroup.GroupB, bestState.Ending, bestState.PriorityIndex, bestState.Residue);
-            bestState.Advance();
-            candidate = result;
+            candidate = new HeuristicDivisorCandidate(
+                divisor,
+                group,
+                (byte)(divisor % 10UL),
+                0,
+                (ushort)(divisor % Wheel210));
             return true;
         }
-    }
 
-    internal struct HeuristicGroupBSequenceState(ushort residue, byte ending, byte priorityIndex)
-	{
-        public ulong Candidate = residue;
-        public byte PriorityIndex = priorityIndex;
-        public byte Ending = ending;
-        public ushort Residue = residue;
+        ulong bestCandidate = ulong.MaxValue;
+        int bestIndex = -1;
 
-		public void TryNormalize()
+        for (int i = 0; i < groupBStates.Length; i++)
         {
-            ulong candidate = Candidate;
-            if (candidate == ulong.MaxValue)
+            ref HeuristicGroupBSequenceState state = ref groupBStates[i];
+            if (!state.TryPeek(maxDivisorSquare, out ulong candidateValue))
             {
-                return;
+                continue;
             }
 
-            while (candidate <= 13UL)
+            if (bestIndex == -1 || candidateValue < bestCandidate ||
+                (candidateValue == bestCandidate && state.PriorityIndex < groupBStates[bestIndex].PriorityIndex))
             {
-                ulong next = candidate + Wheel210;
-                if (next <= candidate)
-                {
-                    Candidate = ulong.MaxValue;
-                    return;
-                }
-
-                candidate = next;
+                bestCandidate = candidateValue;
+                bestIndex = i;
             }
-
-            Candidate = candidate;
         }
 
-        public void Advance()
+        if (bestIndex == -1)
         {
-            ulong candidate = Candidate;
-            if (candidate == ulong.MaxValue)
-            {
-                return;
-            }
-
-            ulong next = candidate + Wheel210;
-            Candidate = next > candidate ? next : ulong.MaxValue;
+            candidate = default;
+            return false;
         }
 
-
+        ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex];
+        bestState.Advance();
+        candidate = new HeuristicDivisorCandidate(
+            bestCandidate,
+            HeuristicDivisorGroup.GroupB,
+            bestState.Ending,
+            bestState.PriorityIndex,
+            (ushort)(bestCandidate % Wheel210));
+        return true;
     }
+}
+
+
+
+internal struct HeuristicGroupBSequenceState
+{
+    public byte Ending;
+    public byte PriorityIndex;
+    public int Index;
+
+    public HeuristicGroupBSequenceState(byte ending, byte priorityIndex)
+    {
+        Ending = ending;
+        PriorityIndex = priorityIndex;
+        Index = 0;
+    }
+
+    public bool TryPeek(ulong maxDivisorSquare, out ulong divisor)
+    {
+        ReadOnlySpan<uint> divisors = GetGroupBDivisors(Ending);
+        ReadOnlySpan<ulong> squares = GetGroupBDivisorSquares(Ending);
+        int index = Index;
+        if ((uint)index >= (uint)divisors.Length)
+        {
+            divisor = 0UL;
+            return false;
+        }
+
+        if (squares[index] > maxDivisorSquare)
+        {
+            Index = divisors.Length;
+            divisor = 0UL;
+            return false;
+        }
+
+        divisor = (ulong)divisors[index];
+        return true;
+    }
+
+    public void Advance()
+    {
+        int next = Index + 1;
+        ReadOnlySpan<uint> divisors = GetGroupBDivisors(Ending);
+        if ((uint)next >= (uint)divisors.Length)
+        {
+            Index = divisors.Length;
+        }
+        else
+        {
+            Index = next;
+        }
+    }
+}
+
 
     internal ref struct MersenneHeuristicDivisorEnumerator
     {

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -411,8 +411,8 @@ public class MersenneDivisorCycles
 		// Atomic.Add(ref _tryFactorIntoCountsInternalHits, 1UL);
 		// Console.WriteLine($"MersenneDivisorCycles.TryFactorIntoCountsInternal hits {Volatile.Read(ref _tryFactorIntoCountsInternalHits)}");
 
-		HeuristicPrimeTester primeTester = _primeTester ??= new();
-		if (primeTester.IsPrimeGpu(remaining))
+		// HeuristicPrimeTester primeTester = _primeTester ??= new();
+		if (PrimeTester.IsPrime(remaining))
 		// if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None))
 		{
 			AddFactor(counts, remaining);
@@ -917,8 +917,8 @@ public class MersenneDivisorCycles
         return CalculateCycleLengthFallback(divisor);
     }
 
-	[ThreadStatic]
-	private static HeuristicPrimeTester? _primeTester;
+	// [ThreadStatic]
+	// private static HeuristicPrimeTester? _primeTester;
 
 	internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false)
 	{
@@ -940,10 +940,10 @@ public class MersenneDivisorCycles
 		// 	Console.WriteLine($"MersenneDivisorCycles.CalculateCycleLengthFallback hits {Volatile.Read(ref _tryCalculateCycleLengthHeuristicHits)}");			
 		// }
 
-		HeuristicPrimeTester primeTester = _primeTester ??= new();
+		// HeuristicPrimeTester primeTester = _primeTester ??= new();
 		// if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None))
-		if (!skipPrimeOrderHeuristic && primeTester.IsPrimeGpu(divisor))
-		// if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeCpu(divisor, CancellationToken.None))
+		// if (!skipPrimeOrderHeuristic && primeTester.IsPrimeGpu(divisor))
+		if (!skipPrimeOrderHeuristic && PrimeTester.IsPrime(divisor))
 		{
 			ulong computedOrder = PrimeOrderCalculator.Calculate(
 					divisor,

PerfectNumbers.Core/MersennePrimeFactorTester.cs

@@ -34,7 +34,7 @@ public static class MersennePrimeFactorTester
 		// Atomic.Add(ref _isPrimeFactorHits, 1UL);
 		// Console.WriteLine($"MersennePrimeFactorTester.IsPrimeFactor hits {Volatile.Read(ref _isPrimeFactorHits)}");
 
-        if (!PrimeTester.IsPrimeCpu(q, ct))
+        if (!PrimeTester.IsPrime(q))
         {
             return false;
         }
@@ -247,7 +247,7 @@ public static class MersennePrimeFactorTester
 			// Console.WriteLine($"MersennePrimeFactorTester.Factor64 hits {Volatile.Read(ref _factor64Hits)}");
 
             // if (HeuristicPrimeTester.Exclusive.IsPrime(m, ct))
-            if (PrimeTester.IsPrimeCpu(m, ct))
+            if (PrimeTester.IsPrime(m))
             {
                 dict[m] = 1;
                 continue;

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs

@@ -1079,7 +1079,7 @@ internal static partial class PrimeOrderCalculator
 		PartialFactorResult result;
 		PrimeTester primeTester = PrimeTester.Exclusive;
 		CancellationToken ct = CancellationToken.None;
-		HeuristicPrimeTester? tester = _tester ??= new();
+		// HeuristicPrimeTester? tester = _tester ??= new();
 
 		if (!gpuFactored)
 		{
@@ -1152,7 +1152,7 @@ internal static partial class PrimeOrderCalculator
 					// Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}");
 
 					// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None);
-					bool isPrime = tester.IsPrimeGpu(composite);
+					bool isPrime = PrimeTester.IsPrime(composite);
 					// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
 
 					if (isPrime)
@@ -1220,7 +1220,7 @@ internal static partial class PrimeOrderCalculator
 				// Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}");
 
 				// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None);
-				bool isPrime = tester.IsPrimeGpu(composite);
+				bool isPrime = PrimeTester.IsPrime(composite);
 				// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
 
 				entry = entry.WithPrimality(isPrime);
@@ -1257,7 +1257,7 @@ internal static partial class PrimeOrderCalculator
 			// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
 			// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
 
-			cofactorIsPrime = tester.IsPrimeGpu(cofactor);
+			cofactorIsPrime = PrimeTester.IsPrime(cofactor);
 			// cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor);
 			// cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor);
 		}
@@ -1782,8 +1782,8 @@ internal static partial class PrimeOrderCalculator
 
 	// private static ulong _factorCompletelyHits;
 	// private static ulong _factorCompletelyPollardRhoHits;
-	[ThreadStatic]
-	private static HeuristicPrimeTester? _tester;
+	// [ThreadStatic]
+	// private static HeuristicPrimeTester? _tester;
 
 	private static void FactorCompletely(ulong value, Dictionary<ulong, int> counts, bool knownComposite)
 	{
@@ -1794,10 +1794,10 @@ internal static partial class PrimeOrderCalculator
 
 		// Atomic.Add(ref _factorCompletelyHits, 1UL);
 		// Console.WriteLine($"Factor completely hits {Volatile.Read(ref _factorCompletelyHits)}");
-		HeuristicPrimeTester tester = _tester ??= new();
+		// HeuristicPrimeTester tester = _tester ??= new();
 
 		// bool isPrime = primeTester.HeuristicIsPrimeGpu(value);
-		bool isPrime = tester.IsPrimeGpu(value);
+		bool isPrime = PrimeTester.IsPrime(value);
 		// bool isPrime = PrimeTester.IsPrimeGpu(value);
 		// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value);
 		if (!knownComposite && isPrime)
@@ -1813,7 +1813,7 @@ internal static partial class PrimeOrderCalculator
 		// Console.WriteLine($"Factor completely after PollardRho hits {Volatile.Read(ref _factorCompletelyPollardRhoHits)}");
 
 		// isPrime = primeTester.HeuristicIsPrimeGpu(factor);
-		isPrime = tester.IsPrimeGpu(factor);
+		isPrime = PrimeTester.IsPrime(factor);
 		// isPrime = PrimeTester.IsPrimeGpu(factor);
 		// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor);
 		if (isPrime)
@@ -1850,7 +1850,7 @@ internal static partial class PrimeOrderCalculator
 			// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
 			// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
 
-			isPrime = tester.IsPrimeGpu(quotient);
+			isPrime = PrimeTester.IsPrime(quotient);
 			// isPrime = PrimeTester.IsPrimeGpu(quotient);
 			// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient);
 		}

PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs

@@ -950,7 +950,7 @@ internal static partial class PrimeOrderCalculator
         if (value <= ulong.MaxValue)
         {
             // return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None);
-            return PrimeTester.IsPrimeCpu((ulong)value, CancellationToken.None);
+            return PrimeTester.IsPrime((ulong)value);
         }
 
         if ((value & UInt128.One) == UInt128.Zero)

PerfectNumbers.Core/PrimeTester.cs

@@ -2,6 +2,7 @@ using System.Buffers;
 using System.Collections.Concurrent;
 using System.Numerics;
 using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
 using ILGPU;
 using ILGPU.Runtime;
 using PerfectNumbers.Core.Gpu;
@@ -19,7 +20,7 @@ public sealed class PrimeTester
 
 	public static PrimeTester Exclusive => _tester ??= new();
 
-	public static bool IsPrime(ulong n, CancellationToken ct)
+	public static bool IsPrime(ulong n)
 	{
 		if (n <= 1UL)
 		{
@@ -89,56 +90,70 @@ public sealed class PrimeTester
 		return result;
 	}
 
-	public static bool IsPrimeGpu(ulong n)
-	{
-		return Exclusive.IsPrimeGpu(n, CancellationToken.None);
-	}
-
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public bool IsPrimeGpu(ulong n, CancellationToken ct)
+	public bool IsPrimeGpu(ulong n)
 	{
 		bool forceCpu = GpuContextPool.ForceCpu;
 		bool belowGpuRange = n < 31UL;
 
 		if (forceCpu || belowGpuRange)
 		{
-			return IsPrimeCpu(n, ct);
+			return IsPrime(n);
 		}
 
-		Span<ulong> one = stackalloc ulong[1];
-		Span<byte> outFlags = stackalloc byte[1];
-		one[0] = n;
-		outFlags[0] = 0;
+		var limiter = GpuPrimeWorkLimiter.Acquire();
+		var gpu = PrimeTesterGpuContextPool.Rent(1);
+		var accelerator = gpu.Accelerator;
+		var kernel = gpu.Kernel;
 
-		IsPrimeBatchGpu(one, outFlags);
+		ulong value = n;
+		byte flag = 0;
 
-		return outFlags[0] != 0;
+		var input = gpu.Input;
+		var output = gpu.Output;
+		var inputView = input.View;
+		var outputView = output.View;
+
+		inputView.CopyFromCPU(ref value, 1);
+		kernel(
+						1,
+						inputView,
+						gpu.DevicePrimesDefault.View,
+						gpu.DevicePrimesLastOne.View,
+						gpu.DevicePrimesLastSeven.View,
+						gpu.DevicePrimesLastThree.View,
+						gpu.DevicePrimesLastNine.View,
+						gpu.DevicePrimesPow2Default.View,
+						gpu.DevicePrimesPow2LastOne.View,
+						gpu.DevicePrimesPow2LastSeven.View,
+						gpu.DevicePrimesPow2LastThree.View,
+						gpu.DevicePrimesPow2LastNine.View,
+						outputView);
+		accelerator.Synchronize();
+		outputView.CopyToCPU(ref flag, 1);
+
+		gpu.Dispose();
+		limiter.Dispose();
+
+		return flag != 0;
 	}
 
-	public static bool IsPrimeCpu(ulong n, CancellationToken ct)
-	{
-		// Preserve the legacy entry point for callers that bypass the thread-local caches.
-		return IsPrime(n, ct);
-	}
-
-
 	public static int GpuBatchSize { get; set; } = 262_144;
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	public static void IsPrimeBatchGpu(ReadOnlySpan<ulong> values, Span<byte> results)
 	{
 		var limiter = GpuPrimeWorkLimiter.Acquire();
-		var gpu = PrimeTesterGpuContextPool.Rent();
-		var state = gpu.State;
+		var gpu = PrimeTesterGpuContextPool.Rent(GpuBatchSize);
 		var accelerator = gpu.Accelerator;
+		var kernel = gpu.Kernel;
 		int totalLength = values.Length;
-		int batchSize = Math.Max(1, GpuBatchSize);
+		int batchSize = GpuBatchSize;
 
-		var scratch = state.RentScratch(batchSize, accelerator);
-		var input = scratch.Input;
-		var output = scratch.Output;
-		ArrayPool<ulong> pool = ThreadStaticPools.UlongPool;
-		ulong[] temp = pool.Rent(batchSize);
+		var input = gpu.Input;
+		var output = gpu.Output;
+		var inputView = input.View;
+		var outputView = output.View;
 
 		int pos = 0;
 		while (pos < totalLength)
@@ -146,290 +161,230 @@ public sealed class PrimeTester
 			int remaining = totalLength - pos;
 			int count = remaining > batchSize ? batchSize : remaining;
 
-			values.Slice(pos, count).CopyTo(temp);
-			input.View.CopyFromCPU(ref temp[0], count);
+			var valueSlice = values.Slice(pos, count);
+			ref ulong valueRef = ref MemoryMarshal.GetReference(valueSlice);
+			inputView.CopyFromCPU(ref valueRef, count);
 
-			state.Kernel(
-				count,
-				input.View,
-				state.DevicePrimesDefault.View,
-				state.DevicePrimesLastOne.View,
-				state.DevicePrimesLastSeven.View,
-				state.DevicePrimesLastThree.View,
-				state.DevicePrimesLastNine.View,
-				state.DevicePrimesPow2Default.View,
-				state.DevicePrimesPow2LastOne.View,
-				state.DevicePrimesPow2LastSeven.View,
-				state.DevicePrimesPow2LastThree.View,
-				state.DevicePrimesPow2LastNine.View,
-				output.View);
+			kernel(
+					count,
+					inputView,
+					gpu.DevicePrimesDefault.View,
+					gpu.DevicePrimesLastOne.View,
+					gpu.DevicePrimesLastSeven.View,
+					gpu.DevicePrimesLastThree.View,
+					gpu.DevicePrimesLastNine.View,
+					gpu.DevicePrimesPow2Default.View,
+					gpu.DevicePrimesPow2LastOne.View,
+					gpu.DevicePrimesPow2LastSeven.View,
+					gpu.DevicePrimesPow2LastThree.View,
+					gpu.DevicePrimesPow2LastNine.View,
+					outputView);
 			accelerator.Synchronize();
-			output.View.CopyToCPU(ref results[pos], count);
+			var resultSlice = results.Slice(pos, count);
+			ref byte resultRef = ref MemoryMarshal.GetReference(resultSlice);
+			outputView.CopyToCPU(ref resultRef, count);
 
 			pos += count;
 		}
 
-		pool.Return(temp, clearArray: false);
-		state.ReturnScratch(scratch);
-
 		gpu.Dispose();
 		limiter.Dispose();
 	}
 
 	internal static class PrimeTesterGpuContextPool
 	{
-		internal sealed class PooledContext
+		private static readonly ConcurrentQueue<PrimeTesterGpuContextLease> Pool = new();
+
+		internal static PrimeTesterGpuContextLease Rent(int minBufferCapacity = 1)
 		{
-			public Context Context { get; }
-
-			public Accelerator Accelerator { get; }
-
-			private KernelState? _kernelState;
-
-			public PooledContext()
+			if (!Pool.TryDequeue(out var lease))
 			{
-				Context = Context.CreateDefault();
-				Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context);
-				_kernelState = GpuKernelState.GetOrCreate(Accelerator);
+				lease = new PrimeTesterGpuContextLease(minBufferCapacity);
+			}
+			else
+			{
+				lease.EnsureCapacity(minBufferCapacity);
 			}
 
-			public KernelState KernelState
-			{
-				get
-				{
-					var state = _kernelState;
-					if (state is null || state.IsDisposed)
-					{
-						state = GpuKernelState.GetOrCreate(Accelerator);
-						_kernelState = state;
-					}
-
-					return state;
-				}
-			}
-
-			public void Dispose()
-			{
-				ClearGpuCaches(Accelerator);
-				_kernelState = null;
-				Accelerator.Dispose();
-				Context.Dispose();
-			}
+			lease.ResetReturnFlag();
+			return lease;
 		}
 
-		private static readonly ConcurrentQueue<PooledContext> Pool = new();
-
-		internal static PrimeTesterGpuContextLease Rent()
+		internal static void Return(PrimeTesterGpuContextLease lease)
 		{
-			if (Pool.TryDequeue(out var ctx))
-			{
-				return new PrimeTesterGpuContextLease(ctx);
-			}
-
-			return new PrimeTesterGpuContextLease(new PooledContext());
+			lease.Accelerator.Synchronize();
+			Pool.Enqueue(lease);
 		}
 
 		internal static void DisposeAll()
 		{
-			while (Pool.TryDequeue(out var ctx))
+			while (Pool.TryDequeue(out var lease))
 			{
-				ctx.Dispose();
+				lease.DisposeResources();
 			}
 		}
 
-		private static void Return(PooledContext ctx)
+		internal static void Clear(Accelerator accelerator)
 		{
-			ctx.Accelerator.Synchronize();
-			Pool.Enqueue(ctx);
-		}
-
-		internal readonly struct PrimeTesterGpuContextLease(PrimeTesterGpuContextPool.PooledContext ctx)
-		{
-			private readonly PooledContext? _ctx = ctx;
-
-			public readonly Accelerator Accelerator => _ctx!.Accelerator;
-
-
-			public readonly KernelState State => _ctx!.KernelState;
-
-			public readonly void Dispose()
-			{
-				Return(_ctx ?? throw new InvalidOperationException("GPU context lease is not initialized."));
-			}
-		}
-	}
-
-	// Per-accelerator GPU state for prime sieve (kernel + uploaded primes).
-	internal sealed class KernelState
-	{
-		public Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> Kernel;
-		public Action<Index1D, ArrayView<ulong>, ulong, ArrayView<byte>> HeuristicTrialDivisionKernel;
-		public MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesDefault;
-		public MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastOne;
-		public MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastSeven;
-		public MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastThree;
-		public MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastNine;
-		public MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2Default;
-		public MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastOne;
-		public MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastSeven;
-		public MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastThree;
-		public MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastNine;
-		public bool IsDisposed;
-		private readonly ConcurrentBag<ScratchBuffers> _scratchPool = [];
-		// TODO: Replace this ConcurrentBag with the lock-free ring buffer variant validated in
-		// GpuModularArithmeticBenchmarks so renting scratch buffers stops contending on the bag's internal locks when
-		// thousands of GPU batches execute per second.
-
-		public KernelState(Accelerator accelerator)
-		{
-			IsDisposed = false;
-			// Compile once per accelerator and upload primes once.
-			Kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterKernels.SmallPrimeSieveKernel);
-			HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ulong, ArrayView<byte>>(PrimeTesterKernels.HeuristicTrialDivisionKernel);
-
-			var primesDefault = DivisorGenerator.SmallPrimes;
-			DevicePrimesDefault = accelerator.Allocate1D<uint>(primesDefault.Length);
-			DevicePrimesDefault.View.CopyFromCPU(primesDefault);
-
-			var primesLastOne = DivisorGenerator.SmallPrimesLastOne;
-			DevicePrimesLastOne = accelerator.Allocate1D<uint>(primesLastOne.Length);
-			DevicePrimesLastOne.View.CopyFromCPU(primesLastOne);
-
-			var primesLastSeven = DivisorGenerator.SmallPrimesLastSeven;
-			DevicePrimesLastSeven = accelerator.Allocate1D<uint>(primesLastSeven.Length);
-			DevicePrimesLastSeven.View.CopyFromCPU(primesLastSeven);
-
-			var primesLastThree = DivisorGenerator.SmallPrimesLastThree;
-			DevicePrimesLastThree = accelerator.Allocate1D<uint>(primesLastThree.Length);
-			DevicePrimesLastThree.View.CopyFromCPU(primesLastThree);
-
-			var primesLastNine = DivisorGenerator.SmallPrimesLastNine;
-			DevicePrimesLastNine = accelerator.Allocate1D<uint>(primesLastNine.Length);
-			DevicePrimesLastNine.View.CopyFromCPU(primesLastNine);
-
-			var primesPow2Default = DivisorGenerator.SmallPrimesPow2;
-			DevicePrimesPow2Default = accelerator.Allocate1D<ulong>(primesPow2Default.Length);
-			DevicePrimesPow2Default.View.CopyFromCPU(primesPow2Default);
-
-			var primesPow2LastOne = DivisorGenerator.SmallPrimesPow2LastOne;
-			DevicePrimesPow2LastOne = accelerator.Allocate1D<ulong>(primesPow2LastOne.Length);
-			DevicePrimesPow2LastOne.View.CopyFromCPU(primesPow2LastOne);
-
-			var primesPow2LastSeven = DivisorGenerator.SmallPrimesPow2LastSeven;
-			DevicePrimesPow2LastSeven = accelerator.Allocate1D<ulong>(primesPow2LastSeven.Length);
-			DevicePrimesPow2LastSeven.View.CopyFromCPU(primesPow2LastSeven);
-
-			var primesPow2LastThree = DivisorGenerator.SmallPrimesPow2LastThree;
-			DevicePrimesPow2LastThree = accelerator.Allocate1D<ulong>(primesPow2LastThree.Length);
-			DevicePrimesPow2LastThree.View.CopyFromCPU(primesPow2LastThree);
-
-			var primesPow2LastNine = DivisorGenerator.SmallPrimesPow2LastNine;
-			DevicePrimesPow2LastNine = accelerator.Allocate1D<ulong>(primesPow2LastNine.Length);
-			DevicePrimesPow2LastNine.View.CopyFromCPU(primesPow2LastNine);
-		}
-
-		internal sealed class ScratchBuffers
-		{
-			public MemoryBuffer1D<ulong, Stride1D.Dense> Input { get; private set; }
-			public MemoryBuffer1D<byte, Stride1D.Dense> Output { get; private set; }
-			public int Capacity { get; private set; }
-
-			public ScratchBuffers(Accelerator accel, int capacity)
-			{
-				Capacity = Math.Max(1, capacity);
-				Input = accel.Allocate1D<ulong>(Capacity);
-				Output = accel.Allocate1D<byte>(Capacity);
-			}
-
-			public void Dispose()
-			{
-				Output.Dispose();
-				Input.Dispose();
-			}
-		}
-
-		[MethodImpl(MethodImplOptions.AggressiveInlining)]
-		public ScratchBuffers RentScratch(int minCapacity, Accelerator accel)
-		{
-			while (_scratchPool.TryTake(out var sb))
-			{
-				if (sb.Capacity >= minCapacity)
-				{
-					return sb;
-				}
-
-				sb.Dispose();
-			}
-
-			return new ScratchBuffers(accel, minCapacity);
-		}
-
-		[MethodImpl(MethodImplOptions.AggressiveInlining)]
-		public void ReturnScratch(ScratchBuffers scratch)
-		{
-			_scratchPool.Add(scratch);
-		}
-
-		public void Clear()
-		{
-			if (IsDisposed)
+			if (accelerator is null)
 			{
 				return;
 			}
 
-			while (_scratchPool.TryTake(out var sb))
+			var retained = new List<PrimeTesterGpuContextLease>();
+			while (Pool.TryDequeue(out var lease))
 			{
-				sb.Dispose();
-			}
-
-			DevicePrimesDefault.Dispose();
-			DevicePrimesLastOne.Dispose();
-			DevicePrimesLastSeven.Dispose();
-			DevicePrimesLastThree.Dispose();
-			DevicePrimesLastNine.Dispose();
-			DevicePrimesPow2Default.Dispose();
-			DevicePrimesPow2LastOne.Dispose();
-			DevicePrimesPow2LastSeven.Dispose();
-			DevicePrimesPow2LastThree.Dispose();
-			DevicePrimesPow2LastNine.Dispose();
-			IsDisposed = true;
-		}
-	}
-
-	internal static class GpuKernelState
-	{
-		// Map accelerator to cached state; use Lazy to serialize kernel creation
-		private static readonly System.Collections.Concurrent.ConcurrentDictionary<Accelerator, Lazy<KernelState>> States = new();
-		// TODO: Prewarm this per-accelerator cache during startup (and reuse a simple array keyed by accelerator index)
-		// once the kernel pool exposes deterministic ordering; the Lazy wrappers showed measurable overhead in the
-		// GpuModularArithmeticBenchmarks hot path.
-
-		[MethodImpl(MethodImplOptions.AggressiveInlining)]
-		public static KernelState GetOrCreate(Accelerator accelerator)
-		{
-			var lazy = States.GetOrAdd(
-				accelerator,
-				acc => new Lazy<KernelState>(() => new KernelState(acc), LazyThreadSafetyMode.ExecutionAndPublication));
-			return lazy.Value;
-		}
-
-		public static void Clear(Accelerator accelerator)
-		{
-			if (States.TryRemove(accelerator, out var lazy))
-			{
-				if (lazy.IsValueCreated)
+				if (lease.Accelerator == accelerator)
 				{
-					var state = lazy.Value;
-					state.Clear();
+					lease.DisposeResources();
+				}
+				else
+				{
+					retained.Add(lease);
 				}
 			}
+
+			int retainedCount = retained.Count;
+			for (int i = 0; i < retainedCount; i++)
+			{
+				Pool.Enqueue(retained[i]);
+			}
+		}
+
+		internal sealed class PrimeTesterGpuContextLease : IDisposable
+		{
+			private bool _returned;
+
+			internal readonly Context AcceleratorContext;
+			public readonly Accelerator Accelerator;
+			public readonly Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> Kernel;
+			public readonly Action<Index1D, ArrayView<ulong>, ulong, ArrayView<int>> HeuristicTrialDivisionKernel;
+			public readonly MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesDefault;
+			public readonly MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastOne;
+			public readonly MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastSeven;
+			public readonly MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastThree;
+			public readonly MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastNine;
+			public readonly MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2Default;
+			public readonly MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastOne;
+			public readonly MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastSeven;
+			public readonly MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastThree;
+			public readonly MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimesPow2LastNine;
+			public MemoryBuffer1D<ulong, Stride1D.Dense> Input;
+			public MemoryBuffer1D<byte, Stride1D.Dense> Output;
+			public MemoryBuffer1D<int, Stride1D.Dense> HeuristicFlag;
+			public int BufferCapacity;
+
+			internal PrimeTesterGpuContextLease(int minBufferCapacity)
+			{
+				AcceleratorContext = Context.CreateDefault();
+				Accelerator = AcceleratorContext.GetPreferredDevice(false).CreateAccelerator(AcceleratorContext);
+				Kernel = Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterKernels.SmallPrimeSieveKernel);
+				HeuristicTrialDivisionKernel = Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ulong, ArrayView<int>>(PrimeTesterKernels.HeuristicTrialDivisionKernel);
+
+				var primesDefault = DivisorGenerator.SmallPrimes;
+				DevicePrimesDefault = Accelerator.Allocate1D<uint>(primesDefault.Length);
+				DevicePrimesDefault.View.CopyFromCPU(primesDefault);
+
+				var primesLastOne = DivisorGenerator.SmallPrimesLastOne;
+				DevicePrimesLastOne = Accelerator.Allocate1D<uint>(primesLastOne.Length);
+				DevicePrimesLastOne.View.CopyFromCPU(primesLastOne);
+
+				var primesLastSeven = DivisorGenerator.SmallPrimesLastSeven;
+				DevicePrimesLastSeven = Accelerator.Allocate1D<uint>(primesLastSeven.Length);
+				DevicePrimesLastSeven.View.CopyFromCPU(primesLastSeven);
+
+				var primesLastThree = DivisorGenerator.SmallPrimesLastThree;
+				DevicePrimesLastThree = Accelerator.Allocate1D<uint>(primesLastThree.Length);
+				DevicePrimesLastThree.View.CopyFromCPU(primesLastThree);
+
+				var primesLastNine = DivisorGenerator.SmallPrimesLastNine;
+				DevicePrimesLastNine = Accelerator.Allocate1D<uint>(primesLastNine.Length);
+				DevicePrimesLastNine.View.CopyFromCPU(primesLastNine);
+
+				var primesPow2Default = DivisorGenerator.SmallPrimesPow2;
+				DevicePrimesPow2Default = Accelerator.Allocate1D<ulong>(primesPow2Default.Length);
+				DevicePrimesPow2Default.View.CopyFromCPU(primesPow2Default);
+
+				var primesPow2LastOne = DivisorGenerator.SmallPrimesPow2LastOne;
+				DevicePrimesPow2LastOne = Accelerator.Allocate1D<ulong>(primesPow2LastOne.Length);
+				DevicePrimesPow2LastOne.View.CopyFromCPU(primesPow2LastOne);
+
+				var primesPow2LastSeven = DivisorGenerator.SmallPrimesPow2LastSeven;
+				DevicePrimesPow2LastSeven = Accelerator.Allocate1D<ulong>(primesPow2LastSeven.Length);
+				DevicePrimesPow2LastSeven.View.CopyFromCPU(primesPow2LastSeven);
+
+				var primesPow2LastThree = DivisorGenerator.SmallPrimesPow2LastThree;
+				DevicePrimesPow2LastThree = Accelerator.Allocate1D<ulong>(primesPow2LastThree.Length);
+				DevicePrimesPow2LastThree.View.CopyFromCPU(primesPow2LastThree);
+
+				var primesPow2LastNine = DivisorGenerator.SmallPrimesPow2LastNine;
+				DevicePrimesPow2LastNine = Accelerator.Allocate1D<ulong>(primesPow2LastNine.Length);
+				DevicePrimesPow2LastNine.View.CopyFromCPU(primesPow2LastNine);
+
+				BufferCapacity = minBufferCapacity < 1 ? 1 : minBufferCapacity;
+				Input = Accelerator.Allocate1D<ulong>(BufferCapacity);
+				Output = Accelerator.Allocate1D<byte>(BufferCapacity);
+				HeuristicFlag = Accelerator.Allocate1D<int>(BufferCapacity);
+				_returned = false;
+			}
+
+			public void EnsureCapacity(int minCapacity)
+			{
+				int required = minCapacity < 1 ? 1 : minCapacity;
+				if (required <= BufferCapacity)
+				{
+					return;
+				}
+
+				Input.Dispose();
+				Output.Dispose();
+				HeuristicFlag.Dispose();
+				Input = Accelerator.Allocate1D<ulong>(required);
+				Output = Accelerator.Allocate1D<byte>(required);
+				HeuristicFlag = Accelerator.Allocate1D<int>(required);
+				BufferCapacity = required;
+			}
+
+			internal void ResetReturnFlag()
+			{
+				_returned = false;
+			}
+
+			public void Dispose()
+			{
+				if (_returned)
+				{
+					return;
+				}
+
+				_returned = true;
+				PrimeTesterGpuContextPool.Return(this);
+			}
+
+			internal void DisposeResources()
+			{
+				Input.Dispose();
+				Output.Dispose();
+				HeuristicFlag.Dispose();
+				DevicePrimesDefault.Dispose();
+				DevicePrimesLastOne.Dispose();
+				DevicePrimesLastSeven.Dispose();
+				DevicePrimesLastThree.Dispose();
+				DevicePrimesLastNine.Dispose();
+				DevicePrimesPow2Default.Dispose();
+				DevicePrimesPow2LastOne.Dispose();
+				DevicePrimesPow2LastSeven.Dispose();
+				DevicePrimesPow2LastThree.Dispose();
+				DevicePrimesPow2LastNine.Dispose();
+				Accelerator.Dispose();
+				AcceleratorContext.Dispose();
+			}
 		}
 	}
-
 	// Expose cache clearing for accelerator disposal coordination
 	public static void ClearGpuCaches(Accelerator accelerator)
 	{
-		GpuKernelState.Clear(accelerator);
+		PrimeTesterGpuContextPool.Clear(accelerator);
 	}
 
 	internal static void DisposeGpuContexts()

docs/HeuristicIsPrimePlan.md

@@ -131,7 +131,7 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels,
    * [done] The CPU heuristic now operates as a pure trial-division sweep without maintaining per-call summaries, relying on `PrepareHeuristicDivisor` only when downstream consumers request Montgomery data.
    * [done] Introduced the temporary `UseHeuristicGroupBTrialDivision` gate so current builds execute Group A locally and then fall back to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, keeping the full Group B implementation available for reactivation.
 7. [done] **Implement `IsPrimeGpu`**
-   * [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `sqrtLimit`/`nMod10` inputs, and falling back to CPU when GPUs are disabled.
+   * [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `maxDivisorSquare`/`nMod10` inputs, and falling back to CPU when GPUs are disabled.
    * [done] Added `PrimeTesterKernels.HeuristicTrialDivisionKernel` so GPU batches respect Group A/B ordering while emitting hit flags for early exits.
    * [done] Kept GPU divisibility in standard modular arithmetic while reserving Montgomery transforms for CPU confirmation logic.
    * [done] When the temporary Group B gate is disabled, delegate GPU calls to the CPU fallback so Group A coverage remains in place while the accelerator path stays available for future re-enablement.