Restore CPU PrimeTester usage and rename fallback helpers

EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs

@@ -12,8 +12,8 @@ namespace EvenPerfectBitScanner.Benchmarks;
 [MemoryDiagnoser]
 public class PrimeTesterBenchmarks
 {
-    private PrimeTester _cpuTester = null!;
-    private PrimeTester _gpuTester = null!;
+    private HeuristicPrimeTester _cpuTester = null!;
+    private HeuristicPrimeTester _gpuTester = null!;
     private ulong[] _candidates = [];
 
     [ParamsSource(nameof(GetBenchmarkCases))]
@@ -22,9 +22,8 @@ public class PrimeTesterBenchmarks
     [GlobalSetup]
     public void GlobalSetup()
     {
-        PrimeTester.EnableHeuristicPrimeTesting = true;
-        _cpuTester = new PrimeTester();
-        _gpuTester = new PrimeTester();
+        _cpuTester = new HeuristicPrimeTester();
+        _gpuTester = new HeuristicPrimeTester();
         _candidates = BenchmarkCase.Candidates;
 
         if (_candidates.Length == 0)
@@ -41,7 +40,7 @@ public class PrimeTesterBenchmarks
     [Benchmark(Baseline = true)]
     public int HeuristicCpu()
     {
-        PrimeTester tester = _cpuTester;
+        HeuristicPrimeTester tester = _cpuTester;
         ulong[] values = _candidates;
         CancellationToken cancellationToken = CancellationToken.None;
         int primeCount = 0;
@@ -60,14 +59,14 @@ public class PrimeTesterBenchmarks
     [Benchmark]
     public int HeuristicGpu()
     {
-        PrimeTester tester = _gpuTester;
+        HeuristicPrimeTester 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]))
+            if (tester.IsPrimeGpu(values[i]))
             {
                 primeCount++;
             }
@@ -87,7 +86,7 @@ public class PrimeTesterBenchmarks
 
         for (int i = 0; i < values.Length; i++)
         {
-            if (PrimeTester.LegacyIsPrimeGpu(values[i], cancellationToken))
+            if (PrimeTester.IsPrimeGpuFallback(values[i], cancellationToken))
             {
                 primeCount++;
             }
@@ -106,7 +105,7 @@ public class PrimeTesterBenchmarks
 
         for (int i = 0; i < values.Length; i++)
         {
-            if (PrimeTester.LegacyIsPrimeInternal(values[i], cancellationToken))
+            if (PrimeTester.IsPrimeInternal(values[i], cancellationToken))
             {
                 primeCount++;
             }

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -106,13 +106,17 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         processedAll = false;
         processedCount = 0UL;
 
+        // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
         var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
 
+        // while (enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate))
         while (enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
         {
             processedCount = enumerator.ProcessedCount;
 
+            // HeuristicPrimeTester.HeuristicDivisorPreparation preparation = HeuristicPrimeTester.PrepareHeuristicDivisor(in candidate);
             PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate);
+            // ulong divisorCycle = HeuristicPrimeTester.ResolveHeuristicCycleLength(
             ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength(
                 prime,
                 in preparation,

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -189,6 +189,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         processedCount = 0UL;
         lastProcessed = 0UL;
 
+        // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
         var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
 
         Span<ulong> divisorSpan = divisors.AsSpan();
@@ -203,12 +204,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             int admissibleCount = 0;
             int exponentIndex = 0;
 
+            // while (admissibleCount < batchCapacity && enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate))
             while (admissibleCount < batchCapacity && enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
             {
                 processedCount = enumerator.ProcessedCount;
                 ulong divisorValue = candidate.Value;
 
+                // HeuristicPrimeTester.HeuristicDivisorPreparation preparation = HeuristicPrimeTester.PrepareHeuristicDivisor(in candidate);
                 PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate);
+                // ulong divisorCycle = HeuristicPrimeTester.ResolveHeuristicCycleLength(
                 ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength(
                     prime,
                     in preparation,

PerfectNumbers.Core/HeuristicPrimeTester.cs

@@ -0,0 +1,968 @@
+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 System.Threading;
+using ILGPU;
+using ILGPU.Algorithms;
+using ILGPU.Runtime;
+using PerfectNumbers.Core.Gpu;
+
+namespace PerfectNumbers.Core;
+
+public sealed class HeuristicPrimeTester
+{
+    static HeuristicPrimeTester()
+    {
+        HeuristicPrimeSieves.EnsureInitialized();
+    }
+
+    // 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];
+
+    // 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 static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B.
+
+
+    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 static bool EnableHeuristicPrimeTesting { get; set; } = true;
+
+    public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096;
+
+    public static int HeuristicDivisorInterleaveBatchSize { get; set; } = 64;
+
+    [ThreadStatic]
+    private static HeuristicPrimeTester? _tester;
+
+    public static HeuristicPrimeTester Exclusive => _tester ??= new();
+
+    public bool IsPrime(ulong n, CancellationToken ct)
+    {
+        if (!EnableHeuristicPrimeTesting)
+        {
+            return PrimeTester.IsPrimeInternal(n, ct);
+        }
+
+        return HeuristicIsPrimeCpu(n);
+    }
+
+    public bool IsPrime(ulong n)
+    {
+        return IsPrime(n, CancellationToken.None);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public bool IsPrimeGpu(ulong n, CancellationToken ct)
+    {
+        if (!EnableHeuristicPrimeTesting)
+        {
+            return PrimeTester.IsPrimeGpuFallback(n, ct);
+        }
+
+        return IsPrimeGpu(n);
+    }
+
+    public bool IsPrimeGpu(ulong n)
+    {
+        byte nMod10 = (byte)n.Mod10();
+        ulong sqrtLimit = ComputeHeuristicSqrt(n);
+        return IsPrimeGpu(n, sqrtLimit, nMod10);
+    }
+
+    public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10)
+    {
+        if (!EnableHeuristicPrimeTesting)
+        {
+            return PrimeTester.IsPrimeGpuFallback(n, CancellationToken.None);
+        }
+
+        return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10);
+    }
+
+    public bool HeuristicIsPrimeGpu(ulong n)
+    {
+        return IsPrimeGpu(n);
+    }
+
+    public static bool HeuristicIsPrimeCpu(ulong n)
+    {
+        return HeuristicIsPrimeCpu(n, 0UL, (byte)n.Mod10());
+    }
+
+    public static bool HeuristicIsPrimeCpu(ulong n, ulong sqrtLimit, byte nMod10)
+    {
+        if (!EnableHeuristicPrimeTesting)
+        {
+            return PrimeTester.IsPrimeInternal(n, CancellationToken.None);
+        }
+
+        if (sqrtLimit == 0UL)
+        {
+            sqrtLimit = ComputeHeuristicSqrt(n);
+        }
+
+        if (sqrtLimit < 3UL)
+        {
+            return EvaluateWithOpenNumericFallback(n);
+        }
+
+        bool includeGroupB = UseHeuristicGroupBTrialDivision;
+        return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB);
+    }
+
+    private bool HeuristicIsPrimeGpuCore(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;
+    }
+
+    private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB)
+    {
+        ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
+        int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize);
+        int groupAIndex = 0;
+
+        if (!includeGroupB)
+        {
+            while (groupAIndex < groupADivisors.Length)
+            {
+                int processed = 0;
+                while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
+                {
+                    ulong divisor = (ulong)groupADivisors[groupAIndex];
+                    if (divisor > sqrtLimit)
+                    {
+                        groupAIndex = groupADivisors.Length;
+                        break;
+                    }
+
+                    if (n % divisor == 0UL)
+                    {
+                        return false;
+                    }
+
+                    groupAIndex++;
+                    processed++;
+                }
+
+                if (groupAIndex >= groupADivisors.Length)
+                {
+                    break;
+                }
+
+                if ((ulong)groupADivisors[groupAIndex] > sqrtLimit)
+                {
+                    groupAIndex = groupADivisors.Length;
+                    break;
+                }
+            }
+
+            return EvaluateWithOpenNumericFallback(n);
+        }
+
+        ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
+        if (endingOrder.IsEmpty)
+        {
+            return EvaluateWithOpenNumericFallback(n);
+        }
+
+        Span<int> indices = endingOrder.Length <= 8
+            ? stackalloc int[endingOrder.Length]
+            : new int[endingOrder.Length];
+
+        bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
+        bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit);
+
+        while (groupAHasMore || groupBHasMore)
+        {
+            if (groupAHasMore)
+            {
+                int processed = 0;
+                while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
+                {
+                    ulong divisor = (ulong)groupADivisors[groupAIndex];
+                    if (divisor > sqrtLimit)
+                    {
+                        groupAIndex = groupADivisors.Length;
+                        groupAHasMore = false;
+                        break;
+                    }
+
+                    if (n % divisor == 0UL)
+                    {
+                        return false;
+                    }
+
+                    groupAIndex++;
+                    processed++;
+                }
+
+                if (groupAIndex >= groupADivisors.Length)
+                {
+                    groupAHasMore = false;
+                }
+                else
+                {
+                    groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
+                    if (!groupAHasMore)
+                    {
+                        groupAIndex = groupADivisors.Length;
+                    }
+                }
+            }
+
+            if (groupBHasMore)
+            {
+                int processed = 0;
+                while (processed < interleaveBatchSize)
+                {
+                    if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor))
+                    {
+                        groupBHasMore = false;
+                        break;
+                    }
+
+                    if (n % divisor == 0UL)
+                    {
+                        return false;
+                    }
+
+                    processed++;
+                }
+
+                if (groupBHasMore)
+                {
+                    groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit);
+                }
+            }
+        }
+
+        return true;
+    }
+
+
+
+    private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10)
+    {
+        int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize);
+        int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize);
+        var divisorPool = ThreadStaticPools.UlongPool;
+        var hitPool = ThreadStaticPools.BytePool;
+
+        ulong[]? divisorArray = null;
+        byte[]? hitFlags = null;
+
+        var limiter = GpuPrimeWorkLimiter.Acquire();
+        var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent();
+
+        bool compositeDetected = false;
+
+        try
+        {
+            divisorArray = divisorPool.Rent(batchCapacity);
+            hitFlags = hitPool.Rent(batchCapacity);
+
+            var accelerator = gpu.Accelerator;
+            var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator);
+
+            lock (gpu.ExecutionLock)
+            {
+                var scratch = state.RentScratch(batchCapacity, accelerator);
+                try
+                {
+                    int count = 0;
+
+                    bool ProcessBatch(int length)
+                    {
+                        if (length <= 0)
+                        {
+                            return false;
+                        }
+
+                        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 = divisorArray![i];
+                            if (divisor > 1UL && n % divisor == 0UL)
+                            {
+                                return true;
+                            }
+                        }
+
+                        return false;
+                    }
+
+                    bool QueueDivisor(ulong divisor)
+                    {
+                        divisorArray![count] = divisor;
+                        count++;
+
+                        if (count == batchCapacity)
+                        {
+                            if (ProcessBatch(count))
+                            {
+                                return true;
+                            }
+
+                            count = 0;
+                        }
+
+                        return false;
+                    }
+
+                    ReadOnlySpan<int> groupADivisors = HeuristicPrimeSieves.GroupADivisors;
+                    int groupAIndex = 0;
+                    ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
+                    Span<int> indices = endingOrder.Length <= 8
+                        ? stackalloc int[endingOrder.Length]
+                        : new int[endingOrder.Length];
+
+                    bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
+                    bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, sqrtLimit);
+
+                    while (!compositeDetected && (groupAHasMore || groupBHasMore))
+                    {
+                        if (groupAHasMore)
+                        {
+                            int processed = 0;
+                            while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length)
+                            {
+                                ulong divisor = (ulong)groupADivisors[groupAIndex];
+                                if (divisor > sqrtLimit)
+                                {
+                                    groupAIndex = groupADivisors.Length;
+                                    groupAHasMore = false;
+                                    break;
+                                }
+
+                                compositeDetected = QueueDivisor(divisor);
+                                if (compositeDetected)
+                                {
+                                    break;
+                                }
+
+                                groupAIndex++;
+                                processed++;
+                            }
+
+                            if (compositeDetected)
+                            {
+                                break;
+                            }
+
+                            if (groupAIndex >= groupADivisors.Length)
+                            {
+                                groupAHasMore = false;
+                            }
+                            else
+                            {
+                                groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit;
+                                if (!groupAHasMore)
+                                {
+                                    groupAIndex = groupADivisors.Length;
+                                }
+                            }
+                        }
+
+                        if (groupBHasMore)
+                        {
+                            int processed = 0;
+                            while (processed < interleaveBatchSize)
+                            {
+                                if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor))
+                                {
+                                    groupBHasMore = false;
+                                    break;
+                                }
+
+                                compositeDetected = QueueDivisor(divisor);
+                                if (compositeDetected)
+                                {
+                                    break;
+                                }
+
+                                processed++;
+                            }
+
+                            if (compositeDetected)
+                            {
+                                break;
+                            }
+
+                            if (groupBHasMore)
+                            {
+                                groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit);
+                            }
+                        }
+                    }
+
+                    if (!compositeDetected && count > 0)
+                    {
+                        compositeDetected = ProcessBatch(count);
+                    }
+                }
+                finally
+                {
+                    state.ReturnScratch(scratch);
+                }
+            }
+        }
+        finally
+        {
+            if (divisorArray is not null)
+            {
+                divisorPool.Return(divisorArray);
+            }
+
+            if (hitFlags is not null)
+            {
+                hitPool.Return(hitFlags);
+            }
+
+            gpu.Dispose();
+            limiter.Dispose();
+        }
+
+        return compositeDetected;
+    }
+
+
+
+    private static bool EvaluateWithOpenNumericFallback(ulong n)
+    {
+        return Prime.Numbers.IsPrime(n);
+    }
+
+    internal static ulong ComputeHeuristicSqrt(ulong n)
+    {
+        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,
+    };
+
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static bool TrySelectNextGroupBDivisor(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong sqrtLimit, out ulong divisor)
+    {
+        ulong bestCandidate = ulong.MaxValue;
+        int bestEndingIndex = -1;
+
+        for (int i = 0; i < endingOrder.Length; i++)
+        {
+            ReadOnlySpan<int> divisors = GetGroupBDivisors(endingOrder[i]);
+            int index = indices[i];
+            if ((uint)index >= (uint)divisors.Length)
+            {
+                continue;
+            }
+
+            ulong candidate = (ulong)divisors[index];
+            if (candidate > sqrtLimit)
+            {
+                indices[i] = divisors.Length;
+                continue;
+            }
+
+            if (bestEndingIndex == -1 || candidate < bestCandidate)
+            {
+                bestCandidate = candidate;
+                bestEndingIndex = i;
+            }
+        }
+
+        if (bestEndingIndex == -1)
+        {
+            divisor = 0UL;
+            return false;
+        }
+
+        divisor = bestCandidate;
+        indices[bestEndingIndex]++;
+        return true;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static bool HasGroupBCandidates(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong sqrtLimit)
+    {
+        for (int i = 0; i < endingOrder.Length; i++)
+        {
+            ReadOnlySpan<int> divisors = GetGroupBDivisors(endingOrder[i]);
+            int index = indices[i];
+            if ((uint)index >= (uint)divisors.Length)
+            {
+                continue;
+            }
+
+            ulong candidate = (ulong)divisors[index];
+            if (candidate > sqrtLimit)
+            {
+                indices[i] = divisors.Length;
+                continue;
+            }
+
+            return true;
+        }
+
+        return false;
+    }
+
+    private static ReadOnlySpan<int> GetGroupBDivisors(byte ending) => ending switch
+    {
+        1 => HeuristicPrimeSieves.GroupBDivisorsEnding1,
+        3 => HeuristicPrimeSieves.GroupBDivisorsEnding3,
+        7 => HeuristicPrimeSieves.GroupBDivisorsEnding7,
+        9 => HeuristicPrimeSieves.GroupBDivisorsEnding9,
+        _ => ReadOnlySpan<int>.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);
+        }
+    }
+}

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -395,6 +395,7 @@ public class MersenneDivisorCycles
         }
 
         // TODO: Benchmark PrimeTester.IsPrimeInternal against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites.
+        // if (HeuristicPrimeTester.Exclusive.IsPrime(remaining, CancellationToken.None))
         if (PrimeTester.IsPrimeInternal(remaining, CancellationToken.None))
         {
             AddFactor(counts, remaining);
@@ -927,6 +928,7 @@ public class MersenneDivisorCycles
             return true;
         }
 
+        // if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrime(divisor, CancellationToken.None))
         if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
         {
             ulong computedOrder = PrimeOrderCalculator.Calculate(

PerfectNumbers.Core/MersennePrimeFactorTester.cs

@@ -25,6 +25,7 @@ public static class MersennePrimeFactorTester
             return false;
         }
 
+        // if (!HeuristicPrimeTester.Exclusive.IsPrime(q, ct))
         if (!PrimeTester.IsPrimeInternal(q, ct))
         {
             return false;
@@ -233,6 +234,7 @@ public static class MersennePrimeFactorTester
                 continue;
             }
 
+            // if (HeuristicPrimeTester.Exclusive.IsPrime(m, ct))
             if (PrimeTester.IsPrimeInternal(m, ct))
             {
                 dict[m] = 1;

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs

@@ -1124,7 +1124,7 @@ internal static partial class PrimeOrderCalculator
 
 					// bool isPrime = primeTester.HeuristicIsPrimeGpu(composite);
 					// bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
-					bool isPrime = PrimeTester.LegacyIsPrimeInternal(composite, CancellationToken.None);
+					bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
 					// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
 
 					if (isPrime)
@@ -1187,7 +1187,7 @@ internal static partial class PrimeOrderCalculator
 			{
 				// bool isPrime = primeTester.HeuristicIsPrimeGpu(composite);
 				// bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
-				bool isPrime = PrimeTester.LegacyIsPrimeInternal(composite, CancellationToken.None);
+				bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None);
 				// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
 				
 				entry = entry.WithPrimality(isPrime);
@@ -1220,7 +1220,7 @@ internal static partial class PrimeOrderCalculator
 		{
 			// cofactorIsPrime = primeTester.HeuristicIsPrimeGpu(cofactor);
 			// cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None);
-			cofactorIsPrime = PrimeTester.LegacyIsPrimeInternal(cofactor, CancellationToken.None);
+			cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None);
 			// cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor);
 		}
 
@@ -1753,7 +1753,7 @@ internal static partial class PrimeOrderCalculator
 
 		// bool isPrime = primeTester.HeuristicIsPrimeGpu(value);
 		// bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None);
-		bool isPrime = PrimeTester.LegacyIsPrimeInternal(value, CancellationToken.None);
+		bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None);
 		// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value);
 		if (!knownComposite && isPrime)
 		{
@@ -1766,7 +1766,7 @@ internal static partial class PrimeOrderCalculator
 
 		// isPrime = primeTester.HeuristicIsPrimeGpu(factor);
 		// isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None);
-		isPrime = PrimeTester.LegacyIsPrimeInternal(factor, CancellationToken.None);
+		isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None);
 		// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor);
 		if (isPrime)
 		{
@@ -1798,7 +1798,7 @@ internal static partial class PrimeOrderCalculator
 		{
 			// isPrime = primeTester.HeuristicIsPrimeGpu(quotient);
 			// isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None);
-			isPrime = PrimeTester.LegacyIsPrimeInternal(quotient, CancellationToken.None);
+			isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None);
 			// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient);
 		}
 

PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs

@@ -955,6 +955,7 @@ internal static partial class PrimeOrderCalculator
     {
         if (value <= ulong.MaxValue)
         {
+            // return HeuristicPrimeTester.Exclusive.IsPrime((ulong)value, CancellationToken.None);
             return PrimeTester.IsPrimeInternal((ulong)value, CancellationToken.None);
         }
 

docs/HeuristicIsPrimePlan.md

@@ -109,7 +109,7 @@ internal static bool IsPrimeInternal(ulong n, CancellationToken ct)
 }
 ```
 
-## 3. Implementation Plan for `PrimeTester.HeuristicIsPrime*`
+## 3. Implementation Plan for `HeuristicPrimeTester`
 
 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]`.
 
@@ -117,12 +117,12 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels,
 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] Added `HeuristicPrimeSieves` to precompute 4M-entry Group A/B divisor tables so heuristic prime checks reuse cached sequences without regenerating wheel steps at runtime; the enumerator remains available for residue-driven consumers.
-   * [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] Surface helpers that accept the precomputed `R = ⌊√n⌋` and `nMod10` values from callers via the new overload of `HeuristicIsPrimeCpu`; HeuristicPrimeTester 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.
+   * [done] Publish wheel iterators that can stream candidates directly into `DivisorCycleCache`, GPU hit buffers, or CPU Montgomery checks without recomputing offsets by exposing `HeuristicPrimeTester.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.
@@ -136,16 +136,16 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels,
    * [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] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.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.
+   * [done] Introduce `HeuristicPrimeTester` so heuristic CPU/GPU checks live alongside the legacy `PrimeTester` implementations while callers opt in explicitly.
+   * [done] Maintain legacy paths for benchmarking via `PrimeTester.IsPrimeInternal`, `PrimeTester.IsPrimeGpuFallback`, and the global `HeuristicPrimeTester.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.
+    * [done] Build BenchmarkDotNet comparisons for `HeuristicPrimeTester.IsPrimeGpu`, `HeuristicPrimeTester.HeuristicIsPrimeCpu`, 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.

docs/HeuristicIsPrimeSpecification.md

@@ -1,6 +1,6 @@
 # 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.
+This document describes the heuristic prime testing procedure implemented by `HeuristicPrimeTester.HeuristicIsPrimeCpu` and `HeuristicPrimeTester.IsPrimeGpu`. 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
 
@@ -61,9 +61,9 @@ No other endings are allowed or reordered.
 
 ## 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.
+* `PrimeTester.IsPrime` and `PrimeTester.IsPrimeGpu` continue to run the legacy trial-division and GPU sieve paths. Use `HeuristicPrimeTester` when the heuristic CPU/GPU flow is required.
+* Toggle `HeuristicPrimeTester.EnableHeuristicPrimeTesting` or fall back to `PrimeTester.IsPrimeInternal` / `PrimeTester.IsPrimeGpuFallback` to control heuristic rollouts in test harnesses.
+* Benchmark harnesses can compare behaviours by calling the legacy helpers on `PrimeTester` alongside the heuristic equivalents on `HeuristicPrimeTester`.
 
 ## Integration with the By-Divisor Pipeline
 
@@ -72,9 +72,9 @@ The by-divisor executors (`MersenneNumberDivisorByDivisorCpuTester` and `Mersenn
 * **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** – `HeuristicPrimeSieves` precomputes the wheel 30/210 sequences once at startup so both the heuristic prime tester and the by-divisor executors consume identical Group A/B streams without regenerating residues on every call.
-* **Divisor streaming API** – `PrimeTester.CreateHeuristicDivisorEnumerator` exposes the same precomputed sieves so by-divisor scans can pull candidates directly into Montgomery preparation and GPU kernels without bespoke residue loops.
-* **GPU batching** – `HeuristicIsPrimeGpu` streams Group A/B candidates from `HeuristicPrimeSieves`, 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.
+* **Divisor streaming API** – `HeuristicPrimeTester.CreateHeuristicDivisorEnumerator` exposes the same precomputed sieves so by-divisor scans can pull candidates directly into Montgomery preparation and GPU kernels without bespoke residue loops.
+* **GPU batching** – `HeuristicIsPrimeGpu` streams Group A/B candidates from `HeuristicPrimeSieves`, packs them into batches controlled by `HeuristicPrimeTester.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. `HeuristicPrimeTester.PrepareHeuristicDivisor` together with `HeuristicPrimeTester.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.
 
@@ -103,5 +103,5 @@ The following summaries cover the first 500,000 odd composite numbers after remo
 
 ## 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.
+Benchmark results comparing `HeuristicPrimeTester.IsPrimeGpu`, `HeuristicPrimeTester.HeuristicIsPrimeCpu`, 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.