Switch GPU pow2mod helpers to plain windowed arithmetic

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -545,10 +545,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
     private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
     {
-        if (prime <= 1UL)
-        {
-            return 0UL;
-        }
+        // The by-divisor search emits candidate exponents starting at the first odd prime, so this guard
+        // would never trigger on the scanner hot path.
+        // if (prime <= 1UL)
+        // {
+        //     return 0UL;
+        // }
 
         if (prime - 1UL >= 64UL)
         {
@@ -562,14 +564,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
     {
         MontgomeryDivisorData divisor = divisors[index];
         ulong modulus = divisor.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            hits[index] = 0;
-            return;
-        }
+        // The divisor staging path filters out non-prime and even moduli before they reach the device, so this guard
+        // would never fire for EvenPerfectBitScanner workloads.
+        // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        // {
+        //     hits[index] = 0;
+        //     return;
+        // }
 
         ulong exponent = exponents[index];
-        hits[index] = exponent.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = exponent.Pow2ModWindowedGpu(modulus) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -780,15 +784,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
     private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
     {
         ulong modulus = divisor.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            results[index] = 0UL;
-            return;
-        }
+        // Kernel inputs are constrained to odd prime moduli, so this defensive branch would never execute on the configured
+        // scanning path.
+        // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        // {
+        //     results[index] = 0UL;
+        //     return;
+        // }
 
         ulong exponent = exponents[index];
 
-        results[index] = exponent.Pow2ModBinaryGpu(modulus);
+        results[index] = exponent.Pow2ModWindowedGpu(modulus);
     }
 
     public ulong DivisorLimit

PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs

@@ -110,12 +110,14 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
 
         ulong exponent = exponents[idx];
         ulong modulus = divisorData.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            results[idx] = 0UL;
-            return;
-        }
+        // Residue batches only include valid odd moduli, so the defensive check stays commented to avoid re-validating
+        // the invariant on every iteration.
+        // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        // {
+        //     results[idx] = 0UL;
+        //     return;
+        // }
 
-        results[idx] = exponent.Pow2ModBinaryGpu(modulus);
+        results[idx] = exponent.Pow2ModWindowedGpu(modulus);
     }
 }

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -805,13 +805,15 @@ internal static partial class PrimeOrderGpuHeuristics
         for (int i = 0; i < factorCount; i++)
         {
             ulong factor = factors[i];
-            if (factor <= 1UL)
-            {
-                continue;
-            }
+            // The factorization routine emits strictly greater-than-one factors for phi on this path, so the
+            // defensive continue would never execute.
+            // if (factor <= 1UL)
+            // {
+            //     continue;
+            // }
 
             ulong reduced = phi / factor;
-            if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL)
+            if (reduced.Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
             {
                 return false;
             }
@@ -827,7 +829,7 @@ internal static partial class PrimeOrderGpuHeuristics
         if (residue == 1UL || residue == 7UL)
         {
             ulong half = phi >> 1;
-            if (Pow2EqualsOneKernel(half, divisor))
+            if ((half).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
             {
                 order = half;
             }
@@ -861,7 +863,7 @@ internal static partial class PrimeOrderGpuHeuristics
                 }
 
                 ulong reduced = order / primeFactor;
-                if (Pow2EqualsOneKernel(reduced, divisor))
+                if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
                 {
                     order = reduced;
                     continue;
@@ -892,7 +894,7 @@ internal static partial class PrimeOrderGpuHeuristics
             return false;
         }
 
-        if (!Pow2EqualsOneKernel(order, divisor))
+        if ((order).Pow2ModWindowedGpu(divisor.Modulus) != 1UL)
         {
             return false;
         }
@@ -936,7 +938,7 @@ internal static partial class PrimeOrderGpuHeuristics
                 }
 
                 reduced /= primeFactor;
-                if (Pow2EqualsOneKernel(reduced, divisor))
+                if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
                 {
                     return false;
                 }
@@ -1027,7 +1029,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             powUsed++;
-            if (!Pow2EqualsOneKernel(candidate, divisor))
+            if ((candidate).Pow2ModWindowedGpu(divisor.Modulus) != 1UL)
             {
                 continue;
             }
@@ -1327,7 +1329,7 @@ internal static partial class PrimeOrderGpuHeuristics
                 }
 
                 powUsed++;
-                if (Pow2EqualsOneKernel(reduced, divisor))
+                if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
                 {
                     return false;
                 }
@@ -1593,11 +1595,6 @@ internal static partial class PrimeOrderGpuHeuristics
         }
     }
 
-    private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
-    {
-        return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL;
-    }
-
     private static ulong CalculateByDoublingKernel(ulong prime)
     {
         ulong order = 1UL;
@@ -1698,9 +1695,11 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void ComputePow2ModCpu(ReadOnlySpan<ulong> exponents, ulong prime, in MontgomeryDivisorData divisorData, Span<ulong> results)
     {
         int length = exponents.Length;
+        ulong modulus = divisorData.Modulus;
+        _ = prime; // modulus and divisorData both describe the same q; retain the parameter to preserve the public signature.
         for (int i = 0; i < length; i++)
         {
-            results[i] = Pow2ModCpu(exponents[i], prime, divisorData);
+            results[i] = exponents[i].Pow2ModWindowedCpu(modulus);
         }
     }
 
@@ -1788,16 +1787,6 @@ internal static partial class PrimeOrderGpuHeuristics
         }
     }
 
-    private static ulong Pow2ModCpu(ulong exponent, ulong modulus, in MontgomeryDivisorData divisorData)
-    {
-        if (modulus <= 1UL)
-        {
-            return 0UL;
-        }
-
-        return exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: false);
-    }
-
     private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
     {
         return Pow2ModKernelCache.GetOrAdd(accelerator, static accel =>
@@ -1811,14 +1800,16 @@ internal static partial class PrimeOrderGpuHeuristics
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong modulus = divisor.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            remainders[index] = 0UL;
-            return;
-        }
+        // The by-divisor heuristics call this kernel only with odd prime moduli, so the defensive branch stays commented to
+        // document the invariant without spending a branch in the kernel.
+        // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        // {
+        //     remainders[index] = 0UL;
+        //     return;
+        // }
 
         ulong exponent = exponents[index];
-        remainders[index] = exponent.Pow2ModBinaryGpu(modulus);
+        remainders[index] = exponent.Pow2ModWindowedGpu(modulus);
     }
 
     private static Action<AcceleratorStream, Index1D, ArrayView1D<GpuUInt128, Stride1D.Dense>, GpuUInt128, ArrayView1D<GpuUInt128, Stride1D.Dense>> GetPow2ModWideKernel(Accelerator accelerator)

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -81,6 +81,8 @@ public class MersenneDivisorCycles
 
         if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
         {
+            // Even though production scans never revisit divisors below this bound, the shared cache still serves
+            // precomputation steps and tests, so keep the fast lookup enabled.
             var small = Shared._smallCycles;
             if (small is not null)
             {
@@ -98,6 +100,7 @@ public class MersenneDivisorCycles
         }
 
         Dictionary<ulong, FactorCacheEntry>? cache = factorCache;
+        // Each by-divisor scan session starts with a null cache and fills it lazily, so keep the allocation guard in place.
         if (cache is null)
         {
             cache = new Dictionary<ulong, FactorCacheEntry>(8);
@@ -121,6 +124,8 @@ public class MersenneDivisorCycles
 
         if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
         {
+            // The by-divisor path only handles divisors far above this cutoff, but Pollard Rho factoring and unit tests
+            // still reuse the cached small-cycle table, so keep the fast path enabled.
             MersenneDivisorCycles shared = Shared;
             ulong[]? small = shared._smallCycles;
             if (small is not null)
@@ -138,12 +143,16 @@ public class MersenneDivisorCycles
 
     private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
     {
+        // The divisor enumeration pipeline has already filtered invalid inputs, but keep the guard so external callers and
+        // tests continue to receive a defensive answer.
         if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
         {
             return false;
         }
 
         UInt128 step = (UInt128)exponent << 1;
+        // Exponents drawn from prime candidates are far below the overflow boundary, yet the broader factoring utilities rely
+        // on this bail-out when validating arbitrary inputs.
         if (step == UInt128.Zero)
         {
             return false;

PerfectNumbers.Core/ULongExtensions.cs

@@ -12,6 +12,8 @@ public static class ULongExtensions
     private const int Pow2WindowSize = 8;
     private const ulong Pow2WindowFallbackThreshold = 32UL;
     private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1);
+    private const int Pow2WindowSizeGpu = 5;
+    private const int Pow2WindowOddPowerCountGpu = 1 << (Pow2WindowSizeGpu - 1);
 
     public static ulong CalculateOrder(this ulong q)
     {
@@ -310,25 +312,26 @@ public static class ULongExtensions
         return state.MulModWithNativeModulo(b, modulus);
     }
 
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2ModBinaryCpu(this ulong exponent, ulong modulus)
-    {
-        return Pow2ModWindowedCpu(exponent, modulus);
-    }
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong Pow2ModWindowedCpu(this ulong exponent, ulong modulus)
     {
-        if (modulus <= 1UL)
-        {
-            return 0UL;
-        }
+        // The by-divisor pipeline only schedules odd prime moduli, so the defensive check below
+        // would never trigger on that path. Leaving it commented documents the invariant without
+        // forcing a branch in the hot loop.
+        // if (modulus <= 1UL)
+        // {
+        //     return 0UL;
+        // }
 
-        if (exponent == 0UL)
-        {
-            return 1UL % modulus;
-        }
+        // Exponents flow from the strictly positive candidate list on the by-divisor scanner, so
+        // this guard would be redundant there.
+        // if (exponent == 0UL)
+        // {
+        //     return 1UL % modulus;
+        // }
 
+        // Even though the by-divisor pipeline only feeds very large exponents, the legacy callers and unit tests
+        // still hit this small-exponent fast path, so keep the fallback in place.
         if (exponent <= Pow2WindowFallbackThreshold)
         {
             return Pow2ModBinaryFallback(exponent, modulus);
@@ -341,11 +344,12 @@ public static class ULongExtensions
         Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
             ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount]
             : new ulong[oddPowerCount];
-        // Reuse oddPowerStorage to expose only the populated slice and avoid another local span.
         oddPowerStorage = oddPowerStorage[..oddPowerCount];
         InitializePlainOddPowers(modulus, oddPowerStorage);
 
-        ulong result = 1UL % modulus;
+        // With modulus guaranteed to exceed one on the scanning path, the modulo operation would
+        // always yield one, so keep the literal to highlight the intended base value.
+        ulong result = 1UL;
         int index = bitLength - 1;
 
         while (index >= 0)
@@ -387,6 +391,8 @@ public static class ULongExtensions
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static ulong Pow2ModBinaryFallback(ulong exponent, ulong modulus)
     {
+        // The scanning path keeps modulus strictly greater than one, but other utility callers still
+        // rely on this fallback for tiny exponents, so leave the modulo in place for safety.
         ulong remainingExponent = exponent;
         ulong result = 1UL % modulus;
         ulong baseValue = 2UL % modulus;
@@ -411,13 +417,77 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus)
+    public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus)
     {
-        if (modulus <= 1UL)
+        // The GPU pipeline feeds only odd prime moduli, so the removed guard would never fire for
+        // EvenPerfectBitScanner workloads.
+        // if (modulus <= 1UL)
+        // {
+        //     return 0UL;
+        // }
+
+        // Exponents fed into the GPU kernels are positive on the scanning path, yet the zero case remains for
+        // completeness when unit tests probe boundary behavior.
+        if (exponent == 0UL)
         {
-            return 0UL;
+            return 1UL;
         }
 
+        // The GPU kernels mostly process large exponents as well, but the binary fallback is still useful when
+        // unit tests probe edge cases around the threshold.
+        if (exponent <= Pow2WindowFallbackThreshold)
+        {
+            return Pow2ModBinaryGpuFallback(exponent, modulus);
+        }
+
+        int bitLength = GetPortableBitLength(exponent);
+        int windowSize = GetGpuWindowSize(bitLength);
+        int oddPowerCount = 1 << (windowSize - 1);
+
+        PlainOddPowerTableGpu oddPowers = new PlainOddPowerTableGpu(modulus, oddPowerCount);
+
+        ulong result = 1UL;
+        int index = bitLength - 1;
+
+        while (index >= 0)
+        {
+            if (((exponent >> index) & 1UL) == 0UL)
+            {
+                result = result.MulMod64GpuCompatible(result, modulus);
+                index--;
+                continue;
+            }
+
+            int windowStart = index - windowSize + 1;
+            if (windowStart < 0)
+            {
+                windowStart = 0;
+            }
+
+            while (((exponent >> windowStart) & 1UL) == 0UL)
+            {
+                windowStart++;
+            }
+
+            int windowLength = index - windowStart + 1;
+            for (int square = 0; square < windowLength; square++)
+            {
+                result = result.MulMod64GpuCompatible(result, modulus);
+            }
+
+            ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL);
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+            result = result.MulMod64GpuCompatible(oddPowers.Get(tableIndex), modulus);
+
+            index = windowStart - 1;
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Pow2ModBinaryGpuFallback(ulong exponent, ulong modulus)
+    {
         ulong remainingExponent = exponent;
         ulong result = 1UL % modulus;
         ulong baseValue = 2UL % modulus;
@@ -678,10 +748,19 @@ public static class ULongExtensions
         return Pow2WindowSize;
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static int GetGpuWindowSize(int bitLength)
+    {
+        int window = GetWindowSize(bitLength);
+        return window > Pow2WindowSizeGpu ? Pow2WindowSizeGpu : window;
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static void InitializePlainOddPowers(ulong modulus, Span<ulong> destination)
     {
         int oddPowerCount = destination.Length;
+        // The by-divisor kernels always request at least one odd power, but other helpers reuse this initializer
+        // for smaller workloads, so the early exit stays in place.
         if (oddPowerCount == 0)
         {
             return;
@@ -703,6 +782,150 @@ public static class ULongExtensions
         }
     }
 
+    private struct PlainOddPowerTableGpu
+    {
+        private readonly int _count;
+        private ulong _value0;
+        private ulong _value1;
+        private ulong _value2;
+        private ulong _value3;
+        private ulong _value4;
+        private ulong _value5;
+        private ulong _value6;
+        private ulong _value7;
+        private ulong _value8;
+        private ulong _value9;
+        private ulong _value10;
+        private ulong _value11;
+        private ulong _value12;
+        private ulong _value13;
+        private ulong _value14;
+        private ulong _value15;
+
+        public PlainOddPowerTableGpu(ulong modulus, int requestedCount)
+        {
+            _count = requestedCount > Pow2WindowOddPowerCountGpu ? Pow2WindowOddPowerCountGpu : requestedCount;
+            _value0 = 0UL;
+            _value1 = 0UL;
+            _value2 = 0UL;
+            _value3 = 0UL;
+            _value4 = 0UL;
+            _value5 = 0UL;
+            _value6 = 0UL;
+            _value7 = 0UL;
+            _value8 = 0UL;
+            _value9 = 0UL;
+            _value10 = 0UL;
+            _value11 = 0UL;
+            _value12 = 0UL;
+            _value13 = 0UL;
+            _value14 = 0UL;
+            _value15 = 0UL;
+
+            if (_count == 0)
+            {
+                return;
+            }
+
+            ulong baseValue = 2UL % modulus;
+            _value0 = baseValue;
+            if (_count == 1)
+            {
+                return;
+            }
+
+            ulong square = baseValue.MulMod64GpuCompatible(baseValue, modulus);
+            ulong current = baseValue;
+            for (int i = 1; i < _count; i++)
+            {
+                current = current.MulMod64GpuCompatible(square, modulus);
+                SetValue(i, current);
+            }
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public ulong Get(int index)
+        {
+            if ((uint)index >= (uint)_count)
+            {
+                return 0UL;
+            }
+
+            return index switch
+            {
+                0 => _value0,
+                1 => _value1,
+                2 => _value2,
+                3 => _value3,
+                4 => _value4,
+                5 => _value5,
+                6 => _value6,
+                7 => _value7,
+                8 => _value8,
+                9 => _value9,
+                10 => _value10,
+                11 => _value11,
+                12 => _value12,
+                13 => _value13,
+                14 => _value14,
+                _ => _value15,
+            };
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private void SetValue(int index, ulong value)
+        {
+            switch (index)
+            {
+                case 1:
+                    _value1 = value;
+                    break;
+                case 2:
+                    _value2 = value;
+                    break;
+                case 3:
+                    _value3 = value;
+                    break;
+                case 4:
+                    _value4 = value;
+                    break;
+                case 5:
+                    _value5 = value;
+                    break;
+                case 6:
+                    _value6 = value;
+                    break;
+                case 7:
+                    _value7 = value;
+                    break;
+                case 8:
+                    _value8 = value;
+                    break;
+                case 9:
+                    _value9 = value;
+                    break;
+                case 10:
+                    _value10 = value;
+                    break;
+                case 11:
+                    _value11 = value;
+                    break;
+                case 12:
+                    _value12 = value;
+                    break;
+                case 13:
+                    _value13 = value;
+                    break;
+                case 14:
+                    _value14 = value;
+                    break;
+                case 15:
+                    _value15 = value;
+                    break;
+            }
+        }
+    }
+
     private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> destination)
     {
         int oddPowerCount = destination.Length;