Merge pull request #128 from mlemanczyk/codex/modify-evenperfectbitscanner-execution-path

Remove CPU fallbacks when GPU pow2 fails

EvenPerfectBitScanner/Program.cs

@@ -515,7 +515,9 @@ internal static class Program
         {
             _byDivisorTester = mersenneOnGpu
                     ? new MersenneNumberDivisorByDivisorGpuTester()
-                    : new MersenneNumberDivisorByDivisorCpuTester();
+                    : new MersenneNumberDivisorByDivisorCpuTester(
+                        useGpuPow2Mod: !GpuContextPool.ForceCpu,
+                        useGpuCycleComputation: orderOnGpu);
             _byDivisorTester.BatchSize = scanBatchSize;
         }
 

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -3,6 +3,7 @@ using System.Collections.Concurrent;
 using System.Collections.Generic;
 using System.Numerics;
 using PerfectNumbers.Core;
+using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core.Cpu;
 
@@ -10,11 +11,18 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 {
     private readonly object _sync = new();
     private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
+    private readonly bool _useGpuPow2Mod;
+    private readonly bool _useGpuCycleComputation;
     private ulong _divisorLimit;
     private ulong _lastStatusDivisor;
     private bool _isConfigured;
     private int _batchSize = 1_024;
 
+    public MersenneNumberDivisorByDivisorCpuTester(bool useGpuPow2Mod = false, bool useGpuCycleComputation = false)
+    {
+        _useGpuPow2Mod = useGpuPow2Mod;
+        _useGpuCycleComputation = useGpuCycleComputation;
+    }
 
     public int BatchSize
     {
@@ -80,9 +88,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             return true;
         }
 
-        ulong processedCount = 0UL;
+        ulong processedCount;
-        ulong lastProcessed = 0UL;
+        ulong lastProcessed;
-        bool processedAll = false;
+        bool processedAll;
 
         bool composite = CheckDivisors(
             prime,
@@ -190,6 +198,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         }
 
         Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
+        Dictionary<ulong, ulong>? cycleRemainderCache = null;
+        CycleHint cycleHint = default;
         DivisorCycleCache cycleCache = DivisorCycleCache.Shared;
 
         byte step10 = (byte)(step % 10UL);
@@ -223,34 +233,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0)
             {
                 MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(candidate);
-                ulong divisorCycle;
+                if (TryResolveDivisorCycle(candidate, prime, divisorData, cycleCache, ref factorCache, ref cycleHint, out ulong divisorCycle) && divisorCycle != 0UL)
-
-                if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles)
                 {
-                    divisorCycle = cycleCache.GetCycleLength(candidate);
+                    if (CheckDivisor(prime, candidate, divisorCycle, divisorData, ref cycleRemainderCache, ref cycleHint) != 0)
-                }
-                else
-                {
-                    factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
-                    if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL)
                     {
-                        divisorCycle = cycleCache.GetCycleLength(candidate);
+                        processedAll = true;
+                        return true;
                     }
                 }
-
-                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.
-                    processedAll = true;
-                    return true;
-                }
-
-                if (divisorCycle == 0UL)
-                {
-                    Console.WriteLine($"Divisor cycle was not calculated for {prime}");
-                }
             }
 
             divisor += step;
@@ -266,12 +256,152 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         return false;
     }
 
-    private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
+    private byte CheckDivisor(
+        ulong prime,
+        ulong divisor,
+        ulong divisorCycle,
+        in MontgomeryDivisorData divisorData,
+        ref Dictionary<ulong, ulong>? cycleRemainders,
+        ref CycleHint cycleHint)
     {
-        ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData);
+        if (divisorCycle == 0UL)
+        {
+            throw new InvalidOperationException($"Cycle length resolution returned zero for divisor {divisor}.");
+        }
+
+        Dictionary<ulong, ulong>? remainderCache = cycleRemainders;
+        if (_useGpuPow2Mod && remainderCache is null)
+        {
+            remainderCache = new Dictionary<ulong, ulong>(8);
+            cycleRemainders = remainderCache;
+        }
+
+        ulong reducedExponent;
+        if (remainderCache is not null && remainderCache.TryGetValue(divisorCycle, out ulong cachedRemainder))
+        {
+            reducedExponent = cachedRemainder;
+        }
+        else if (cycleHint.HasCycle && cycleHint.LastCycleLength == divisorCycle)
+        {
+            reducedExponent = cycleHint.LastReducedExponent;
+        }
+        else
+        {
+            reducedExponent = prime % divisorCycle;
+            if (remainderCache is not null)
+            {
+                remainderCache[divisorCycle] = reducedExponent;
+            }
+        }
+
+        cycleHint.LastCycleLength = divisorCycle;
+        cycleHint.LastReducedExponent = reducedExponent;
+        cycleHint.HasCycle = true;
+
+        if (_useGpuPow2Mod && remainderCache is not null)
+        {
+            GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(reducedExponent, divisor, out ulong remainder, divisorData);
+            if (status == GpuPow2ModStatus.Success)
+            {
+                return remainder == 1UL ? (byte)1 : (byte)0;
+            }
+
+            PrimeOrderGpuHeuristics.ReportPow2Failure(
+                "by-divisor residue scan",
+                divisor,
+                status);
+        }
+
+        ulong residue = reducedExponent.Pow2MontgomeryModFromCycleRemainder(divisorData);
         return residue == 1UL ? (byte)1 : (byte)0;
     }
 
+    private bool TryResolveDivisorCycle(
+        ulong divisor,
+        ulong prime,
+        in MontgomeryDivisorData divisorData,
+        DivisorCycleCache cycleCache,
+        ref Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache,
+        ref CycleHint cycleHint,
+        out ulong cycleLength)
+    {
+        if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+        {
+            cycleLength = cycleCache.GetCycleLength(divisor);
+            if (cycleLength == 0UL)
+            {
+                throw new InvalidOperationException($"Missing cached cycle length for divisor {divisor}.");
+            }
+
+            cycleHint.LastCycleLength = cycleLength;
+            cycleHint.HasCycle = true;
+            return true;
+        }
+
+        if (_useGpuCycleComputation && TryCalculateCycleLengthGpu(divisor, prime, divisorData, ref cycleHint, out cycleLength))
+        {
+            return cycleLength != 0UL;
+        }
+
+        factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
+        if (MersenneDivisorCycles.TryCalculateCycleLengthForExponent(
+                divisor,
+                prime,
+                divisorData,
+                factorCache,
+                out cycleLength,
+                allowGpuPow2: _useGpuPow2Mod) && cycleLength != 0UL)
+        {
+            cycleHint.LastCycleLength = cycleLength;
+            cycleHint.HasCycle = true;
+            return true;
+        }
+
+        cycleLength = cycleCache.GetCycleLength(divisor);
+        if (cycleLength == 0UL)
+        {
+            return false;
+        }
+
+        cycleHint.LastCycleLength = cycleLength;
+        cycleHint.HasCycle = true;
+        return true;
+    }
+
+    private bool TryCalculateCycleLengthGpu(
+        ulong divisor,
+        ulong prime,
+        in MontgomeryDivisorData divisorData,
+        ref CycleHint cycleHint,
+        out ulong cycleLength)
+    {
+        PrimeOrderCalculator.PrimeOrderSearchConfig config = PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault;
+        ulong? previousOrder = cycleHint.HasCycle ? cycleHint.LastCycleLength : null;
+
+        try
+        {
+            if (!PrimeOrderGpuHeuristics.TryCalculateOrder(divisor, previousOrder, config, divisorData, out PrimeOrderCalculator.PrimeOrderResult result))
+            {
+                throw new InvalidOperationException("GPU requested fallback.");
+            }
+
+            cycleLength = result.Order;
+            if (cycleLength == 0UL)
+            {
+                return false;
+            }
+
+            cycleHint.LastCycleLength = cycleLength;
+            cycleHint.HasCycle = true;
+            return true;
+        }
+        catch (InvalidOperationException ex)
+        {
+            Console.Error.WriteLine($"GPU cycle length computation failed for divisor {divisor}: {ex.Message}");
+            throw;
+        }
+    }
+
     private static byte AddMod(byte value, byte delta, byte modulus)
     {
         int sum = value + delta;
@@ -298,16 +428,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         }
 
         ulong total = _lastStatusDivisor + processedCount;
-        if (total >= interval)
+        // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
-        {
+        // status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly.
-            do
+        _lastStatusDivisor = total % interval;
-            {
-                total -= interval;
-            }
-            while (total >= interval);
-        }
-
-        _lastStatusDivisor = total;
     }
 
     private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
@@ -325,6 +448,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         return (1UL << (int)(maxPrime - 1UL)) - 1UL;
     }
 
+    private struct CycleHint
+    {
+        public ulong LastCycleLength;
+        public ulong LastReducedExponent;
+        public bool HasCycle;
+    }
+
     private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
     {
         if (prime <= 1UL)

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -4,10 +4,12 @@ using System.Collections.Concurrent;
 using ILGPU;
 using ILGPU.Runtime;
 using System.Numerics;
+using System.Globalization;
 using PerfectNumbers.Core;
 using System.Runtime.InteropServices;
 using System.Runtime.CompilerServices;
 using ILGPU.Runtime.OpenCL;
+
 namespace PerfectNumbers.Core.Gpu;
 
 internal enum GpuPow2ModStatus
@@ -27,6 +29,25 @@ internal static partial class PrimeOrderGpuHeuristics
     private static readonly ConcurrentDictionary<Accelerator, OrderKernelLauncher> OrderKernelCache = new();
     private static readonly ConcurrentDictionary<Accelerator, SmallPrimeDeviceCache> SmallPrimeDeviceCaches = new();
 
+    private static void ReportPow2FailureCore(string context, string modulusDescription, GpuPow2ModStatus status)
+    {
+        string message = $"GPU pow2mod failed with status {status} for divisor {modulusDescription} during {context}.";
+        Console.Error.WriteLine(message);
+        throw new InvalidOperationException(message);
+    }
+
+    internal static void ReportPow2Failure(string context, ulong modulus, GpuPow2ModStatus status)
+    {
+        string modulusDescription = modulus.ToString(CultureInfo.InvariantCulture);
+        ReportPow2FailureCore(context, modulusDescription, status);
+    }
+
+    internal static void ReportPow2Failure(string context, UInt128 modulus, GpuPow2ModStatus status)
+    {
+        string modulusDescription = modulus.ToString(CultureInfo.InvariantCulture);
+        ReportPow2FailureCore(context, modulusDescription, status);
+    }
+
     public readonly struct OrderKernelConfig(ulong previousOrder, byte hasPreviousOrder, uint smallFactorLimit, int maxPowChecks, int mode)
     {
         public readonly ulong PreviousOrder = previousOrder;

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -269,7 +269,8 @@ public class MersenneDivisorCycles
         ulong exponent,
         in MontgomeryDivisorData divisorData,
         Dictionary<ulong, FactorCacheEntry>? factorCache,
-        out ulong cycleLength)
+        out ulong cycleLength,
+        bool allowGpuPow2 = false)
     {
         cycleLength = 0UL;
 
@@ -290,7 +291,8 @@ public class MersenneDivisorCycles
             previousOrder: null,
             divisorData,
             PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
-            PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
+            PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu,
+            allowGpuPow2);
         if (orderResult.Order != 0UL)
         {
             cycleLength = orderResult.Order;
@@ -327,7 +329,7 @@ public class MersenneDivisorCycles
             return false;
         }
 
-        cycleLength = ReduceOrder(divisorData, phi, factorCounts);
+        cycleLength = ReduceOrder(divisorData, phi, factorCounts, allowGpuPow2);
         return true;
     }
 
@@ -471,7 +473,7 @@ public class MersenneDivisorCycles
         }
     }
 
-    private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts)
+    private static ulong ReduceOrder(in MontgomeryDivisorData divisorData, ulong initialOrder, Dictionary<ulong, int> factorCounts, bool allowGpuPow2)
     {
         if (factorCounts.Count == 0)
         {
@@ -505,7 +507,7 @@ public class MersenneDivisorCycles
                     }
 
                     ulong candidate = order / prime;
-                    if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+                    if (Pow2EqualsOne(candidate, divisorData, allowGpuPow2))
                     {
                         order = candidate;
                         continue;
@@ -523,6 +525,26 @@ public class MersenneDivisorCycles
         }
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static bool Pow2EqualsOne(ulong exponent, in MontgomeryDivisorData divisorData, bool allowGpuPow2)
+    {
+        if (allowGpuPow2)
+        {
+            GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, divisorData.Modulus, out ulong remainder, divisorData);
+            if (status == GpuPow2ModStatus.Success)
+            {
+                return remainder == 1UL;
+            }
+
+            PrimeOrderGpuHeuristics.ReportPow2Failure(
+                "divisor cycle reduction",
+                divisorData.Modulus,
+                status);
+        }
+
+        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+    }
+
     private static ulong PollardRho64(ulong n)
     {
         if ((n & 1UL) == 0UL)
@@ -924,7 +946,7 @@ public class MersenneDivisorCycles
         return CalculateCycleLengthFallback(divisor);
     }
 
-    internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength)
+    internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool allowGpuPow2 = false)
     {
         if ((divisor & (divisor - 1UL)) == 0UL)
         {
@@ -945,7 +967,8 @@ public class MersenneDivisorCycles
                     previousOrder: null,
                     divisorData,
                     PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault,
-                    PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu);
+                    PrimeOrderCalculator.PrimeOrderHeuristicDevice.Cpu,
+                    allowGpuPow2);
             if (orderResult.Order != 0UL)
             {
                 cycleLength = orderResult.Order;

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -31,7 +31,8 @@ internal static partial class PrimeOrderCalculator
         PrimeOrderSearchConfig config,
         PrimeOrderHeuristicDevice device)
     {
-        using var scope = UsePow2Mode(device);
+        bool enableGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu;
+        using var scope = UsePow2Mode(device, enableGpuPow2);
         MontgomeryDivisorData divisorData;
         if (prime <= ulong.MaxValue)
         {
@@ -51,7 +52,7 @@ internal static partial class PrimeOrderCalculator
 
             ulong prime64 = (ulong)prime;
             divisorData = MontgomeryDivisorData.FromModulus(prime64);
-            PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device);
+            PrimeOrderResult result = Calculate(prime64, previous, divisorData, config, device, enableGpuPow2);
             return new PrimeOrderResultWide(result.Status, result.Order);
         }
         else
@@ -1030,6 +1031,11 @@ internal static partial class PrimeOrderCalculator
             {
                 return wideRemainder;
             }
+
+            PrimeOrderGpuHeuristics.ReportPow2Failure(
+                "wide prime order pow2 check",
+                modulus,
+                wideStatus);
         }
 
         return exponent.Pow2MontgomeryModWindowed(modulus);

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -90,13 +90,13 @@ internal static partial class PrimeOrderCalculator
         private readonly bool _previousAllow;
         private readonly PrimeOrderHeuristicDevice _previousDevice;
 
-        public Pow2ModeScope(PrimeOrderHeuristicDevice device)
+        public Pow2ModeScope(PrimeOrderHeuristicDevice device, bool allowGpuPow2)
         {
             _previousInitialized = s_pow2ModeInitialized;
             _previousAllow = s_allowGpuPow2;
             _previousDevice = s_deviceMode;
             s_pow2ModeInitialized = true;
-            s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu;
+            s_allowGpuPow2 = allowGpuPow2;
             s_deviceMode = device;
         }
 
@@ -124,7 +124,7 @@ internal static partial class PrimeOrderCalculator
         }
     }
 
-    private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device);
+    private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device, bool allowGpuPow2) => new(device, allowGpuPow2);
 
     private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled);
 
@@ -154,9 +154,11 @@ internal static partial class PrimeOrderCalculator
         ulong? previousOrder,
         in MontgomeryDivisorData divisorData,
         PrimeOrderSearchConfig config,
-        PrimeOrderHeuristicDevice device)
+        PrimeOrderHeuristicDevice device,
+        bool allowGpuPow2 = false)
     {
-        using var scope = UsePow2Mode(device);
+        bool enableGpuPow2 = allowGpuPow2 || device == PrimeOrderHeuristicDevice.Gpu;
+        using var scope = UsePow2Mode(device, enableGpuPow2);
         return CalculateInternal(prime, previousOrder, divisorData, config);
     }
 
@@ -443,7 +445,6 @@ internal static partial class PrimeOrderCalculator
         int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks;
         int powUsed = 0;
         int candidateCount = candidates.Count;
-        bool allowGpuBatch = true;
         Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
 
         DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
@@ -467,9 +468,11 @@ internal static partial class PrimeOrderCalculator
             bool gpuSuccess = false;
             bool gpuStackRemainders = false;
             GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
+            bool attemptedGpu = false;
 
-            if (allowGpuBatch && IsGpuPow2Allowed)
+            if (IsGpuPow2Allowed)
             {
+                attemptedGpu = true;
                 if (batchSize <= StackGpuBatchSize)
                 {
                     Span<ulong> localRemainders = stackGpuRemainders.Slice(0, batchSize);
@@ -496,11 +499,14 @@ internal static partial class PrimeOrderCalculator
                         gpuPool = null;
                     }
                 }
+            }
 
-                if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable))
+            if (attemptedGpu && !gpuSuccess)
-                {
+            {
-                    allowGpuBatch = false;
+                PrimeOrderGpuHeuristics.ReportPow2Failure(
-                }
+                    "prime order candidate scan",
+                    prime,
+                    status);
             }
 
             for (int i = 0; i < batchSize && powUsed < powBudget; i++)
@@ -780,6 +786,11 @@ internal static partial class PrimeOrderCalculator
             {
                 return remainder == 1UL;
             }
+
+            PrimeOrderGpuHeuristics.ReportPow2Failure(
+                "prime order pow2 check",
+                prime,
+                status);
         }
 
         return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;