Add shift-first multiply helper and document benchmarks

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -0,0 +1,67 @@
+using System;
+using System.Collections.Concurrent;
+
+namespace PerfectNumbers.Core.Gpu;
+
+internal enum GpuPow2ModStatus
+{
+    Success,
+    Overflow,
+    Unavailable,
+}
+
+internal static class PrimeOrderGpuHeuristics
+{
+    private static readonly ConcurrentDictionary<ulong, byte> OverflowedPrimes = new();
+
+    internal static void MarkOverflow(ulong prime) => OverflowedPrimes[prime] = 0;
+
+    internal static void ClearOverflow(ulong prime) => OverflowedPrimes.TryRemove(prime, out _);
+
+    internal static void ClearAllOverflowForTesting() => OverflowedPrimes.Clear();
+
+    public static GpuPow2ModStatus TryPow2Mod(ulong exponent, ulong prime, out ulong remainder)
+    {
+        remainder = 0UL;
+        if (prime <= 1UL)
+        {
+            return GpuPow2ModStatus.Unavailable;
+        }
+
+        if (OverflowedPrimes.ContainsKey(prime))
+        {
+            return GpuPow2ModStatus.Overflow;
+        }
+
+        // TODO: Implement the GPU heuristic order calculator so this path mirrors PrimeOrderCalculator before enabling device execution.
+        return GpuPow2ModStatus.Unavailable;
+    }
+
+    public static GpuPow2ModStatus TryPow2ModBatch(ReadOnlySpan<ulong> exponents, ulong prime, Span<ulong> remainders)
+    {
+        if (exponents.Length == 0)
+        {
+            return GpuPow2ModStatus.Success;
+        }
+
+        if (remainders.Length < exponents.Length)
+        {
+            throw new ArgumentException("Remainder span is shorter than the exponent span.", nameof(remainders));
+        }
+
+        remainders.Slice(0, exponents.Length).Clear();
+
+        if (prime <= 1UL)
+        {
+            return GpuPow2ModStatus.Unavailable;
+        }
+
+        if (OverflowedPrimes.ContainsKey(prime))
+        {
+            return GpuPow2ModStatus.Overflow;
+        }
+
+        // TODO: Implement the GPU heuristic order calculator so this path mirrors PrimeOrderCalculator before enabling device execution.
+        return GpuPow2ModStatus.Unavailable;
+    }
+}

PerfectNumbers.Core/HeuristicArithmetic.cs

@@ -0,0 +1,48 @@
+using System;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+internal static class HeuristicArithmetic
+{
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong MultiplyShiftRight(ulong value, ulong multiplier, int shift)
+    {
+        if (shift < 0)
+        {
+            throw new ArgumentOutOfRangeException(nameof(shift), "Shift must be non-negative.");
+        }
+
+        UInt128 product = (UInt128)value * multiplier;
+        return (ulong)(product >> shift);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift)
+    {
+        if (shift < 0)
+        {
+            throw new ArgumentOutOfRangeException(nameof(shift), "Shift must be non-negative.");
+        }
+
+        if (shift == 0)
+        {
+            return MultiplyShiftRight(value, multiplier, 0);
+        }
+
+        if (shift >= 64)
+        {
+            return 0UL;
+        }
+
+        ulong high = value >> shift;
+        ulong mask = (1UL << shift) - 1UL;
+        ulong low = value & mask;
+
+        UInt128 highContribution = (UInt128)high * multiplier;
+        UInt128 lowContribution = (UInt128)low * multiplier;
+
+        UInt128 combined = highContribution + (lowContribution >> shift);
+        return (ulong)combined;
+    }
+}

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -189,6 +189,16 @@ public class MersenneDivisorCycles
             return true;
         }
 
+        PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate(
+            divisor,
+            previousOrder: null,
+            PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault);
+        if (orderResult.Order != 0UL)
+        {
+            cycleLength = orderResult.Order;
+            return true;
+        }
+
         ulong phi = divisor - 1UL;
         if (phi == 0UL)
         {

PerfectNumbers.Core/PerfectNumbers.Core.csproj

@@ -1,4 +1,4 @@
-<Project Sdk="Microsoft.NET.Sdk">
+<Project Sdk="Microsoft.NET.Sdk">
 
   <PropertyGroup>
     <TargetFramework>net8.0</TargetFramework>
@@ -17,6 +17,7 @@
   <ItemGroup>
     <InternalsVisibleTo Include="EvenPerfectBitScanner.Tests" />
     <InternalsVisibleTo Include="PerfectNumbers.Core.Tests" />
+    <InternalsVisibleTo Include="HeuristicMultiplyShiftTiming" />
     <InternalsVisibleTo Include="EvenPerfectBitScanner.Benchmarks" />
   </ItemGroup>
 

PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs

@@ -3,6 +3,7 @@ using System.Collections.Generic;
 using System.Diagnostics;
 using System.Numerics;
 using System.Threading;
+using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core;
 
@@ -954,7 +955,17 @@ internal static partial class PrimeOrderCalculator
         {
             return UInt128.Zero;
         }
-
+        if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue)
+        {
+            ulong prime64 = (ulong)modulus;
+            ulong exponent64 = (ulong)exponent;
+            // TODO: Extend PrimeOrderGpuHeuristics to handle 128-bit moduli so this fast path can remain on the device.
+            GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent64, prime64, out ulong remainder);
+            if (status == GpuPow2ModStatus.Success)
+            {
+                return remainder;
+            }
+        }
         BigInteger result = BigInteger.ModPow(2, (BigInteger)exponent, (BigInteger)modulus);
         return (UInt128)result;
     }

PerfectNumbers.Core/PrimeOrderCalculator.cs

@@ -4,7 +4,9 @@ using System.Collections.Generic;
 using System.Diagnostics;
 using System.Numerics;
 using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
 using System.Security.Cryptography;
+using PerfectNumbers.Core.Gpu;
 using System.Threading;
 
 namespace PerfectNumbers.Core;
@@ -134,7 +136,7 @@ internal static partial class PrimeOrderCalculator
         {
             ulong factor = factorSpan[i].Value;
             ulong reduced = phi / factor;
-            if (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+            if (Pow2EqualsOne(reduced, prime, divisorData))
             {
                 return false;
             }
@@ -149,7 +151,7 @@ internal static partial class PrimeOrderCalculator
         if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL)
         {
             ulong half = phi >> 1;
-            if (half.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+            if (Pow2EqualsOne(half, prime, divisorData))
             {
                 order = half;
             }
@@ -191,7 +193,7 @@ internal static partial class PrimeOrderCalculator
                     }
 
                     ulong reduced = order / primeFactor;
-                    if (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+                    if (Pow2EqualsOne(reduced, prime, divisorData))
                     {
                         order = reduced;
                         continue;
@@ -223,7 +225,7 @@ internal static partial class PrimeOrderCalculator
 #if DEBUG
         Console.WriteLine("Verifying a^order ≡ 1 (mod p)");
 #endif
-        if (order.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) != 1UL)
+        if (!Pow2EqualsOne(order, prime, divisorData))
         {
             return false;
         }
@@ -283,7 +285,7 @@ internal static partial class PrimeOrderCalculator
                 }
 
                 reduced /= primeFactor;
-                if (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+                if (Pow2EqualsOne(reduced, prime, divisorData))
                 {
                     return false;
                 }
@@ -355,32 +357,108 @@ 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);
 
 #if DEBUG
         Console.WriteLine($"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
 #endif
-        for (int i = 0; i < candidateCount; i++)
+        int index = 0;
+        const int MaxGpuBatchSize = 256;
+        const int StackGpuBatchSize = 32;
+        Span<ulong> stackGpuRemainders = stackalloc ulong[StackGpuBatchSize];
+        while (index < candidateCount && powUsed < powBudget)
         {
-            if (powUsed >= powBudget)
+            int remaining = candidateCount - index;
+            int budgetRemaining = powBudget - powUsed;
+            int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize));
+            if (batchSize <= 0)
             {
                 break;
             }
 
-            ulong candidate = candidates[i];
-            powUsed++;
+            ReadOnlySpan<ulong> batch = candidateSpan.Slice(index, batchSize);
+            ulong[]? gpuPool = null;
+            Span<ulong> pooledGpuRemainders = default;
+            bool gpuSuccess = false;
+            bool gpuStackRemainders = false;
+            GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable;
 
-            if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) != 1UL)
+            if (allowGpuBatch)
             {
-                continue;
+                if (batchSize <= StackGpuBatchSize)
+                {
+                    Span<ulong> localRemainders = stackGpuRemainders.Slice(0, batchSize);
+                    status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders);
+                    if (status == GpuPow2ModStatus.Success)
+                    {
+                        gpuSuccess = true;
+                        gpuStackRemainders = true;
+                    }
+                }
+                else
+                {
+                    gpuPool = ArrayPool<ulong>.Shared.Rent(batchSize);
+                    Span<ulong> pooledRemainders = gpuPool.AsSpan(0, batchSize);
+                    status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders);
+                    if (status == GpuPow2ModStatus.Success)
+                    {
+                        pooledGpuRemainders = pooledRemainders;
+                        gpuSuccess = true;
+                    }
+                    else
+                    {
+                        ArrayPool<ulong>.Shared.Return(gpuPool, clearArray: false);
+                        gpuPool = null;
+                    }
+                }
+
+                if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable))
+                {
+                    allowGpuBatch = false;
+                }
             }
 
-            if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget))
+            for (int i = 0; i < batchSize && powUsed < powBudget; i++)
             {
-                continue;
+                ulong candidate = batch[i];
+                powUsed++;
+
+                bool equalsOne;
+                if (gpuSuccess)
+                {
+                    ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i];
+                    equalsOne = remainderValue == 1UL;
+                }
+                else
+                {
+                    equalsOne = Pow2EqualsOne(candidate, prime, divisorData);
+                }
+                if (!equalsOne)
+                {
+                    continue;
+                }
+
+                if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget))
+                {
+                    continue;
+                }
+
+                if (gpuPool is not null)
+                {
+                    ArrayPool<ulong>.Shared.Return(gpuPool, clearArray: false);
+                }
+
+                result = candidate;
+                return true;
             }
 
-            result = candidate;
-            return true;
+            if (gpuPool is not null)
+            {
+                ArrayPool<ulong>.Shared.Return(gpuPool, clearArray: false);
+            }
+
+            index += batchSize;
         }
 
 #if DEBUG
@@ -601,7 +679,7 @@ internal static partial class PrimeOrderCalculator
                 }
 
                 powUsed++;
-                if (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+                if (Pow2EqualsOne(reduced, prime, divisorData))
                 {
                     return false;
                 }
@@ -611,6 +689,18 @@ internal static partial class PrimeOrderCalculator
         return true;
     }
 
+    private static bool Pow2EqualsOne(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData)
+    {
+        // TODO: Implement the GPU heuristic path so TryPow2Mod can verify candidates on the device instead of immediately falling back to the CPU Montgomery reducer.
+        GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder);
+        if (status == GpuPow2ModStatus.Success)
+        {
+            return remainder == 1UL;
+        }
+
+        return exponent.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL;
+    }
+
     private static PartialFactorResult PartialFactor(ulong value, PrimeOrderSearchConfig config)
     {
         if (value <= 1UL)