Restore legacy by-divisor scans for Mersenne testers

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -1,5 +1,6 @@
 using System.Collections.Generic;
 using System.Runtime.CompilerServices;
+using PerfectNumbers.Core;
 using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core.Cpu;
@@ -12,6 +13,16 @@ 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;
@@ -103,41 +114,152 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         out bool processedAll,
         out ulong processedCount)
     {
-        processedAll = false;
         processedCount = 0UL;
+        processedAll = false;
 
-        // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
-        var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
+        // 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;
+        // }
 
-        // while (enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate))
-        while (enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
+        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)
         {
-            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,
-                out _,
-                out _,
-                out _);
-
-            if (divisorCycle == prime)
-            {
-                processedAll = true;
-                return true;
-            }
-
-            if (divisorCycle == 0UL)
-            {
-                Console.WriteLine($"Divisor cycle was not calculated for {prime}");
-            }
+            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;
-        processedAll = enumerator.Exhausted;
+        // 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;
+
+        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.
+        LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
+        bool lastIsSeven = lastDigit == LastDigit.Seven;
+        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)
+            {
+                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}");
+                }
+            }
+
+            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);
+        }
+
+        processedAll = divisor.CompareTo(limit) > 0;
+        divisor.High = 0UL;
+        divisor.Low = 0UL;
+        step.High = 0UL;
+        step.Low = 0UL;
+        limit.High = 0UL;
+        limit.Low = 0UL;
         return false;
     }
 
@@ -147,6 +269,17 @@ 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)

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -112,6 +112,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             resources.HitIndexBuffer,
             resources.Divisors,
             resources.Exponents,
+            resources.FilteredDivisors,
             resources.DivisorData,
             resources.Offsets,
             resources.Counts,
@@ -176,6 +177,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         MemoryBuffer1D<int, Stride1D.Dense> hitIndexBuffer,
         ulong[] divisors,
         ulong[] exponents,
+        ulong[] filteredDivisors,
         GpuDivisorPartialData[] divisorData,
         int[] offsets,
         int[] counts,
@@ -186,40 +188,131 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
     {
         int batchCapacity = (int)divisorDataBuffer.Length;
         bool composite = false;
+        bool processedAll = false;
         processedCount = 0UL;
         lastProcessed = 0UL;
 
-        // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
-        var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax);
+        int chunkCapacity = Math.Max(1, batchCapacity);
 
-        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)
+        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)
         {
+            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));
+        LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One;
+
+        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, lastDigit))
+                {
+                    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;
-
-            // while (admissibleCount < batchCapacity && enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate))
-            while (admissibleCount < batchCapacity && enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate))
+            for (int i = 0; i < filteredCount; i++)
             {
-                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,
-                    out _,
-                    out _,
-                    out _);
-
+                ulong divisorValue = filteredDivisorsSpan[i];
+                MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue);
+                ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData);
                 if (divisorCycle != prime)
                 {
                     continue;
@@ -236,16 +329,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
                 exponentIndex++;
             }
 
-            processedCount = enumerator.ProcessedCount;
-            lastProcessed = enumerator.LastDivisor;
-
             if (admissibleCount == 0)
             {
-                if (enumerator.Exhausted)
-                {
-                    break;
-                }
-
                 continue;
             }
 
@@ -277,10 +362,74 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex];
         }
 
-        coveredRange = composite || enumerator.Exhausted;
+        coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128);
         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,
+        LastDigit lastDigit)
+    {
+        const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
+        const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
+
+        ushort decimalMask = lastDigit == LastDigit.Seven ? 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,
@@ -571,6 +720,8 @@ 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!;
@@ -645,6 +796,7 @@ 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);
@@ -683,12 +835,14 @@ 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>();