diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs
new file mode 100644
index 0000000..e1e9f9c
--- /dev/null
+++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs
@@ -0,0 +1,246 @@
+using System;
+using System.Collections.Generic;
+using System.Threading;
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Jobs;
+using Open.Numeric.Primes;
+using PerfectNumbers.Core;
+
+namespace EvenPerfectBitScanner.Benchmarks;
+
+[ShortRunJob(RuntimeMoniker.Net80)]
+[MemoryDiagnoser]
+public class PrimeTesterBenchmarks
+{
+    private PrimeTester _cpuTester = null!;
+    private PrimeTester _gpuTester = null!;
+    private ulong[] _candidates = [];
+
+    [ParamsSource(nameof(GetBenchmarkCases))]
+    public PrimeBenchmarkCase BenchmarkCase { get; set; }
+
+    [GlobalSetup]
+    public void GlobalSetup()
+    {
+        PrimeTester.EnableHeuristicPrimeTesting = true;
+        _cpuTester = new PrimeTester();
+        _gpuTester = new PrimeTester();
+        _candidates = BenchmarkCase.Candidates;
+
+        if (_candidates.Length == 0)
+        {
+            return;
+        }
+
+        CancellationToken cancellationToken = CancellationToken.None;
+        _cpuTester.IsPrime(_candidates[0], cancellationToken);
+        _gpuTester.IsPrimeGpu(_candidates[0], cancellationToken);
+        _ = Prime.Numbers.IsPrime(_candidates[0]);
+    }
+
+    [Benchmark(Baseline = true)]
+    public int HeuristicCpu()
+    {
+        PrimeTester tester = _cpuTester;
+        ulong[] values = _candidates;
+        CancellationToken cancellationToken = CancellationToken.None;
+        int primeCount = 0;
+
+        for (int i = 0; i < values.Length; i++)
+        {
+            if (tester.IsPrime(values[i], cancellationToken))
+            {
+                primeCount++;
+            }
+        }
+
+        return primeCount;
+    }
+
+    [Benchmark]
+    public int HeuristicGpu()
+    {
+        PrimeTester tester = _gpuTester;
+        ulong[] values = _candidates;
+        CancellationToken cancellationToken = CancellationToken.None;
+        int primeCount = 0;
+
+        for (int i = 0; i < values.Length; i++)
+        {
+            if (tester.HeuristicIsPrimeGpu(values[i]))
+            {
+                primeCount++;
+            }
+        }
+
+        return primeCount;
+    }
+
+    // This benchmark is very slow. No point in comparing it.
+    // [Benchmark]
+    public int NonHeuristicGpu()
+    {
+        // PrimeTester tester = _gpuTester;
+        ulong[] values = _candidates;
+        CancellationToken cancellationToken = CancellationToken.None;
+        int primeCount = 0;
+
+        for (int i = 0; i < values.Length; i++)
+        {
+            if (PrimeTester.LegacyIsPrimeGpu(values[i], cancellationToken))
+            {
+                primeCount++;
+            }
+        }
+
+        return primeCount;
+    }
+
+    [Benchmark]
+    public int NonHeuristicCpu()
+    {
+        // PrimeTester tester = _gpuTester;
+        ulong[] values = _candidates;
+        CancellationToken cancellationToken = CancellationToken.None;
+        int primeCount = 0;
+
+        for (int i = 0; i < values.Length; i++)
+        {
+            if (PrimeTester.LegacyIsPrimeInternal(values[i], cancellationToken))
+            {
+                primeCount++;
+            }
+        }
+
+        return primeCount;
+    }
+
+    [Benchmark]
+    public int OpenNumericPrimes()
+    {
+        ulong[] values = _candidates;
+        int primeCount = 0;
+
+        for (int i = 0; i < values.Length; i++)
+        {
+            if (Prime.Numbers.IsPrime(values[i]))
+            {
+                primeCount++;
+            }
+        }
+
+        return primeCount;
+    }
+
+    public static IEnumerable<PrimeBenchmarkCase> GetBenchmarkCases()
+    {
+        yield return PrimeBenchmarkCase.FromRange(
+            "≤100",
+            "Admissible odd values ≤ 100 (k ≤ 5, not divisible by 5)",
+            3UL,
+            100UL,
+            int.MaxValue);
+
+        yield return PrimeBenchmarkCase.FromRange(
+            "≤4_000_000",
+            "Sampled admissible odd values ≤ 4,000,000 (k ≤ 5, not divisible by 5)",
+            3UL,
+            4_000_000UL,
+            4_096);
+
+        yield return PrimeBenchmarkCase.FromRange(
+            "≥138_000_000",
+            "Sampled admissible odd values ≥ 138,000,000 (k ≤ 5, not divisible by 5)",
+            138_000_001UL,
+            138_500_000UL,
+            4_096);
+    }
+
+    public readonly record struct PrimeBenchmarkCase(string Name, string Description, ulong[] Candidates)
+    {
+        public override string ToString()
+        {
+            return Name;
+        }
+
+        public static PrimeBenchmarkCase FromRange(
+            string name,
+            string description,
+            ulong startInclusive,
+            ulong endInclusive,
+            int maxCount)
+        {
+            ulong[] values = CollectRange(startInclusive, endInclusive, maxCount);
+            return new PrimeBenchmarkCase(name, description, values);
+        }
+    }
+
+    private static ulong[] CollectRange(ulong startInclusive, ulong endInclusive, int maxCount)
+    {
+        maxCount = maxCount < 0 ? 0 : maxCount;
+        if (maxCount == 0)
+        {
+            return [];
+        }
+
+        int estimatedLength = (int)Math.Min((ulong)maxCount, endInclusive - startInclusive + 1UL);
+        ulong[] buffer = new ulong[estimatedLength];
+        int count = 0;
+
+        for (ulong value = startInclusive; value <= endInclusive && count < maxCount; value++)
+        {
+            if (!IsAdmissibleMersenneDivisorCandidate(value))
+            {
+                continue;
+            }
+
+            if (count == buffer.Length)
+            {
+                Array.Resize(ref buffer, buffer.Length * 2);
+            }
+
+            buffer[count++] = value;
+        }
+
+        if (count == buffer.Length)
+        {
+            return buffer;
+        }
+
+        ulong[] result = new ulong[count];
+        Array.Copy(buffer, result, count);
+        return result;
+    }
+
+    private static bool IsAdmissibleMersenneDivisorCandidate(ulong value)
+    {
+        if ((value & 1UL) == 0UL || value % 5UL == 0UL || value <= 1UL)
+        {
+            return false;
+        }
+
+        ulong qMinusOne = value - 1UL;
+
+        for (ulong k = 1; k <= 5; k++)
+        {
+            ulong denominator = k << 1;
+            if (qMinusOne % denominator != 0UL)
+            {
+                continue;
+            }
+
+            ulong p = qMinusOne / denominator;
+            if (p < 2UL)
+            {
+                continue;
+            }
+
+            if (Prime.Numbers.IsPrime(p))
+            {
+                return true;
+            }
+        }
+
+        return false;
+    }
+}
diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
index 5ba86d0..e4f4924 100644
--- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
+++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs
@@ -12,16 +12,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
     private ulong _lastStatusDivisor;
     private int _batchSize = 1_024;
 
-    [ThreadStatic]
-    private static GpuUInt128WorkSet _divisorScanGpuWorkSet;
-
-    private struct GpuUInt128WorkSet
-    {
-        public GpuUInt128 Step;
-        public GpuUInt128 Divisor;
-        public GpuUInt128 Limit;
-    }
-
     public int BatchSize
     {
         get => _batchSize;
@@ -113,151 +103,37 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         out bool processedAll,
         out ulong processedCount)
     {
-        processedCount = 0UL;
         processedAll = false;
+        processedCount = 0UL;
 
-        // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs.
-        // Keeping the guard commented out documents the reasoning for benchmarks and tests.
-        // if (allowedMax < 3UL)
-        // {
-        //     return false;
-        // }
+        var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
 
-        ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet;
-
-        ref GpuUInt128 step = ref workSet.Step;
-        step.High = 0UL;
-        step.Low = prime;
-        step.ShiftLeft(1);
-
-        ref GpuUInt128 limit = ref workSet.Limit;
-        limit.High = 0UL;
-        limit.Low = allowedMax;
-
-        ref GpuUInt128 divisor = ref workSet.Divisor;
-        divisor.High = step.High;
-        divisor.Low = step.Low;
-        divisor.Add(1UL);
-        if (divisor.CompareTo(limit) > 0)
+        while (enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
         {
-            processedAll = true;
-            divisor.High = 0UL;
-            divisor.Low = 0UL;
-            step.High = 0UL;
-            step.Low = 0UL;
-            limit.High = 0UL;
-            limit.Low = 0UL;
-            return false;
-        }
+            processedCount = enumerator.ProcessedCount;
 
-        // Intentionally recomputes factorizations without a per-thread cache.
-        // The previous factor cache recorded virtually no hits and only slowed down the scan.
-        DivisorCycleCache cycleCache = DivisorCycleCache.Shared;
+            PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate);
+            ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength(
+                prime,
+                in preparation,
+                out _,
+                out _,
+                out _);
 
-        ulong stepHigh = step.High;
-        ulong stepLow = step.Low;
-        byte step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL);
-        byte step8 = (byte)(stepLow % 8UL);
-        byte step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL);
-        byte step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL);
-        byte step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL);
-        byte step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL);
-
-        ulong divisorHigh = divisor.High;
-        ulong divisorLow = divisor.Low;
-        byte remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL);
-        byte remainder8 = (byte)(divisorLow % 8UL);
-        byte remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL);
-        byte remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL);
-        byte remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL);
-        byte remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL);
-
-        // Keep the divisibility filters aligned with the divisor-cycle generator so the
-        // CPU path never requests cycles that were skipped during cache creation.
-        bool lastIsSeven = (prime & 3UL) == 3UL;
-        const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
-        const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
-
-        while (divisor.CompareTo(limit) <= 0)
-        {
-            ulong candidate = divisor.Low;
-            processedCount++;
-
-            ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
-            bool admissible = ((decimalMask >> remainder10) & 1) != 0;
-
-            if (admissible
-                && (remainder8 == 1 || remainder8 == 7)
-                && remainder3 != 0
-                && remainder5 != 0
-                && remainder7 != 0
-                && remainder11 != 0)
+            if (divisorCycle == prime)
             {
-                MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
-                ulong divisorCycle;
-                // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs.
-                // if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
-                // {
-                //     divisorCycle = cycleCache.GetCycleLength(candidate);
-                // }
-                // else
-                {
-                    if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
-                            candidate,
-                            prime,
-                            divisorData,
-                            out ulong computedCycle,
-                            out bool primeOrderFailed) || computedCycle == 0UL)
-                    {
-                        // Divisors produced by 2 * k * p + 1 always exceed PerfectNumberConstants.MaxQForDivisorCycles
-                        // for the exponents scanned here, so skip the unused cache fallback and compute directly.
-                        divisorCycle = MersenneDivisorCycles.CalculateCycleLength(
-                            candidate,
-                            divisorData,
-                            skipPrimeOrderHeuristic: primeOrderFailed);
-                    }
-                    else
-                    {
-                        divisorCycle = computedCycle;
-                    }
-                }
-
-                if (divisorCycle == prime)
-                {
-                    // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides
-                    // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p.
-                    divisor.High = 0UL;
-                    divisor.Low = 0UL;
-                    step.High = 0UL;
-                    step.Low = 0UL;
-                    limit.High = 0UL;
-                    limit.Low = 0UL;
-                    processedAll = true;
-                    return true;
-                }
-
-                if (divisorCycle == 0UL)
-                {
-                    Console.WriteLine($"Divisor cycle was not calculated for {prime}");
-                }
+                processedAll = true;
+                return true;
             }
 
-            divisor.Add(step);
-            remainder10 = AddMod(remainder10, step10, (byte)10);
-            remainder8 = AddMod(remainder8, step8, (byte)8);
-            remainder5 = AddMod(remainder5, step5, (byte)5);
-            remainder3 = AddMod(remainder3, step3, (byte)3);
-            remainder7 = AddMod(remainder7, step7, (byte)7);
-            remainder11 = AddMod(remainder11, step11, (byte)11);
+            if (divisorCycle == 0UL)
+            {
+                Console.WriteLine($"Divisor cycle was not calculated for {prime}");
+            }
         }
 
-        processedAll = divisor.CompareTo(limit) > 0;
-        divisor.High = 0UL;
-        divisor.Low = 0UL;
-        step.High = 0UL;
-        step.Low = 0UL;
-        limit.High = 0UL;
-        limit.Low = 0UL;
+        processedCount = enumerator.ProcessedCount;
+        processedAll = enumerator.Exhausted;
         return false;
     }
 
@@ -267,17 +143,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         return residue == 1UL ? (byte)1 : (byte)0;
     }
 
-    private static byte AddMod(byte value, byte delta, byte modulus)
-    {
-        int sum = value + delta;
-        if (sum >= modulus)
-        {
-            sum -= modulus;
-        }
-
-        return (byte)sum;
-    }
-
     private void UpdateStatusUnsafe(ulong processedCount)
     {
         if (processedCount == 0UL)
diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs
index 1897d8f..1d284e6 100644
--- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs
+++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs
@@ -41,6 +41,17 @@ internal static class PrimeTesterKernels
         results[index] = composite ? (byte)0 : (byte)1;
     }
 
+    public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView<ulong> divisors, ulong n, ArrayView<byte> results)
+    {
+        ulong divisor = divisors[index];
+        if (divisor <= 1UL)
+        {
+            results[index] = 0;
+            return;
+        }
+
+        results[index] = n % divisor == 0UL ? (byte)1 : (byte)0;
+    }
 
     public static void SharesFactorKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<byte> results)
     {
diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
index f363126..9e9b729 100644
--- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
@@ -112,7 +112,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             resources.HitIndexBuffer,
             resources.Divisors,
             resources.Exponents,
-            resources.FilteredDivisors,
             resources.DivisorData,
             resources.Offsets,
             resources.Counts,
@@ -177,7 +176,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         MemoryBuffer1D<int, Stride1D.Dense> hitIndexBuffer,
         ulong[] divisors,
         ulong[] exponents,
-        ulong[] filteredDivisors,
         GpuDivisorPartialData[] divisorData,
         int[] offsets,
         int[] counts,
@@ -188,131 +186,36 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
     {
         int batchCapacity = (int)divisorDataBuffer.Length;
         bool composite = false;
-        bool processedAll = false;
         processedCount = 0UL;
         lastProcessed = 0UL;
 
-        int chunkCapacity = Math.Max(1, batchCapacity);
+        var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
 
-        UInt128 twoP128 = (UInt128)prime << 1;
-        UInt128 allowedMax128 = allowedMax;
-        UInt128 firstDivisor128 = twoP128 + UInt128.One;
-        bool invalidStride = twoP128 == UInt128.Zero;
-        bool outOfRange = firstDivisor128 > allowedMax128;
-        UInt128 numerator = allowedMax128 - UInt128.One;
-        UInt128 maxK128 = invalidStride || outOfRange ? UInt128.Zero : numerator / twoP128;
-        bool hasCandidates = !invalidStride && !outOfRange && maxK128 != UInt128.Zero;
-        if (!hasCandidates)
+        Span<ulong> divisorSpan = divisors.AsSpan();
+        Span<ulong> exponentSpan = exponents.AsSpan();
+        Span<GpuDivisorPartialData> divisorDataSpan = divisorData.AsSpan();
+        Span<int> offsetSpan = offsets.AsSpan();
+        Span<int> countSpan = counts.AsSpan();
+        Span<ulong> cycleSpan = cycles.AsSpan();
+
+        while (!composite)
         {
-            coveredRange = true;
-            return false;
-        }
-
-        ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128;
-        ulong currentK = 1UL;
-
-        byte step10 = (byte)(twoP128 % 10UL);
-        byte step8 = (byte)(twoP128 % 8UL);
-        byte step5 = (byte)(twoP128 % 5UL);
-        byte step3 = (byte)(twoP128 % 3UL);
-        byte step7 = (byte)(twoP128 % 7UL);
-        byte step11 = (byte)(twoP128 % 11UL);
-
-        UInt128 currentDivisor128 = firstDivisor128;
-        byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL));
-        byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL));
-        byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL));
-        byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL));
-        byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL));
-        byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL));
-        bool lastIsSeven = (prime & 3UL) == 3UL;
-
-        Span<ulong> filteredStorage = filteredDivisors.AsSpan();
-        Span<ulong> divisorStorage = divisors.AsSpan();
-        Span<ulong> exponentStorage = exponents.AsSpan();
-        Span<GpuDivisorPartialData> divisorDataStorage = divisorData.AsSpan();
-        Span<int> offsetStorage = offsets.AsSpan();
-        Span<int> countStorage = counts.AsSpan();
-        Span<ulong> cycleStorage = cycles.AsSpan();
-        while (currentK <= maxK && !composite)
-        {
-            int chunkCount = Math.Min(chunkCapacity, batchCapacity);
-            ulong remainingK = maxK - currentK + 1UL;
-            if ((ulong)chunkCount > remainingK)
-            {
-                chunkCount = (int)remainingK;
-            }
-
-            if (chunkCount <= 0)
-            {
-                processedAll = true;
-                break;
-            }
-
-            int filteredCount = 0;
-            UInt128 nextDivisor128 = currentDivisor128;
-
-            byte localRemainder10 = remainder10;
-            byte localRemainder8 = remainder8;
-            byte localRemainder5 = remainder5;
-            byte localRemainder3 = remainder3;
-            byte localRemainder7 = remainder7;
-            byte localRemainder11 = remainder11;
-
-            for (int i = 0; i < chunkCount; i++)
-            {
-                ulong candidate = (ulong)nextDivisor128;
-
-                if (CandidatePassesHeuristics(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastIsSeven))
-                {
-                    filteredStorage[filteredCount++] = candidate;
-                }
-
-                nextDivisor128 += twoP128;
-                localRemainder10 = AddMod(localRemainder10, step10, 10);
-                localRemainder8 = AddMod(localRemainder8, step8, 8);
-                localRemainder5 = AddMod(localRemainder5, step5, 5);
-                localRemainder3 = AddMod(localRemainder3, step3, 3);
-                localRemainder7 = AddMod(localRemainder7, step7, 7);
-                localRemainder11 = AddMod(localRemainder11, step11, 11);
-            }
-
-            remainder10 = localRemainder10;
-            remainder8 = localRemainder8;
-            remainder5 = localRemainder5;
-            remainder3 = localRemainder3;
-            remainder7 = localRemainder7;
-            remainder11 = localRemainder11;
-
-            processedCount += (ulong)chunkCount;
-
-            UInt128 lastDivisor128 = nextDivisor128 - twoP128;
-            lastProcessed = (ulong)lastDivisor128;
-            currentDivisor128 = nextDivisor128;
-
-            currentK += (ulong)chunkCount;
-            processedAll = currentK > maxK;
-
-            if (filteredCount == 0)
-            {
-                continue;
-            }
-
-            Span<ulong> filteredDivisorsSpan = filteredStorage[..filteredCount];
-
-            Span<ulong> divisorSpan = divisorStorage;
-            Span<ulong> exponentSpan = exponentStorage;
-            Span<GpuDivisorPartialData> divisorDataSpan = divisorDataStorage;
-            Span<int> offsetSpan = offsetStorage;
-            Span<int> countSpan = countStorage;
-            Span<ulong> cycleSpan = cycleStorage;
             int admissibleCount = 0;
             int exponentIndex = 0;
-            for (int i = 0; i < filteredCount; i++)
+
+            while (admissibleCount < batchCapacity && enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
             {
-                ulong divisorValue = filteredDivisorsSpan[i];
-                MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue);
-                ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData);
+                processedCount = enumerator.ProcessedCount;
+                ulong divisorValue = candidate.Value;
+
+                PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate);
+                ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength(
+                    prime,
+                    in preparation,
+                    out _,
+                    out _,
+                    out _);
+
                 if (divisorCycle != prime)
                 {
                     continue;
@@ -329,8 +232,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                 exponentIndex++;
             }
 
+            processedCount = enumerator.ProcessedCount;
+            lastProcessed = enumerator.LastDivisor;
+
             if (admissibleCount == 0)
             {
+                if (enumerator.Exhausted)
+                {
+                    break;
+                }
+
                 continue;
             }
 
@@ -362,74 +273,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex];
         }
 
-        coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128);
+        coveredRange = composite || enumerator.Exhausted;
         return composite;
     }
 
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong ResolveDivisorCycle(ulong divisor, ulong prime, in MontgomeryDivisorData divisorData)
-    {
-        if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL)
-        {
-            return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
-        }
-
-        return computedCycle;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static bool CandidatePassesHeuristics(
-        byte remainder10,
-        byte remainder8,
-        byte remainder3,
-        byte remainder5,
-        byte remainder7,
-        byte remainder11,
-        bool lastIsSeven)
-    {
-        const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
-        const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
-
-        ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
-        if (((decimalMask >> remainder10) & 1) == 0)
-        {
-            return false;
-        }
-
-        if (remainder8 != 1 && remainder8 != 7)
-        {
-            return false;
-        }
-
-        if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0)
-        {
-            return false;
-        }
-
-        return true;
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static byte AddMod(byte value, byte delta, byte modulus)
-    {
-        int result = value + delta;
-        if (result >= modulus)
-        {
-            result -= modulus;
-            if (result >= modulus)
-            {
-                result -= modulus;
-            }
-        }
-
-        return (byte)result;
-    }
-
-
-
-
-
     private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime)
     {
         // The --mersenne=bydivisor flow in EvenPerfectBitScanner only calls ConfigureFromMaxPrime with primes greater than 1,
@@ -720,8 +567,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         internal ulong[] Exponents { get; private set; } = null!;
 
-        internal ulong[] FilteredDivisors { get; private set; } = null!;
-
         internal GpuDivisorPartialData[] DivisorData { get; private set; } = null!;
 
         internal int[] Offsets { get; private set; } = null!;
@@ -796,7 +641,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
             Divisors = ulongPool.Rent(capacity);
             Exponents = ulongPool.Rent(capacity);
-            FilteredDivisors = ulongPool.Rent(capacity);
             DivisorData = partialPool.Rent(capacity);
             Offsets = intPool.Rent(capacity);
             Counts = intPool.Rent(capacity);
@@ -835,14 +679,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             ArrayPool<GpuDivisorPartialData> partialPool = ThreadStaticPools.GpuDivisorPartialDataPool;
             ulongPool.Return(Divisors, clearArray: false);
             ulongPool.Return(Exponents, clearArray: false);
-            ulongPool.Return(FilteredDivisors, clearArray: false);
             partialPool.Return(DivisorData, clearArray: false);
             intPool.Return(Offsets, clearArray: false);
             intPool.Return(Counts, clearArray: false);
             ulongPool.Return(Cycles, clearArray: false);
             Divisors = Array.Empty<ulong>();
             Exponents = Array.Empty<ulong>();
-            FilteredDivisors = Array.Empty<ulong>();
             DivisorData = Array.Empty<GpuDivisorPartialData>();
             Offsets = Array.Empty<int>();
             Counts = Array.Empty<int>();
diff --git a/PerfectNumbers.Core/HeuristicPrimeSieves.cs b/PerfectNumbers.Core/HeuristicPrimeSieves.cs
new file mode 100644
index 0000000..0f26646
--- /dev/null
+++ b/PerfectNumbers.Core/HeuristicPrimeSieves.cs
@@ -0,0 +1,209 @@
+using System;
+
+namespace PerfectNumbers.Core;
+
+public static class HeuristicPrimeSieves
+{
+    private const int SieveLength = PerfectNumberConstants.MaxQForDivisorCycles;
+    private const int Wheel210 = 210;
+
+    private static readonly int[] GroupAConstants = [3, 7, 11, 13];
+    private static readonly int[] GroupAIncrementPattern = [20, 10];
+
+    private static readonly ushort[] GroupBResiduesEnding1 = [1, 11, 31, 41, 61, 71, 101, 121, 131, 151, 181, 191];
+    private static readonly ushort[] GroupBResiduesEnding3 = [13, 23, 43, 53, 73, 83, 103, 113, 143, 163, 173, 193];
+    private static readonly ushort[] GroupBResiduesEnding7 = [17, 37, 47, 67, 97, 107, 127, 137, 157, 167, 187, 197];
+    private static readonly ushort[] GroupBResiduesEnding9 = [19, 29, 59, 79, 89, 109, 139, 149, 169, 179, 199, 209];
+
+    private static int[]? s_groupADivisors;
+    private static ulong[]? s_groupADivisorSquares;
+
+    private static int[]? s_groupBDivisorsEnding1;
+    private static ulong[]? s_groupBDivisorSquaresEnding1;
+
+    private static int[]? s_groupBDivisorsEnding3;
+    private static ulong[]? s_groupBDivisorSquaresEnding3;
+
+    private static int[]? s_groupBDivisorsEnding7;
+    private static ulong[]? s_groupBDivisorSquaresEnding7;
+
+    private static int[]? s_groupBDivisorsEnding9;
+    private static ulong[]? s_groupBDivisorSquaresEnding9;
+
+    private static readonly object InitLock = new();
+
+    public static ReadOnlySpan<int> GroupADivisors
+    {
+        get
+        {
+            EnsureGroupAInitialized();
+            return s_groupADivisors!;
+        }
+    }
+
+    public static ReadOnlySpan<ulong> GroupADivisorSquares
+    {
+        get
+        {
+            EnsureGroupAInitialized();
+            return s_groupADivisorSquares!;
+        }
+    }
+
+    public static ReadOnlySpan<int> GroupBDivisorsEnding1
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding1, ref s_groupBDivisorSquaresEnding1, GroupBResiduesEnding1);
+            return s_groupBDivisorsEnding1!;
+        }
+    }
+
+    public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding1
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding1, ref s_groupBDivisorSquaresEnding1, GroupBResiduesEnding1);
+            return s_groupBDivisorSquaresEnding1!;
+        }
+    }
+
+    public static ReadOnlySpan<int> GroupBDivisorsEnding3
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding3, ref s_groupBDivisorSquaresEnding3, GroupBResiduesEnding3);
+            return s_groupBDivisorsEnding3!;
+        }
+    }
+
+    public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding3
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding3, ref s_groupBDivisorSquaresEnding3, GroupBResiduesEnding3);
+            return s_groupBDivisorSquaresEnding3!;
+        }
+    }
+
+    public static ReadOnlySpan<int> GroupBDivisorsEnding7
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding7, ref s_groupBDivisorSquaresEnding7, GroupBResiduesEnding7);
+            return s_groupBDivisorsEnding7!;
+        }
+    }
+
+    public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding7
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding7, ref s_groupBDivisorSquaresEnding7, GroupBResiduesEnding7);
+            return s_groupBDivisorSquaresEnding7!;
+        }
+    }
+
+    public static ReadOnlySpan<int> GroupBDivisorsEnding9
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding9, ref s_groupBDivisorSquaresEnding9, GroupBResiduesEnding9);
+            return s_groupBDivisorsEnding9!;
+        }
+    }
+
+    public static ReadOnlySpan<ulong> GroupBDivisorSquaresEnding9
+    {
+        get
+        {
+            EnsureGroupBInitialized(ref s_groupBDivisorsEnding9, ref s_groupBDivisorSquaresEnding9, GroupBResiduesEnding9);
+            return s_groupBDivisorSquaresEnding9!;
+        }
+    }
+
+    private static void EnsureGroupAInitialized()
+    {
+        if (s_groupADivisors is not null)
+        {
+            return;
+        }
+
+        lock (InitLock)
+        {
+            if (s_groupADivisors is not null)
+            {
+                return;
+            }
+
+            var divisors = new int[SieveLength];
+            var squares = new ulong[SieveLength];
+
+            int index = 0;
+            for (; index < GroupAConstants.Length && index < divisors.Length; index++)
+            {
+                int value = GroupAConstants[index];
+                divisors[index] = value;
+                squares[index] = (ulong)value * (ulong)value;
+            }
+
+            int candidate = 23;
+            int incrementIndex = 0;
+            while (index < divisors.Length)
+            {
+                divisors[index] = candidate;
+                squares[index] = (ulong)candidate * (ulong)candidate;
+
+                candidate += GroupAIncrementPattern[incrementIndex];
+                incrementIndex ^= 1;
+                index++;
+            }
+
+            s_groupADivisors = divisors;
+            s_groupADivisorSquares = squares;
+        }
+    }
+
+    private static void EnsureGroupBInitialized(ref int[]? divisorsField, ref ulong[]? squaresField, ReadOnlySpan<ushort> residues)
+    {
+        if (divisorsField is not null)
+        {
+            return;
+        }
+
+        lock (InitLock)
+        {
+            if (divisorsField is not null)
+            {
+                return;
+            }
+
+            var divisors = new int[SieveLength];
+            var squares = new ulong[SieveLength];
+
+            int index = 0;
+            ulong baseValue = 0UL;
+
+            while (index < divisors.Length)
+            {
+                for (int i = 0; i < residues.Length && index < divisors.Length; i++)
+                {
+                    ulong candidate = baseValue + residues[i];
+                    if (candidate <= 13UL)
+                    {
+                        continue;
+                    }
+
+                    divisors[index] = (int)candidate;
+                    squares[index] = candidate * candidate;
+                    index++;
+                }
+
+                baseValue += (ulong)Wheel210;
+            }
+
+            divisorsField = divisors;
+            squaresField = squares;
+        }
+    }
+}
diff --git a/PerfectNumbers.Core/MersenneDivisorResidueStepper.cs b/PerfectNumbers.Core/MersenneDivisorResidueStepper.cs
new file mode 100644
index 0000000..e4a1bf6
--- /dev/null
+++ b/PerfectNumbers.Core/MersenneDivisorResidueStepper.cs
@@ -0,0 +1,104 @@
+using System.Runtime.CompilerServices;
+using PerfectNumbers.Core.Gpu;
+
+namespace PerfectNumbers.Core;
+
+internal struct MersenneDivisorResidueStepper
+{
+    private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
+    private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
+
+    private readonly ushort decimalMask;
+    private readonly byte step10;
+    private readonly byte step8;
+    private readonly byte step5;
+    private readonly byte step3;
+    private readonly byte step7;
+    private readonly byte step11;
+
+    public byte Remainder10 { get; private set; }
+
+    public byte Remainder8 { get; private set; }
+
+    public byte Remainder5 { get; private set; }
+
+    public byte Remainder3 { get; private set; }
+
+    public byte Remainder7 { get; private set; }
+
+    public byte Remainder11 { get; private set; }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public MersenneDivisorResidueStepper(ulong prime, in GpuUInt128 step, in GpuUInt128 firstDivisor)
+    {
+        decimalMask = (prime & 3UL) == 3UL ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
+
+        step10 = ComputeModulo(step, 10, multiplier: 6);
+        step8 = (byte)(step.Low % 8UL);
+        step5 = ComputeModulo(step, 5, multiplier: 1);
+        step3 = ComputeModulo(step, 3, multiplier: 1);
+        step7 = ComputeModulo(step, 7, multiplier: 2);
+        step11 = ComputeModulo(step, 11, multiplier: 5);
+
+        Remainder10 = ComputeModulo(firstDivisor, 10, multiplier: 6);
+        Remainder8 = (byte)(firstDivisor.Low % 8UL);
+        Remainder5 = ComputeModulo(firstDivisor, 5, multiplier: 1);
+        Remainder3 = ComputeModulo(firstDivisor, 3, multiplier: 1);
+        Remainder7 = ComputeModulo(firstDivisor, 7, multiplier: 2);
+        Remainder11 = ComputeModulo(firstDivisor, 11, multiplier: 5);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public readonly bool IsAdmissible()
+    {
+        if (((decimalMask >> Remainder10) & 1) == 0)
+        {
+            return false;
+        }
+
+        if (Remainder8 != 1 && Remainder8 != 7)
+        {
+            return false;
+        }
+
+        if (Remainder3 == 0 || Remainder5 == 0 || Remainder7 == 0 || Remainder11 == 0)
+        {
+            return false;
+        }
+
+        return true;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public void Advance()
+    {
+        Remainder10 = AddMod(Remainder10, step10, 10);
+        Remainder8 = AddMod(Remainder8, step8, 8);
+        Remainder5 = AddMod(Remainder5, step5, 5);
+        Remainder3 = AddMod(Remainder3, step3, 3);
+        Remainder7 = AddMod(Remainder7, step7, 7);
+        Remainder11 = AddMod(Remainder11, step11, 11);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static byte ComputeModulo(in GpuUInt128 value, byte modulus, byte multiplier)
+    {
+        ulong high = value.High % modulus;
+        ulong low = value.Low % modulus;
+        int result = (int)((high * multiplier) + low);
+        result %= modulus;
+        return (byte)result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static byte AddMod(byte value, byte delta, byte modulus)
+    {
+        int sum = value + delta;
+        if (sum >= modulus)
+        {
+            sum -= modulus;
+        }
+
+        return (byte)sum;
+    }
+}
diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
index ac93691..3e3d5f6 100644
--- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
+++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
@@ -1055,6 +1055,8 @@ internal static partial class PrimeOrderCalculator
 		List<PendingEntry> pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2);
 		Stack<ulong>? compositeStack = null;
 		PartialFactorResult result;
+		PrimeTester primeTester = PrimeTester.Exclusive;
+
 		if (!gpuFactored)
 		{
 			// We don't need to worry about leftovers, because we always use indexes within the calculated counts
@@ -1120,7 +1122,7 @@ internal static partial class PrimeOrderCalculator
 					// 	continue;
 					// }
 
-					bool isPrime = PrimeTester.IsPrimeGpu(composite);
+					bool isPrime = primeTester.HeuristicIsPrimeGpu(composite);
 					// bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
 					// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
 
@@ -1182,9 +1184,10 @@ internal static partial class PrimeOrderCalculator
 
 			if (!entry.HasKnownPrimality)
 			{
-				bool isPrime = PrimeTester.IsPrimeGpu(composite);
+				bool isPrime = primeTester.HeuristicIsPrimeGpu(composite);
 				// bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
 				// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
+				
 				entry = entry.WithPrimality(isPrime);
 				pending[index] = entry;
 			}
@@ -1213,7 +1216,7 @@ internal static partial class PrimeOrderCalculator
 		}
 		else
 		{
-			cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor);
+			cofactorIsPrime = primeTester.HeuristicIsPrimeGpu(cofactor);
 			// cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None);
 			// cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor);
 		}
@@ -1743,7 +1746,9 @@ internal static partial class PrimeOrderCalculator
 			return;
 		}
 
-		bool isPrime = PrimeTester.IsPrimeGpu(value);
+		PrimeTester primeTester = PrimeTester.Exclusive;
+
+		bool isPrime = primeTester.HeuristicIsPrimeGpu(value);
 		// bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None);
 		// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value);
 		if (!knownComposite && isPrime)
@@ -1755,7 +1760,7 @@ internal static partial class PrimeOrderCalculator
 		ulong factor = PollardRhoStrict(value);
 		ulong quotient = value / factor;
 
-		isPrime = PrimeTester.IsPrimeGpu(factor);
+		isPrime = primeTester.HeuristicIsPrimeGpu(factor);
 		// isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None);
 		// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor);
 		if (isPrime)
@@ -1786,7 +1791,7 @@ internal static partial class PrimeOrderCalculator
 		}
 		else
 		{
-			isPrime = PrimeTester.IsPrimeGpu(quotient);
+			isPrime = primeTester.HeuristicIsPrimeGpu(quotient);
 			// isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None);
 			// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient);
 		}
diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs
index 15af30e..e82d0eb 100644
--- a/PerfectNumbers.Core/PrimeTester.cs
+++ b/PerfectNumbers.Core/PrimeTester.cs
@@ -1,5 +1,8 @@
+using System;
 using System.Buffers;
 using System.Collections.Concurrent;
+using Open.Numeric.Primes;
+using System.Collections.Generic;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using ILGPU;
@@ -11,15 +14,96 @@ namespace PerfectNumbers.Core;
 
 public sealed class PrimeTester(bool useInternal = false)
 {
-    public bool IsPrime(ulong n, CancellationToken ct) =>
-        // TODO: Eliminate this wrapper once callers can reach IsPrimeInternal directly; the extra
-        // indirection shows up when we sieve millions of candidates per second.
-        IsPrimeInternal(n, ct);
+    // 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];
+    private static readonly byte[] GroupBEndingOrderMod9 = [9, 7, 1];
 
-    public static bool IsPrimeGpu(ulong n) => new PrimeTester().IsPrimeGpu(n, CancellationToken.None);
+    // 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;
+    internal const int MaxGroupBSequences = 36;
+
+    private const bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B.
+
+    private readonly bool _useLegacyPrimeTester = useInternal;
+
+    private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
+
+    internal enum HeuristicDivisorGroup : byte
+    {
+        None = 0,
+        GroupAConstant = 1,
+        GroupAWheel = 2,
+        GroupB = 3,
+    }
+
+    internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue)
+    {
+        public readonly ulong Value = value;
+        public readonly HeuristicDivisorGroup Group = group;
+        public readonly byte Ending = ending;
+        public readonly byte PriorityIndex = priorityIndex;
+        public readonly ushort WheelResidue = wheelResidue;
+    }
+
+    internal readonly struct HeuristicDivisorPreparation(
+        in HeuristicDivisorCandidate candidate,
+        in MontgomeryDivisorData divisorData,
+        ulong cycleLengthHint,
+        bool hasCycleLengthHint)
+    {
+        public readonly HeuristicDivisorCandidate Candidate = candidate;
+        public readonly MontgomeryDivisorData DivisorData = divisorData;
+        public readonly ulong CycleLengthHint = cycleLengthHint;
+        public readonly bool HasCycleLengthHint = hasCycleLengthHint;
+        public readonly bool RequiresCycleComputation = !hasCycleLengthHint;
+    }
+
+    public bool IsPrime(ulong n, CancellationToken ct)
+    {
+        if (_useLegacyPrimeTester)
+        {
+            return LegacyIsPrimeInternal(n, ct);
+        }
+
+        return IsPrimeInternal(n, ct);
+    }
+
+    [ThreadStatic]
+    private static PrimeTester? _tester;
+
+    public static PrimeTester Exclusive => _tester ??= new();
+
+    public static bool IsPrimeGpu(ulong n)
+    {
+        return Exclusive.IsPrimeGpu(n, CancellationToken.None);
+    }
+
+    public static bool IsPrimeGpu(ulong n, ulong limit, byte nMod10)
+    {
+        return Exclusive.HeuristicIsPrimeGpu(n, limit, nMod10);
+    }
 
     // Optional GPU-assisted primality: batched small-prime sieve on device.
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public bool IsPrimeGpu(ulong n, CancellationToken ct)
+    {
+        if (_useLegacyPrimeTester || !EnableHeuristicPrimeTesting)
+        {
+            return LegacyIsPrimeGpu(n, ct);
+        }
+
+        return HeuristicIsPrimeGpu(n);
+    }
+
+    public static bool LegacyIsPrimeGpu(ulong n, CancellationToken ct)
     {
         // EvenPerfectBitScanner never enqueues exponents <= 3 on the GPU path; keep the legacy guard commented out
         // so ad-hoc callers that expect that behavior must re-enable it intentionally.
@@ -34,11 +118,11 @@ public sealed class PrimeTester(bool useInternal = false)
 
         // EvenPerfectBitScanner seeds p at 136,279,841 and advances monotonically, so production scans never hit the GPU with
         // exponents below 31. Leave this CPU redirect commented out to keep the wrapper branchless while documenting the guard.
-        // The ternary fallback below still routes sub-31 values through IsPrimeInternal for correctness when ad-hoc callers bypass
+        // The ternary fallback below still routes sub-31 values through LegacyIsPrimeInternal for correctness when ad-hoc callers bypass
         // the scanner pipeline.
         // if (n < 31UL)
         // {
-        //     return IsPrimeInternal(n, ct);
+        //     return LegacyIsPrimeInternal(n, ct);
         // }
 
         bool forceCpu = GpuContextPool.ForceCpu;
@@ -58,16 +142,611 @@ public sealed class PrimeTester(bool useInternal = false)
         bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0;
         bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime;
 
-        // Defensive fallback: GPU sieve produced a composite verdict or execution was skipped.
-        // Confirm with CPU logic to prevent false negatives observed on some accelerators.
-        return requiresCpuFallback ? IsPrimeInternal(n, ct) : true;
+        return requiresCpuFallback ? LegacyIsPrimeInternal(n, ct) : true;
+    }
+
+    public bool HeuristicIsPrimeGpu(ulong n)
+    {
+        byte nMod10 = (byte)n.Mod10();
+        ulong sqrtLimit = ComputeHeuristicSqrt(n);
+        return HeuristicIsPrimeGpu(n, sqrtLimit, nMod10);
+    }
+
+    private bool HeuristicIsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10)
+    {
+        // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner?
+        if (sqrtLimit < 3UL)
+        {
+            return EvaluateWithOpenNumericFallback(n);
+        }
+
+        if (!UseHeuristicGroupBTrialDivision)
+        {
+            return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false);
+        }
+
+        if (GpuContextPool.ForceCpu)
+        {
+            return HeuristicIsPrimeCpu(n, sqrtLimit, nMod10);
+        }
+
+        bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10);
+        return !compositeDetected;
+    }
+
+    internal static bool HeuristicIsPrimeCpu(ulong n)
+    {
+        return HeuristicIsPrimeCpu(n, 0UL, (byte)n.Mod10());
+    }
+
+    internal static bool HeuristicIsPrimeCpu(ulong n, ulong sqrtLimit, byte nMod10)
+    {
+        if (sqrtLimit == 0UL)
+        {
+            sqrtLimit = ComputeHeuristicSqrt(n);
+        }
+
+        if (sqrtLimit < 3UL)
+        {
+            return EvaluateWithOpenNumericFallback(n);
+        }
+
+        bool includeGroupB = UseHeuristicGroupBTrialDivision;
+        return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB);
+    }
+
+    private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB)
+    {
+        Span<HeuristicGroupBSequenceState> groupBBuffer = stackalloc HeuristicGroupBSequenceState[MaxGroupBSequences];
+        var enumerator = new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
+
+        while (enumerator.TryGetNext(out HeuristicDivisorCandidate candidate))
+        {
+            if (!includeGroupB && candidate.Group == HeuristicDivisorGroup.GroupB)
+            {
+                return EvaluateWithOpenNumericFallback(n);
+            }
+
+            ulong divisor = candidate.Value;
+            if (divisor <= 1UL)
+            {
+                continue;
+            }
+
+            if (n % divisor == 0UL)
+            {
+                return false;
+            }
+        }
+
+        return includeGroupB ? true : EvaluateWithOpenNumericFallback(n);
+    }
+
+    private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10)
+    {
+        Span<HeuristicGroupBSequenceState> groupBBuffer = stackalloc HeuristicGroupBSequenceState[MaxGroupBSequences];
+        var enumerator = new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
+
+        int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize);
+        var candidatePool = ArrayPool<HeuristicDivisorCandidate>.Shared;
+        var divisorPool = ThreadStaticPools.UlongPool;
+        var hitPool = ThreadStaticPools.BytePool;
+
+        HeuristicDivisorCandidate[]? candidateArray = null;
+        ulong[]? divisorArray = null;
+        byte[]? hitFlags = null;
+
+        var limiter = GpuPrimeWorkLimiter.Acquire();
+        var gpu = PrimeTesterGpuContextPool.Rent();
+
+        bool compositeDetected = false;
+
+        try
+        {
+            candidateArray = candidatePool.Rent(batchCapacity);
+            divisorArray = divisorPool.Rent(batchCapacity);
+            hitFlags = hitPool.Rent(batchCapacity);
+
+            var accelerator = gpu.Accelerator;
+            var state = GpuKernelState.GetOrCreate(accelerator);
+
+            lock (gpu.ExecutionLock)
+            {
+                var scratch = state.RentScratch(batchCapacity, accelerator);
+                try
+                {
+                    int count = 0;
+
+                    bool ProcessBatch(int length)
+                    {
+                        scratch.Input.View.CopyFromCPU(ref divisorArray![0], length);
+                        state.HeuristicTrialDivisionKernel(length, scratch.Input.View, n, scratch.Output.View);
+                        accelerator.Synchronize();
+                        scratch.Output.View.CopyToCPU(ref hitFlags![0], length);
+
+                        for (int i = 0; i < length; i++)
+                        {
+                            if (hitFlags![i] == 0)
+                            {
+                                continue;
+                            }
+
+                            ulong divisor = candidateArray![i].Value;
+                            if (divisor > 1UL && n % divisor == 0UL)
+                            {
+                                return true;
+                            }
+                        }
+
+                        return false;
+                    }
+
+                    while (enumerator.TryGetNext(out HeuristicDivisorCandidate candidate))
+                    {
+                        ulong divisor = candidate.Value;
+                        if (divisor <= 1UL)
+                        {
+                            continue;
+                        }
+
+                        candidateArray[count] = candidate;
+                        divisorArray[count] = divisor;
+                        count++;
+
+                        if (count == batchCapacity)
+                        {
+                            if (ProcessBatch(count))
+                            {
+                                compositeDetected = true;
+                                break;
+                            }
+
+                            count = 0;
+                        }
+                    }
+
+                    if (!compositeDetected && count > 0)
+                    {
+                        if (ProcessBatch(count))
+                        {
+                            compositeDetected = true;
+                        }
+                    }
+                }
+                finally
+                {
+                    state.ReturnScratch(scratch);
+                }
+            }
+        }
+        finally
+        {
+            if (hitFlags is not null)
+            {
+                hitPool.Return(hitFlags);
+            }
+
+            if (divisorArray is not null)
+            {
+                divisorPool.Return(divisorArray);
+            }
+
+            if (candidateArray is not null)
+            {
+                candidatePool.Return(candidateArray);
+            }
+
+            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;
+    }
+
+    private static ReadOnlySpan<byte> GetGroupBEndingOrder(byte nMod10) => nMod10 switch
+    {
+        1 => GroupBEndingOrderMod1,
+        3 => GroupBEndingOrderMod3,
+        7 => GroupBEndingOrderMod7,
+        9 => GroupBEndingOrderMod9,
+        _ => ReadOnlySpan<byte>.Empty,
+    };
+
+    private static ReadOnlySpan<ushort> GetWheel210ResiduesForEnding(byte ending) => ending switch
+    {
+        1 => Wheel210ResiduesEnding1,
+        7 => Wheel210ResiduesEnding7,
+        9 => Wheel210ResiduesEnding9,
+        _ => ReadOnlySpan<ushort>.Empty,
+    };
+
+    private static int InitializeGroupBStates(byte nMod10, Span<HeuristicGroupBSequenceState> buffer)
+    {
+        ReadOnlySpan<byte> endings = GetGroupBEndingOrder(nMod10);
+        int count = 0;
+
+        for (int priorityIndex = 0; priorityIndex < endings.Length; priorityIndex++)
+        {
+            ReadOnlySpan<ushort> residues = GetWheel210ResiduesForEnding(endings[priorityIndex]);
+            for (int residueIndex = 0; residueIndex < residues.Length; residueIndex++)
+            {
+                buffer[count++] = new HeuristicGroupBSequenceState(residues[residueIndex], endings[priorityIndex], (byte)priorityIndex);
+            }
+        }
+
+        return count;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
+    {
+        return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor)
+    {
+        return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate)
+    {
+        ulong divisor = candidate.Value;
+        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor);
+        bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength);
+        return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static ulong ResolveHeuristicCycleLength(
+        ulong exponent,
+        in HeuristicDivisorPreparation preparation,
+        out bool cycleFromHint,
+        out bool cycleComputed,
+        out bool primeOrderFailed)
+    {
+        if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL)
+        {
+            cycleFromHint = true;
+            cycleComputed = true;
+            primeOrderFailed = false;
+            return preparation.CycleLengthHint;
+        }
+
+        ulong divisor = preparation.Candidate.Value;
+        MontgomeryDivisorData divisorData = preparation.DivisorData;
+
+        bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
+            divisor,
+            exponent,
+            divisorData,
+            out ulong computedCycle,
+            out bool primeOrderFailedLocal);
+
+        if (trySuccess && computedCycle != 0UL)
+        {
+            cycleFromHint = false;
+            cycleComputed = true;
+            primeOrderFailed = primeOrderFailedLocal;
+            return computedCycle;
+        }
+
+        primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL;
+        ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed);
+        cycleFromHint = false;
+        cycleComputed = resolvedCycle != 0UL;
+        return resolvedCycle;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength)
+    {
+        ulong[] snapshot = HeuristicSmallCycleSnapshot;
+        if (divisor < (ulong)snapshot.Length)
+        {
+            cycleLength = snapshot[(int)divisor];
+            if (cycleLength != 0UL)
+            {
+                return true;
+            }
+        }
+
+        cycleLength = 0UL;
+        return false;
+    }
+
+    internal ref struct HeuristicDivisorEnumerator
+    {
+        private readonly ulong sqrtLimit;
+        private int groupAConstantIndex;
+        private ulong groupAWheelCandidate;
+        private int groupAIncrementIndex;
+        private bool groupAWheelActive;
+        private Span<HeuristicGroupBSequenceState> groupBStates;
+
+        public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span<HeuristicGroupBSequenceState> groupBBuffer)
+        {
+            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++)
+            {
+                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;
+            }
+
+            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;
+
+                    ulong next = value + increment;
+                    if (next <= value)
+                    {
+                        groupAWheelActive = false;
+                    }
+                    else
+                    {
+                        groupAWheelCandidate = next;
+                    }
+
+                    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;
+            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;
+
+		public void TryNormalize()
+        {
+            ulong candidate = Candidate;
+            if (candidate == ulong.MaxValue)
+            {
+                return;
+            }
+
+            while (candidate <= 13UL)
+            {
+                ulong next = candidate + Wheel210;
+                if (next <= candidate)
+                {
+                    Candidate = ulong.MaxValue;
+                    return;
+                }
+
+                candidate = next;
+            }
+
+            Candidate = candidate;
+        }
+
+        public void Advance()
+        {
+            ulong candidate = Candidate;
+            if (candidate == ulong.MaxValue)
+            {
+                return;
+            }
+
+            ulong next = candidate + Wheel210;
+            Candidate = next > candidate ? next : ulong.MaxValue;
+        }
+
+
+    }
+
+    internal ref struct MersenneHeuristicDivisorEnumerator
+    {
+        private readonly GpuUInt128 step;
+        private readonly GpuUInt128 limit;
+        private GpuUInt128 current;
+        private MersenneDivisorResidueStepper residueStepper;
+        private bool active;
+        private ulong processedCount;
+        private ulong lastDivisor;
+
+        public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor)
+        {
+            var stepLocal = new GpuUInt128(exponent);
+            stepLocal.ShiftLeft(1);
+            step = stepLocal;
+
+            var limitLocal = new GpuUInt128(maxDivisor);
+            limit = limitLocal;
+
+            var firstDivisor = stepLocal;
+            firstDivisor.Add(1UL);
+
+            bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0;
+
+            current = hasCandidates ? firstDivisor : GpuUInt128.Zero;
+            residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default;
+            active = hasCandidates;
+            processedCount = 0UL;
+            lastDivisor = 0UL;
+        }
+
+        public bool TryGetNext(out HeuristicDivisorCandidate candidate)
+        {
+            while (active)
+            {
+                ulong value = current.Low;
+                processedCount++;
+                lastDivisor = value;
+
+                bool admissible = residueStepper.IsAdmissible();
+
+                Advance();
+
+                if (admissible)
+                {
+                    candidate = CreateCandidate(value);
+                    return true;
+                }
+            }
+
+            candidate = default;
+            return false;
+        }
+
+        private void Advance()
+        {
+            GpuUInt128 next = current;
+            next.Add(step);
+
+            if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0)
+            {
+                active = false;
+                return;
+            }
+
+            current = next;
+            residueStepper.Advance();
+        }
+
+        public ulong ProcessedCount => processedCount;
+
+        public ulong LastDivisor => lastDivisor;
+
+        public bool Exhausted => !active;
+
+        private static HeuristicDivisorCandidate CreateCandidate(ulong value)
+        {
+            byte ending = (byte)(value % 10UL);
+            HeuristicDivisorGroup group;
+
+            if (value == 3UL || value == 7UL || value == 11UL || value == 13UL)
+            {
+                group = HeuristicDivisorGroup.GroupAConstant;
+            }
+            else if (ending == 3)
+            {
+                group = HeuristicDivisorGroup.GroupAWheel;
+            }
+            else
+            {
+                group = HeuristicDivisorGroup.GroupB;
+            }
+
+            ushort residue = (ushort)(value % Wheel210);
+            return new HeuristicDivisorCandidate(value, group, ending, 0, residue);
+        }
     }
 
 
-    // Note: GPU-backed primality path is implemented via IsPrimeGpu/IsPrimeBatchGpu and is routed
-    // from EvenPerfectBitScanner based on --primes-device.
-
     internal static bool IsPrimeInternal(ulong n, CancellationToken ct)
+    {
+        if (!EnableHeuristicPrimeTesting)
+        {
+            return LegacyIsPrimeInternal(n, ct);
+        }
+
+        return HeuristicIsPrimeCpu(n);
+    }
+
+    public static bool LegacyIsPrimeInternal(ulong n, CancellationToken ct)
     {
         // EvenPerfectBitScanner streams monotonically increasing odd exponents that start at 136,279,841 and already exclude
         // multiples of five. Factorization helpers (PrimeOrderCalculator, Pollard routines, and divisor-cycle warmups) reuse
@@ -113,8 +792,13 @@ public sealed class PrimeTester(bool useInternal = false)
         return result;
     }
 
+    public static bool EnableHeuristicPrimeTesting { get; set; } = true;
+
     public static int GpuBatchSize { get; set; } = 262_144;
 
+    public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096;
+
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static void IsPrimeBatchGpu(ReadOnlySpan<ulong> values, Span<byte> results)
     {
@@ -150,7 +834,7 @@ public sealed class PrimeTester(bool useInternal = false)
                         int count = remaining > batchSize ? batchSize : remaining;
 
                         values.Slice(pos, count).CopyTo(temp);
-                        input.View.CopyFromCPU(ref temp[0], count);
+						input.View.CopyFromCPU(ref temp[0], count);
 
                         state.Kernel(count, input.View, state.DevicePrimes.View, output.View);
                         accelerator.Synchronize();
@@ -263,6 +947,7 @@ public sealed class PrimeTester(bool useInternal = false)
     private sealed class KernelState
     {
         public Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<byte>> Kernel { get; }
+        public Action<Index1D, ArrayView<ulong>, ulong, ArrayView<byte>> HeuristicTrialDivisionKernel { get; }
         public MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimes { get; }
         private readonly Accelerator _accel;
         private readonly System.Collections.Concurrent.ConcurrentBag<ScratchBuffers> _scratchPool = [];
@@ -275,6 +960,7 @@ public sealed class PrimeTester(bool useInternal = false)
             _accel = accelerator;
             // Compile once per accelerator and upload primes once.
             Kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<byte>>(PrimeTesterKernels.SmallPrimeSieveKernel);
+            HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ulong, ArrayView<byte>>(PrimeTesterKernels.HeuristicTrialDivisionKernel);
 
             var primes = PrimesGenerator.SmallPrimes;
             DevicePrimes = accelerator.Allocate1D<uint>(primes.Length);
diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md
new file mode 100644
index 0000000..fd2f1ca
--- /dev/null
+++ b/docs/HeuristicIsPrimePlan.md
@@ -0,0 +1,158 @@
+# EvenPerfectBitScanner `--mersenne=bydivisor` Flow and Heuristic Prime Tester Plan
+
+## 1. Execution Flow in `--mersenne=bydivisor` Mode
+
+The following map captures every branch that executes when the scanner runs with `--mersenne=bydivisor`, highlighting interactions with prime testing, divisor heuristics, residue reuse, and factor verification.
+
+1. **CLI processing and guards**
+   * `CliArguments.Parse` fills `_cliArguments`, then `Main` validates thread count, help flag, and shared GPU settings.
+   * `useByDivisor` is set from CLI. If enabled without `--filter-p` (and not in test mode), the scanner prints "--mersenne=bydivisor requires --filter-p=<path>." and exits early.【F:EvenPerfectBitScanner/Program.cs†L32-L132】
+
+2. **Global configuration**
+   * GPU settings (`GpuPrimeWorkLimiter`, `PrimeTester.GpuBatchSize`, transforms) are configured for both CPU and GPU runs.
+   * By-divisor specific state records `_byDivisorStartPrime`, `_byDivisorTester`, and `_byDivisorPreviousResults` for later reuse.【F:EvenPerfectBitScanner/Program.cs†L134-L232】
+
+3. **Tester initialization**
+   * When `useByDivisor` is active, `_byDivisorTester` is instantiated as GPU or CPU implementation, with `BatchSize` set from `--scan-batch`. `PrimeTester`/`MersenneNumberTester` pools are not created for this mode because the flow exits before the generic candidate loop.【F:EvenPerfectBitScanner/Program.cs†L200-L232】【F:EvenPerfectBitScanner/Program.cs†L400-L460】
+
+4. **Results bootstrap**
+   * The scanner builds the results filename, initializes `CalculationResultHandler`, and loads previous results. Previous entries populate `_byDivisorPreviousResults` so the by-divisor runner can skip already processed primes.【F:EvenPerfectBitScanner/Program.cs†L232-L318】
+
+5. **Candidate acquisition**
+   * `byDivisorCandidates` is populated either from deterministic test data (test mode) or by reloading the filter file with `CalculationResultsFile.LoadCandidatesWithinRange`. Limits (`--max-prime`) and skipped candidates are reported to the console.【F:EvenPerfectBitScanner/Program.cs†L318-L362】
+   * In contrast to other modes, generic residue filters and candidate queues are bypassed; only the provided prime list is considered.【F:EvenPerfectBitScanner/Program.cs†L318-L362】
+
+6. **Previous-result pruning**
+   * `MersenneNumberDivisorByDivisorTester.Run` sorts the input primes, removes entries already present in `_byDivisorPreviousResults`, applies `--start-prime`, and converts the remaining list to spans for batch preparation.【F:PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs†L11-L138】
+
+7. **Tester configuration per maximum prime**
+   * `IMersenneNumberDivisorByDivisorTester.ConfigureFromMaxPrime` and `PrepareCandidates` pre-compute divisor limits for every prime. All primes with `allowedMax >= 3` are retained (current implementations never filter any out in production because limits are huge).【F:PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs†L140-L201】
+
+8. **Residue-sieved divisor scanning (CPU)**
+   * `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` initialises Montgomery stepping data (`GpuUInt128`), computes rolling residues for mod 10/8/5/3/7/11, and advances divisors by `2p` while applying decimal masks tied to the exponent’s last digit.【F:PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs†L63-L213】
+   * When a candidate passes residue filters, the CPU session obtains cycle lengths from `DivisorCycleCache` or `MersenneDivisorCycles`, validates them via Montgomery exponentiation, and reports hits while updating `_lastStatusDivisor`.【F:PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs†L213-L318】
+
+9. **Residue-sieved divisor scanning (GPU)**
+   * `MersenneNumberDivisorByDivisorGpuTester.CheckDivisors` reuses the same wheel offsets and residue tests but performs them inside `DivisorByDivisorKernels.CheckKernel`, writing hits to pooled GPU buffers. CPU fallbacks honour the same residue, decimal mask, and cycle rules when GPU context is unavailable.【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L1-L240】【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L240-L520】
+
+10. **Prime decision path**
+    * By-divisor sessions internally reuse divisor remainders across primes, but they **do not call `PrimeTester.IsPrime*`**. The primality decision depends entirely on divisor scans up to the limit configured for each prime candidate.
+    * If `IsPrime` returns `false`, the result file records `(detailedCheck=true, passedAllTests=false)`. When `true`, `divisorsExhausted` reports whether the scan reached the configured divisor bound (driving the `DetailedCheck` flag).【F:PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs†L203-L272】
+
+11. **Shutdown**
+    * After `MersenneNumberDivisorByDivisorTester.Run` completes, the main loop flushes results and exits before any of the generic candidate enumeration logic runs.【F:EvenPerfectBitScanner/Program.cs†L400-L460】
+
+## 2. Baseline Prime Tester Methods
+
+The heuristics will replace the internal logic of the current GPU-assisted and CPU-only paths. The snapshots below capture the existing implementations before modification (copied verbatim for reference).
+
+```csharp
+public bool IsPrimeGpu(ulong n, CancellationToken ct)
+{
+    bool forceCpu = GpuContextPool.ForceCpu;
+    Span<ulong> one = stackalloc ulong[1];
+    Span<byte> outFlags = stackalloc byte[1];
+    one[0] = n;
+    outFlags[0] = 0;
+
+    if (!forceCpu)
+    {
+        IsPrimeBatchGpu(one, outFlags);
+    }
+
+    bool belowGpuRange = n < 31UL;
+    bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0;
+    bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime;
+
+    return requiresCpuFallback ? IsPrimeInternal(n, ct) : true;
+}
+```
+
+```csharp
+internal static bool IsPrimeInternal(ulong n, CancellationToken ct)
+{
+    bool isTwo = n == 2UL;
+    bool isOdd = (n & 1UL) != 0UL;
+    ulong mod5 = n % 5UL;
+    bool divisibleByFive = n > 5UL && mod5 == 0UL;
+
+    bool result = n >= 2UL && (isTwo || isOdd) && !divisibleByFive;
+    bool requiresTrialDivision = result && n >= 7UL && !isTwo;
+
+    if (requiresTrialDivision)
+    {
+        bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n);
+        result &= !sharesMaxExponentFactor;
+
+        if (result)
+        {
+            var smallPrimeDivisorsLength = PrimesGenerator.SmallPrimes.Length;
+            uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes;
+            ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2;
+            for (int i = 0; i < smallPrimeDivisorsLength; i++)
+            {
+                if (smallPrimeDivisorsMul[i] > n)
+                {
+                    break;
+                }
+
+                if (n % smallPrimeDivisors[i] == 0)
+                {
+                    result = false;
+                    break;
+                }
+            }
+        }
+    }
+
+    return result;
+}
+```
+
+## 3. Implementation Plan for `PrimeTester.HeuristicIsPrime*`
+
+The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, residue priorities) and connects them to both CPU and GPU execution, including every by-divisor integration point. Items already completed are marked with `[done]`.
+
+1. [done] **Snapshot current implementations** – preserve `IsPrimeGpu` and `IsPrimeInternal` as references for regression comparisons and migration checkpoints (see Section 2).
+2. [done] **Introduce placeholder heuristic entry points** – add `HeuristicIsPrimeCpu` and `HeuristicIsPrimeGpu` by copying the current internal and GPU methods so subsequent commits can iterate without disrupting legacy behavior.
+3. [done] **Design shared heuristic scaffolding**
+   * [done] Introduced `HeuristicDivisorEnumerator` and stackalloc-backed Group B sequence states so CPU/GPU heuristics share the Group A/B wheel metadata and candidate ordering.
+   * [done] Surface helpers that accept the precomputed `R = ⌊√n⌋` and `nMod10` values from callers via the new overload of `HeuristicIsPrimeCpu`; PrimeTester still computes them locally until call sites forward the parameters.
+   * [done] Simplified the heuristic execution path so it runs without cancellation tokens or callback plumbing while keeping the shared enumerator available to other components.
+4. **Share residue and wheel state with by-divisor sessions**
+   * [done] Extracted the rolling-modulus helpers (mod 10/8/5/3/7/11) and decimal masks from `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` into the shared `MersenneDivisorResidueStepper` so heuristic prime tests and by-divisor sessions consume the same filters.
+   * [done] Audited how `MersenneCpuDivisorScanSession` advances residues via `ExponentRemainderStepperCpu`/`CycleRemainderStepper` and updated the GPU by-divisor path to reuse `MersenneDivisorResidueStepper`, keeping CPU/GPU residue stepping aligned.
+   * [done] Publish wheel iterators that can stream candidates directly into `DivisorCycleCache`, GPU hit buffers, or CPU Montgomery checks without recomputing offsets by exposing `PrimeTester.CreateHeuristicDivisorEnumerator` and the reusable `HeuristicGroupBSequenceState` buffer.
+5. **Expose cycle-preparation helpers**
+   * [done] Extend the heuristic scaffolding so it returns `MontgomeryDivisorData` descriptors and cycle-length hints for each candidate divisor via `HeuristicDivisorPreparation` and `PrepareHeuristicDivisor`.
+   * Wire these outputs into the existing cycle calculators to avoid duplicate cache lookups inside the CPU/GPU divisor scans.
+6. [done] **Implement `HeuristicIsPrimeCpu`**
+   * Baseline Group A/B enumeration now executes inside `HeuristicIsPrimeCpu` via `HeuristicDivisorEnumerator`, computing `⌊√n⌋` locally and short-circuiting on the first divisor.
+   * [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 `HeuristicIsPrimeGpu`**
+   * [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, reusing caller-provided `sqrtLimit`/`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.
+8. [done] **Refactor by-divisor CPU/GPU testers to reuse heuristics**
+   * [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `PrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop while preserving the existing status counters.
+   * [done] Streamed GPU batches from the same enumerator and Montgomery preparation data, removing the filtered-divisor scratch arrays while continuing to feed the existing kernel and hit bookkeeping.
+9. [done] **Route legacy prime checks through heuristics**
+   * [done] Update `IsPrimeInternal` and `IsPrimeGpu` to delegate to the heuristic implementations when `EnableHeuristicPrimeTesting` is enabled, while allowing `PrimeTester(useInternal: true)` to force the legacy CPU/GPU paths.
+   * [done] Maintain legacy paths for benchmarking via `LegacyIsPrimeInternal`, `LegacyIsPrimeGpu`, and the global `EnableHeuristicPrimeTesting` switch so benchmark harnesses can compare behaviours.
+10. **Validation and benchmarking**
+    * Create targeted unit tests covering boundary cases for Group A/B transitions, wheel enumerations, residue deltas, and the stop-on-hit behaviour.
+    * Expand GPU test harnesses to assert equivalence between CPU/GPU heuristics on sampled composite and prime inputs and to verify the residue-stream contract with by-divisor sessions.
+    * Benchmark against historical data (hundreds of thousands of composites) to confirm that heuristic ordering improves early divisor detection and that residue sharing reduces duplicate computation.
+    * [done] Build BenchmarkDotNet comparisons for `PrimeTester.IsPrimeGpu`, `PrimeTester.IsPrimeInternal`, and `Open.Numeric.Primes.Prime.Numbers.IsPrime` across representative ranges (≤100, ≤4,000,000, and ≥138,000,000) before routing heuristic code into production. Implemented in `EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs` to exercise both CPU and GPU heuristic paths alongside the external library.
+11. **Documentation and rollout**
+    * Update README/docs to explain heuristic prerequisites (odd, not divisible by 5) and how configuration toggles affect execution across prime testing and by-divisor scans.
+    * Provide migration guidance for operators toggling between legacy and heuristic modes and describe how telemetry/metrics change after integration.
+    * Publish benchmark results in the heuristic specification once measurements are available so the documented guidance reflects the observed CPU/GPU/external comparisons.
+
+## 4. Remaining Work Summary
+
+* Complete the remaining rollout tasks in steps 10–11, including legacy-call-site integration, validation, and documentation updates.
+* Develop comprehensive validation suites (unit tests and performance benchmarks) to compare heuristic behaviour with current implementations across both prime testing and by-divisor divisor scans.
+* Document configuration switches, runtime metrics, and operational playbooks once the heuristics are in place.
diff --git a/docs/HeuristicIsPrimeSpecification.md b/docs/HeuristicIsPrimeSpecification.md
new file mode 100644
index 0000000..f365ad9
--- /dev/null
+++ b/docs/HeuristicIsPrimeSpecification.md
@@ -0,0 +1,107 @@
+# Heuristic Prime Tester Specification
+
+This document describes the heuristic prime testing procedure implemented by `PrimeTester.HeuristicIsPrimeCpu` and `PrimeTester.HeuristicIsPrimeGpu`. The guidance consolidates the residue analysis and empirical statistics for odd inputs that are not divisible by five, based on the first 500,000 odd composite numbers after removing multiples of five. It also records the integration contracts required by the by-divisor pipeline so the heuristic covers every path that participates in divisor searches, residue sieves, and factorisation checks.
+
+## Input Assumptions
+
+* `n` is odd, `n > 2`, and `n` is not divisible by 5.
+* The goal is to **prove compositeness** by finding a single proper divisor.
+* Trial division always stops at `⌊√n⌋`, which the caller computes and supplies alongside the candidate.
+
+## Candidate Partitioning
+
+### Group A (special case)
+
+```
+A = {3, 7, 11, 13} ∪ { d : d ≡ 3 (mod 10), d ≤ ⌊√n⌋ }
+```
+
+Group A always executes first.
+
+### Group B (remaining candidates)
+
+Group B contains every other proper divisor `≤ ⌊√n⌋` that does **not** end with digit 3.
+
+#### Allowed endings for Group B (after removing A)
+
+The allowed endings depend on `n mod 10`:
+
+* `n ≡ 1 (mod 10)` → `{1, 9}`
+* `n ≡ 3 (mod 10)` → `{7, 9}`
+* `n ≡ 7 (mod 10)` → `{1, 7}`
+* `n ≡ 9 (mod 10)` → `{1, 7, 9}`
+
+These sets are exhaustive. Do **not** add other endings.
+
+## Wheel Strategy
+
+* **Group A** uses the wheel based on 30 (residues coprime with 2 and 5: `{1, 7, 11, 13, 17, 19, 23, 29}`) combined with the `d ≡ 3 (mod 10)` filter. The constants `{3, 7, 11, 13}` are checked explicitly before wheel generation. The implementation uses wheel 30; a wheel of 210 may also be used for full uniformity, but 30 is lighter and already filtered to the `≡ 3 (mod 10)` residue.
+* **Group B** uses a wheel based on 210 (coprime with 2, 3, 5, and 7). Only numbers whose final digit belongs to the allowed-ending set for the current `n mod 10` are retained.
+
+## Testing Order
+
+After Group A finishes, Group B enumerates candidates in strictly increasing order up to `⌊√n⌋`. Within that ascending order, endings are prioritised based on empirical hit rates **after excluding Group A**:
+
+* `n ≡ 1 (mod 10)` → try endings `9` first, then `1`.
+* `n ≡ 3 (mod 10)` → endings `9` then `7`.
+* `n ≡ 7 (mod 10)` → endings `7` then `1`.
+* `n ≡ 9 (mod 10)` → endings `9`, then `7`, then `1`.
+
+No other endings are allowed or reordered.
+
+## Algorithm Steps
+
+1. Accept `R = ⌊√n⌋` and `nMod10` as parameters computed by the caller to avoid redundant square-root and modulo work and to enable reuse in other calculations.
+2. **Group A constants**: test `{3, 7, 11, 13}` when they are `≤ R`. Any hit marks `n` as composite and terminates the search.
+3. **Group A wheel**: generate candidates using the wheel-30 sequence, filter to `d ≡ 3 (mod 10)`, stop at `R`. Any hit marks `n` as composite and terminates the search. 
+4. **Group B**: read the precomputed `n mod 10` (made available by the caller configuration such as `--gpu-prime-batch`). Iterate candidates via wheel 210, maintaining ascending order and applying the ending-priority sequence above. Any hit marks `n` as composite and terminates the search.
+5. **Decision**: if no divisors were found up to `R`, report `n` as “prime by trial”. This outcome is final—no further confirmation runs are performed.
+
+> **Temporary status:** Current builds disable Group B via `UseHeuristicGroupBTrialDivision`. After exhausting Group A, `HeuristicIsPrime*` delegates to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, preserving the Group B implementation for future use while ensuring production checks mirror the external library once Group A filters pass.
+
+## Default Prime Tester Entry Points
+
+* `PrimeTester.IsPrime` and `PrimeTester.IsPrimeGpu` now delegate to the heuristic CPU/GPU implementations by default.
+* Set `PrimeTester.EnableHeuristicPrimeTesting = false` or instantiate `new PrimeTester(useInternal: true)` to force the legacy trial-division and GPU-sieve logic without removing the heuristic helpers.
+* Benchmark harnesses can call `PrimeTester.LegacyIsPrimeInternal` (and `PrimeTester(useInternal: true).IsPrimeGpu`) to compare the legacy behaviour with the heuristic flow when collecting performance baselines.
+
+## Integration with the By-Divisor Pipeline
+
+The by-divisor executors (`MersenneNumberDivisorByDivisorCpuTester` and `MersenneNumberDivisorByDivisorGpuTester`) reuse prime heuristics to seed their divisor scans, residue sieves, and cycle computations. The heuristic implementations must therefore expose the following data and behaviours:
+
+* **Residue filters** – The CPU scan maintains rolling residues modulo 10/8/5/3/7/11 to match the filters in `MersenneCpuDivisorScanSession.CheckDivisors`. `HeuristicIsPrime*` must either share the residue update helpers or produce identical stepping deltas so the divisor scanners can reuse the pre-filtered candidate stream without recomputing moduli. The shared `MersenneDivisorResidueStepper` now centralises those modulo calculations so both paths reuse the same state transitions, and the GPU by-divisor tester instantiates the same stepper before filtering candidates.
+* **Stepper integration** – CPU sessions advance residues and unity checks through `ExponentRemainderStepperCpu` together with `CycleRemainderStepper`; GPU kernels mirror this contract via `ExponentRemainderStepperGpu`. The heuristic implementation must call into these steppers (or emit the same remainder deltas) so the shared stepping state remains aligned across CPU/GPU scans and factor validation.
+* **Wheel alignment** – The wheel 30/210 sequences emitted by `PrimeTester.CreateHeuristicDivisorEnumerator` now drive both the CPU and GPU by-divisor testers, keeping the divisor cache and kernel inputs in lockstep without bespoke residue loops.
+* **Divisor streaming API** – `PrimeTester.CreateHeuristicDivisorEnumerator` now feeds `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` and `MersenneNumberDivisorByDivisorGpuTester.CheckDivisors`, replacing bespoke residue loops while supplying candidates directly to Montgomery preparation and GPU kernels.
+* **GPU batching** – `HeuristicIsPrimeGpu` streams Group A/B candidates through the same enumerator, packs them into batches controlled by `PrimeTester.HeuristicGpuDivisorBatchSize`, and launches `PrimeTesterKernels.HeuristicTrialDivisionKernel` to test each divisor with standard `%` arithmetic. When GPUs are disabled, the method delegates to the CPU heuristic so fallbacks reuse the same candidate ordering without additional plumbing.
+* **Cycle preparation** – When the heuristic discovers a divisor candidate, it should populate `MontgomeryDivisorData` and expose the computed cycle length (or the need to recompute it) so that `DivisorCycleCache` and `MersenneDivisorCycles` do not duplicate work for the same `(prime, divisor)` pair during factor verification. `PrimeTester.PrepareHeuristicDivisor` together with `PrimeTester.ResolveHeuristicCycleLength` now provides this payload to both CPU and GPU by-divisor scans, and the trial-division summaries capture the data whenever a divisor is found.
+
+* **Caching contract** – The by-divisor GPU implementation pools `DivisorScanSession` instances and reuses GPU buffers. The heuristic must accept externally allocated spans/buffers for divisor candidates, cycle lengths, and hits so the existing resource pools remain valid.
+
+## Statistical Guidance
+
+The following summaries cover the first 500,000 odd composite numbers after removing multiples of five (≈381,651 values). Among those, roughly 67,190 numbers had **no** divisor in Group A.
+
+### Distribution of the smallest Group B divisor (when Group A contains a divisor)
+
+* `n ≡ 1 (mod 10)` → ending `7`: **40.30%**, ending `9`: **38.38%**, ending `1`: **21.32%`. In these cases, `min(B) < A_max(n)` occurred **60.77%** of the time, while `min(B) > A_max(n)` occurred **39.23%** of the time.
+* `n ≡ 3 (mod 10)` → ending `9`: **40.91%**, ending `1`: **38.92%**, ending `7`: **20.17%`. Here, `min(B) < A_max(n)` occurred **54.92%** of the time, while `min(B) > A_max(n)` occurred **45.08%** of the time.
+* `n ≡ 7 (mod 10)` → ending `9`: **56.98%**, ending `7`: **21.52%**, ending `1`: **21.50%`. The relation `min(B) < A_max(n)` held **56.09%** of the time versus **43.91%** for `min(B) > A_max(n)`.
+* `n ≡ 9 (mod 10)` → ending `9`: **49.64%**, ending `7`: **27.61%**, ending `1`: **22.75%`. In this class, `min(B) < A_max(n)` occurred **76.20%** of the time, with **23.80%** for `min(B) > A_max(n)`.
+
+### Behaviour when no Group A divisor exists
+
+* With full Group A processing, the minimum Group B divisor satisfied `min(B) < A_max(n)` about **99.67%** of the time, `min(B) = A_max(n)` essentially never, and `min(B) > A_max(n)` about **0.33%** of the time.
+* The ratio `r = min(B) / A_max(n)` showed a median of ≈0.115, with quartiles `Q1 ≈ 0.044` and `Q3 ≈ 0.335`, confirming that the first Group B hit typically appears well before the tail of the Group A sweep.
+* Minimum Group B divisors appeared **only** with the endings listed in the “Allowed endings” table above once Group A candidates were excluded.
+
+### Heuristic cheat sheet
+
+* Group A: test `3, 7, 11, 13` and every `d ≡ 3 (mod 10)` up to `⌊√n⌋` (wheel 30).
+* Group B: test only the allowed endings determined by `n mod 10`, scanning in ascending order (wheel 210) and applying the per-class priority order: `n ≡ 1` → `9` then `1`; `n ≡ 3` → `9` then `7`; `n ≡ 7` → `7` then `1`; `n ≡ 9` → `9`, `7`, `1`.
+* Because `min(B)` is overwhelmingly smaller than `A_max(n)` for the “no Group A divisor” cohort, always start Group B from its smallest candidates instead of resuming above `A_max(n)`.
+
+## Benchmark Reporting
+
+Benchmark results comparing `PrimeTester.IsPrimeGpu`, `PrimeTester.IsPrimeInternal`, and `Open.Numeric.Primes.Prime.Numbers.IsPrime` will be documented here once the dedicated benchmarks covering ≤100, ≤4,000,000, and ≥138,000,000 candidates are executed. The specification will be updated with the measured timings and analysis so operators can evaluate heuristic rollouts against legacy and external implementations. Benchmark harnesses are available via `EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs`, which runs the heuristic CPU and GPU paths alongside `Open.Numeric.Primes` for the three target ranges.
+