Windowed Pow2Mod implemented and integrated across repository.

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -1,5 +1,6 @@
 using System;
 using System.Collections.Concurrent;
+using System.Collections.Generic;
 using System.Numerics;
 using PerfectNumbers.Core;
 
@@ -173,6 +174,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             return false;
         }
 
+        Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
         ulong divisor = 3UL;
         while (divisor <= allowedMax)
         {
@@ -183,7 +185,19 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             // TODO: Hoist MontgomeryDivisorData acquisition into the divisor-cycle cache so we reuse the
             // staged ProcessEightBitWindows-ready metadata instead of reloading the slower Montgomery
             // structs for every divisor.
-            ulong divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
+            ulong divisorCycle;
+            if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+            {
+                divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
+            }
+            else
+            {
+                factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
+                if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, prime, factorCache, out divisorCycle) || divisorCycle == 0UL)
+                {
+                    divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
+                }
+            }
             if (divisorCycle != 0UL)
             {
                 byte hit = CheckDivisor(prime, divisorCycle, divisorData);
@@ -220,12 +234,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
     private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData)
     {
-        ulong modulus = divisorData.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL || divisorCycle == 0UL)
-        {
-            return 0;
-        }
-
         ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData);
         return residue == 1UL ? (byte)1 : (byte)0;
     }
@@ -308,12 +316,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 return;
             }
 
-            if (divisor <= 1UL || (divisor & 1UL) == 0UL)
-            {
-                hits.Clear();
-                return;
-            }
-
             if (divisorCycle == 0UL)
             {
                 divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);

PerfectNumbers.Core/ExponentRemainderStepper.cs

@@ -47,7 +47,7 @@ internal struct ExponentRemainderStepper
             // TODO: Route these state resets through the ProcessEightBitWindows helper once the scalar
             // Pow2MontgomeryMod implementation adopts it so delta stepping inherits the benchmarked
             // 2× gains for large exponents.
-            _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
             _previousExponent = exponent;
             _hasState = true;
             return ReduceCurrent();
@@ -59,7 +59,7 @@ internal struct ExponentRemainderStepper
         // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
         // TODO: Replace this per-delta powmod with the upcoming windowed ladder so incremental
         // updates stop paying the single-bit cost highlighted in GpuPow2ModBenchmarks.
-        ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor);
+        ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
         _previousExponent = exponent;
         return ReduceCurrent();
@@ -77,7 +77,7 @@ internal struct ExponentRemainderStepper
         {
             // TODO: Switch this reload to the shared windowed pow2 helper once available so CPU
             // residue checks align with the optimized ProcessEightBitWindows timings.
-            _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
             _previousExponent = exponent;
             _hasState = true;
             return _currentMontgomery == _montgomeryOne;
@@ -88,7 +88,7 @@ internal struct ExponentRemainderStepper
         // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor.
         // TODO: Use the windowed delta pow2 helper here as well to avoid the single-bit ladder that
         // currently lags behind the benchmarked implementation for huge divisor cycles.
-        ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor);
+        ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
         _previousExponent = exponent;
         return _currentMontgomery == _montgomeryOne;
@@ -123,7 +123,7 @@ internal struct ExponentRemainderStepper
         {
             // TODO: Replace this fallback path with the upcoming ProcessEightBitWindows helper so
             // fresh Montgomery states also benefit from the faster pow2 ladder measured on CPUs.
-            _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor);
+            _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true);
             _previousExponent = exponent;
             _hasState = true;
             isUnity = _currentMontgomery == _montgomeryOne;

PerfectNumbers.Core/Gpu/GpuUInt128.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 
@@ -15,6 +16,10 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     private const ulong NativeModuloChunkMask = (1UL << NativeModuloChunkBits) - 1UL;
     private static readonly ulong[] NativeModuloBitMasks = CreateNativeModuloBitMasks();
 
+    private const int Pow2WindowSizeBits = 8;
+    private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1);
+    private const int Pow2WindowMaxExponent = (1 << Pow2WindowSizeBits) - 1;
+
     public static readonly GpuUInt128 Zero = new();
     public static readonly GpuUInt128 One = new(1UL);
 
@@ -256,122 +261,24 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
             return Zero;
         }
 
-        // TODO: Replace this single-bit ladder with the ProcessEightBitWindows strategy measured
-        // in GpuPow2ModBenchmarks once the GPU kernels can share the windowed helper. The windowed
-        // version held 21.5 µs for 33-bit moduli versus 51.0 µs here when scanning large divisors.
-        GpuUInt128 result = Zero;
-        ulong i = 0UL;
-        while (i < exponent)
+        GpuUInt128 pow = Pow2Mod(exponent, modulus);
+        if (pow.IsZero)
         {
-            // result = (result * 2) % modulus
-            result += result;
-            if (result.CompareTo(modulus) >= 0)
-            {
-                result -= modulus;
-            }
-
-            // result = (result + 1) % modulus
-            result += One;
-            if (result.CompareTo(modulus) >= 0)
-            {
-                result -= modulus;
-            }
-
-            i++;
+            pow = modulus;
         }
 
-        return result;
+        pow.Sub(One);
+        return pow;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static void Pow2Minus1ModBatch(GpuUInt128 modulus, ReadOnlySpan<ulong> exponents, Span<GpuUInt128> results)
     {
-        GpuUInt128 r0 = Zero, r1 = Zero, r2 = Zero, r3 = Zero, r4 = Zero, r5 = Zero, r6 = Zero, r7 = Zero;
-        ulong e0 = 0UL, e1 = 0UL, e2 = 0UL, e3 = 0UL, e4 = 0UL, e5 = 0UL, e6 = 0UL, e7 = 0UL;
-        bool a0 = false, a1 = false, a2 = false, a3 = false, a4 = false, a5 = false, a6 = false, a7 = false;
-
-        // TODO: Mirror the eight-bit window batching once the scalar helper above switches; the
-        // current bit-serial loop lags by ~3× on 8k–33-bit moduli compared to the windowed path.
-        int len = exponents.Length;
-        if (len > 0) { e0 = exponents[0]; a0 = true; }
-        if (len > 1) { e1 = exponents[1]; a1 = true; }
-        if (len > 2) { e2 = exponents[2]; a2 = true; }
-        if (len > 3) { e3 = exponents[3]; a3 = true; }
-        if (len > 4) { e4 = exponents[4]; a4 = true; }
-        if (len > 5) { e5 = exponents[5]; a5 = true; }
-        if (len > 6) { e6 = exponents[6]; a6 = true; }
-        if (len > 7) { e7 = exponents[7]; a7 = true; }
-
-        ulong max = 0UL;
-        if (a0 && e0 > max) max = e0;
-        if (a1 && e1 > max) max = e1;
-        if (a2 && e2 > max) max = e2;
-        if (a3 && e3 > max) max = e3;
-        if (a4 && e4 > max) max = e4;
-        if (a5 && e5 > max) max = e5;
-        if (a6 && e6 > max) max = e6;
-        if (a7 && e7 > max) max = e7;
-
-        for (ulong i = 0UL; i < max; i++)
+        int count = exponents.Length;
+        for (int i = 0; i < count; i++)
         {
-            if (a0)
-            {
-                r0 += r0; if (r0.CompareTo(modulus) >= 0) r0 -= modulus;
-                r0 += One; if (r0.CompareTo(modulus) >= 0) r0 -= modulus;
-                if (i + 1UL == e0) a0 = false;
-            }
-            if (a1)
-            {
-                r1 += r1; if (r1.CompareTo(modulus) >= 0) r1 -= modulus;
-                r1 += One; if (r1.CompareTo(modulus) >= 0) r1 -= modulus;
-                if (i + 1UL == e1) a1 = false;
-            }
-            if (a2)
-            {
-                r2 += r2; if (r2.CompareTo(modulus) >= 0) r2 -= modulus;
-                r2 += One; if (r2.CompareTo(modulus) >= 0) r2 -= modulus;
-                if (i + 1UL == e2) a2 = false;
-            }
-            if (a3)
-            {
-                r3 += r3; if (r3.CompareTo(modulus) >= 0) r3 -= modulus;
-                r3 += One; if (r3.CompareTo(modulus) >= 0) r3 -= modulus;
-                if (i + 1UL == e3) a3 = false;
-            }
-            if (a4)
-            {
-                r4 += r4; if (r4.CompareTo(modulus) >= 0) r4 -= modulus;
-                r4 += One; if (r4.CompareTo(modulus) >= 0) r4 -= modulus;
-                if (i + 1UL == e4) a4 = false;
-            }
-            if (a5)
-            {
-                r5 += r5; if (r5.CompareTo(modulus) >= 0) r5 -= modulus;
-                r5 += One; if (r5.CompareTo(modulus) >= 0) r5 -= modulus;
-                if (i + 1UL == e5) a5 = false;
-            }
-            if (a6)
-            {
-                r6 += r6; if (r6.CompareTo(modulus) >= 0) r6 -= modulus;
-                r6 += One; if (r6.CompareTo(modulus) >= 0) r6 -= modulus;
-                if (i + 1UL == e6) a6 = false;
-            }
-            if (a7)
-            {
-                r7 += r7; if (r7.CompareTo(modulus) >= 0) r7 -= modulus;
-                r7 += One; if (r7.CompareTo(modulus) >= 0) r7 -= modulus;
-                if (i + 1UL == e7) a7 = false;
-            }
+            results[i] = Pow2Minus1Mod(exponents[i], modulus);
         }
-
-        if (len > 0) results[0] = r0;
-        if (len > 1) results[1] = r1;
-        if (len > 2) results[2] = r2;
-        if (len > 3) results[3] = r3;
-        if (len > 4) results[4] = r4;
-        if (len > 5) results[5] = r5;
-        if (len > 6) results[6] = r6;
-        if (len > 7) results[7] = r7;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -483,7 +390,6 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         }
     }
 
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static GpuUInt128 Pow2Mod(ulong exponent, GpuUInt128 modulus)
     {
@@ -492,36 +398,71 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
             return new();
         }
 
-        const int WindowSizeBits = 8;
-        const ulong WindowMask = (1UL << WindowSizeBits) - 1UL;
-
-        GpuUInt128 result = new(1UL);
-        GpuUInt128 baseVal = new(2UL);
-
-        ulong e = exponent;
-        while (e != 0UL)
+        if (exponent == 0UL)
         {
-            int remainingBits = BitOperations.Log2(e) + 1;
-            int bitsToProcess = remainingBits >= WindowSizeBits ? WindowSizeBits : remainingBits;
-            ulong chunk = e & (WindowMask >> (WindowSizeBits - bitsToProcess));
+            return new GpuUInt128(1UL);
+        }
 
-            for (int bitIndex = 0; bitIndex < bitsToProcess; bitIndex++)
+        Span<GpuUInt128> oddPowers = stackalloc GpuUInt128[Pow2WindowOddPowerCount];
+        GpuUInt128 power = new GpuUInt128(1UL);
+        int oddIndex = 0;
+
+        for (int bit = 1; bit <= Pow2WindowMaxExponent; bit++)
+        {
+            power += power;
+            if (power.CompareTo(modulus) >= 0)
             {
-                if ((chunk & 1UL) != 0UL)
-                {
-                    result.MulMod(baseVal, modulus);
-                }
-
-                baseVal.MulMod(baseVal, modulus);
-                chunk >>= 1;
+                power.Sub(modulus);
             }
 
-            e >>= bitsToProcess;
+            if ((bit & 1) != 0)
+            {
+                oddPowers[oddIndex++] = power;
+            }
+        }
+
+        GpuUInt128 result = new GpuUInt128(1UL);
+        int bitLength = BitOperations.Log2(exponent) + 1;
+        int index = bitLength - 1;
+
+        while (index >= 0)
+        {
+            if (((exponent >> index) & 1UL) == 0UL)
+            {
+                result.SquareMod(modulus);
+                index--;
+                continue;
+            }
+
+            int windowStart = index - Pow2WindowSizeBits + 1;
+            if (windowStart < 0)
+            {
+                windowStart = 0;
+            }
+
+            while (((exponent >> windowStart) & 1UL) == 0UL)
+            {
+                windowStart++;
+            }
+
+            int windowBitCount = index - windowStart + 1;
+            for (int square = 0; square < windowBitCount; square++)
+            {
+                result.SquareMod(modulus);
+            }
+
+            ulong mask = (1UL << windowBitCount) - 1UL;
+            ulong windowValue = (exponent >> windowStart) & mask;
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+
+            GpuUInt128 factor = oddPowers[tableIndex];
+            result.MulMod(factor, modulus);
+
+            index = windowStart - 1;
         }
 
         return result;
     }
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static int TrailingZeroCount(GpuUInt128 value)
     {
@@ -1258,7 +1199,6 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         p3 = hHH + c2;
     }
 
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static ulong MulMod64(ulong a, ulong b, ulong modulus)
     {

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -1,6 +1,7 @@
 using System;
 using System.Buffers;
 using System.Collections.Concurrent;
+using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
 using System.Numerics;
@@ -371,7 +372,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulong exponent = exponents[index];
-        hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; // TODO: Replace Pow2MontgomeryMod with the ProcessEightBitWindows helper once available so GPU by-divisor scans reuse the benchmarked fast pow2 ladder.
+        hits[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -384,6 +385,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         private MemoryBuffer1D<ulong, Stride1D.Dense> _resultsBuffer;
         private ulong[] _hostBuffer;
         private int[]? _positionBuffer;
+        private Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? _factorCache;
         private int _capacity;
         private bool _disposed;
 
@@ -403,6 +405,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         internal void Reset()
         {
             _disposed = false;
+            _factorCache?.Clear();
         }
 
         public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
@@ -428,11 +431,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
             if (divisorCycle == 0UL)
             {
-                divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
-                if (divisorCycle == 0UL)
+                Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry> cache = _factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(16);
+                if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, primes[0], cache, out divisorCycle) || divisorCycle == 0UL)
                 {
-                    hits.Clear();
-                    return;
+                    divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
+                    if (divisorCycle == 0UL)
+                    {
+                        hits.Clear();
+                        return;
+                    }
                 }
             }
 
@@ -569,6 +576,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             }
 
             _disposed = true;
+            _factorCache?.Clear();
             _owner._sessionPool.Add(this);
         }
     }
@@ -584,7 +592,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         ulong exponent = exponents[index];
 
-        results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); // TODO: Swap this Montgomery-domain pow2 with the windowed helper slated to replace the slower ladder in Pow2Minus1Mod benchmarks.
+        results[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: true);
     }
 
     public ulong DivisorLimit

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -1,9 +1,14 @@
 using System.Buffers;
+using System.Buffers.Binary;
+using System.Collections.Generic;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Security.Cryptography;
 using System.Text;
+using System.Threading;
 using PerfectNumbers.Core.Gpu;
 using ILGPU;
 using ILGPU.Runtime;
-using System.Runtime.CompilerServices;
 
 namespace PerfectNumbers.Core;
 
@@ -146,6 +151,342 @@ public class MersenneDivisorCycles
 		return order;
 	}
 
+    public static bool TryCalculateCycleLengthForExponent(
+        ulong divisor,
+        ulong exponent,
+        Dictionary<ulong, FactorCacheEntry>? factorCache,
+        out ulong cycleLength)
+    {
+        cycleLength = 0UL;
+
+        if ((divisor & (divisor - 1UL)) == 0UL)
+        {
+            cycleLength = 1UL;
+            return true;
+        }
+
+        if (divisor <= 3UL || exponent <= 1UL)
+        {
+            cycleLength = CalculateCycleLength(divisor);
+            return true;
+        }
+
+        ulong phi = divisor - 1UL;
+        if (phi == 0UL)
+        {
+            return false;
+        }
+
+        int twoCount = BitOperations.TrailingZeroCount(phi);
+        ulong reducedPhi = phi >> twoCount;
+        if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
+        {
+            return false;
+        }
+
+        ulong k = reducedPhi / exponent;
+        if (k == 0UL)
+        {
+            return false;
+        }
+
+        Dictionary<ulong, int> factorCounts = new(8);
+        if (twoCount > 0)
+        {
+            factorCounts[2UL] = twoCount;
+        }
+
+        if (!AccumulateFactors(exponent, factorCache, factorCounts) || !AccumulateFactors(k, factorCache, factorCounts))
+        {
+            return false;
+        }
+
+        cycleLength = ReduceOrder(divisor, phi, factorCounts);
+        return true;
+    }
+
+    private static bool AccumulateFactors(
+        ulong value,
+        Dictionary<ulong, FactorCacheEntry>? cache,
+        Dictionary<ulong, int> counts)
+    {
+        if (value <= 1UL)
+        {
+            return true;
+        }
+
+        if (!TryGetFactorization(value, cache, out FactorCacheEntry factorization))
+        {
+            return false;
+        }
+
+        ulong[] primes = factorization.Primes;
+        byte[] exponents = factorization.Exponents;
+        int length = factorization.Count;
+
+        for (int i = 0; i < length; i++)
+        {
+            AddFactor(counts, primes[i], exponents[i]);
+        }
+
+        return true;
+    }
+
+    private static bool TryGetFactorization(
+        ulong value,
+        Dictionary<ulong, FactorCacheEntry>? cache,
+        out FactorCacheEntry factorization)
+    {
+        if (cache is not null && cache.TryGetValue(value, out factorization))
+        {
+            return true;
+        }
+
+        Dictionary<ulong, int> scratch = new(8);
+        if (!TryFactorIntoCountsInternal(value, scratch))
+        {
+            factorization = default!;
+            return false;
+        }
+
+        int count = scratch.Count;
+        ulong[] primes = new ulong[count];
+        byte[] exponents = new byte[count];
+        int index = 0;
+        foreach (KeyValuePair<ulong, int> entry in scratch)
+        {
+            primes[index] = entry.Key;
+            exponents[index] = checked((byte)entry.Value);
+            index++;
+        }
+
+        Array.Sort(primes, exponents);
+        factorization = new FactorCacheEntry(primes, exponents, count);
+
+        if (cache is not null)
+        {
+            cache[value] = factorization;
+        }
+
+        return true;
+    }
+
+    private static bool TryFactorIntoCountsInternal(ulong value, Dictionary<ulong, int> counts)
+    {
+        if (value <= 1UL)
+        {
+            return true;
+        }
+
+        ulong remaining = value;
+        uint[] smallPrimes = PrimesGenerator.SmallPrimes;
+        ulong[] smallPrimesSquared = PrimesGenerator.SmallPrimesPow2;
+        int smallLength = smallPrimes.Length;
+
+        for (int i = 0; i < smallLength; i++)
+        {
+            ulong prime = smallPrimes[i];
+            if (smallPrimesSquared[i] > remaining)
+            {
+                break;
+            }
+
+            while ((remaining % prime) == 0UL)
+            {
+                AddFactor(counts, prime);
+                remaining /= prime;
+            }
+        }
+
+        if (remaining == 1UL)
+        {
+            return true;
+        }
+
+        if (PrimeTester.IsPrimeInternal(remaining, CancellationToken.None))
+        {
+            AddFactor(counts, remaining);
+            return true;
+        }
+
+        ulong factor = PollardRho64(remaining);
+        if (factor == 0UL || factor == remaining)
+        {
+            return false;
+        }
+
+        ulong quotient = remaining / factor;
+        return TryFactorIntoCountsInternal(factor, counts) && TryFactorIntoCountsInternal(quotient, counts);
+    }
+
+    private static void AddFactor(Dictionary<ulong, int> counts, ulong factor)
+    {
+        if (counts.TryGetValue(factor, out int existing))
+        {
+            counts[factor] = existing + 1;
+        }
+        else
+        {
+            counts[factor] = 1;
+        }
+    }
+
+    private static void AddFactor(Dictionary<ulong, int> counts, ulong factor, int multiplicity)
+    {
+        if (counts.TryGetValue(factor, out int existing))
+        {
+            counts[factor] = existing + multiplicity;
+        }
+        else
+        {
+            counts[factor] = multiplicity;
+        }
+    }
+
+    private static ulong ReduceOrder(ulong divisor, ulong initialOrder, Dictionary<ulong, int> factorCounts)
+    {
+        if (factorCounts.Count == 0)
+        {
+            return initialOrder;
+        }
+
+        MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor);
+        int count = factorCounts.Count;
+        ulong[] primes = ArrayPool<ulong>.Shared.Rent(count);
+        try
+        {
+            int index = 0;
+            foreach (KeyValuePair<ulong, int> entry in factorCounts)
+            {
+                primes[index] = entry.Key;
+                index++;
+            }
+
+            Array.Sort(primes, 0, count);
+
+            ulong order = initialOrder;
+            for (int i = 0; i < count; i++)
+            {
+                ulong prime = primes[i];
+                int multiplicity = factorCounts[prime];
+
+                for (int iteration = 0; iteration < multiplicity; iteration++)
+                {
+                    if (order % prime != 0UL)
+                    {
+                        break;
+                    }
+
+                    ulong candidate = order / prime;
+                    if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL)
+                    {
+                        order = candidate;
+                        continue;
+                    }
+
+                    break;
+                }
+            }
+
+            return order;
+        }
+        finally
+        {
+            ArrayPool<ulong>.Shared.Return(primes, clearArray: false);
+        }
+    }
+
+    private static ulong PollardRho64(ulong n)
+    {
+        if ((n & 1UL) == 0UL)
+        {
+            return 2UL;
+        }
+
+        while (true)
+        {
+            ulong c = (NextRandomUInt64() % (n - 1UL)) + 1UL;
+            ulong x = (NextRandomUInt64() % (n - 2UL)) + 2UL;
+            ulong y = x;
+            ulong d = 1UL;
+
+            while (d == 1UL)
+            {
+                x = AdvancePolynomial(x, c, n);
+                y = AdvancePolynomial(y, c, n);
+                y = AdvancePolynomial(y, c, n);
+
+                ulong diff = x > y ? x - y : y - x;
+                d = BinaryGcd(diff, n);
+            }
+
+            if (d != n)
+            {
+                return d;
+            }
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus)
+    {
+        UInt128 value = (UInt128)x * x + c;
+        return (ulong)(value % modulus);
+    }
+
+    private static ulong BinaryGcd(ulong a, ulong b)
+    {
+        if (a == 0UL)
+        {
+            return b;
+        }
+
+        if (b == 0UL)
+        {
+            return a;
+        }
+
+        int shift = BitOperations.TrailingZeroCount(a | b);
+        a >>= BitOperations.TrailingZeroCount(a);
+
+        while (true)
+        {
+            b >>= BitOperations.TrailingZeroCount(b);
+            if (a > b)
+            {
+                (a, b) = (b, a);
+            }
+
+            b -= a;
+            if (b == 0UL)
+            {
+                return a << shift;
+            }
+        }
+    }
+
+    private static ulong NextRandomUInt64()
+    {
+        Span<byte> buffer = stackalloc byte[8];
+        RandomNumberGenerator.Fill(buffer);
+        return BinaryPrimitives.ReadUInt64LittleEndian(buffer);
+    }
+
+    public readonly struct FactorCacheEntry
+    {
+        public FactorCacheEntry(ulong[] primes, byte[] exponents, int count)
+        {
+            Primes = primes;
+            Exponents = exponents;
+            Count = count;
+        }
+        public ulong[] Primes { get; }
+
+        public byte[] Exponents { get; }
+
+        public int Count { get; }
+    }
+
 	public static (long nextPosition, long completeCount) FindLast(string path)
 	{
 		using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false);

PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs

@@ -237,6 +237,7 @@ public static class MersenneNumberDivisorByDivisorTester
                                 Span<byte> hitsSpan = default;
                                 int hitIndex = 0;
                                 int index;
+                                Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
                                 ulong divisorCycle = 0UL;
 
                                 try
@@ -254,7 +255,6 @@ public static class MersenneNumberDivisorByDivisorTester
                                                 divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted, ref localDivisorCursor, ref localDivisorsRemaining);
                                                 if (divisor != 0UL)
                                                 {
-                                                        divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
 
                                                         activeCount = BuildPrimeBuffer(
                                                                 divisor,
@@ -281,6 +281,19 @@ public static class MersenneNumberDivisorByDivisorTester
                                                         {
                                                                 continue;
                                                         }
+                                                        divisorCycle = 0UL;
+                                                        if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+                                                        {
+                                                                divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
+                                                        }
+                                                        else
+                                                        {
+                                                                factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(16);
+                                                                if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, primeBuffer[0], factorCache, out divisorCycle) || divisorCycle == 0UL)
+                                                                {
+                                                                        divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor);
+                                                                }
+                                                        }
 
                                                         hitsSpan = hitsBuffer.AsSpan(0, activeCount);
                                                         hitsSpan.Clear();

PerfectNumbers.Core/MontgomeryDivisorData.cs

@@ -4,42 +4,42 @@ using System.Runtime.CompilerServices;
 
 namespace PerfectNumbers.Core;
 
-public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo)
+public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared)
 {
-    public ulong Modulus { get; } = modulus;
+    public readonly ulong Modulus = modulus;
 
-    public ulong NPrime { get; } = nPrime;
+    public readonly ulong NPrime = nPrime;
 
-    public ulong MontgomeryOne { get; } = montgomeryOne;
+    public readonly ulong MontgomeryOne = montgomeryOne;
 
-    public ulong MontgomeryTwo { get; } = montgomeryTwo;
+    public readonly ulong MontgomeryTwo = montgomeryTwo;
+    public readonly ulong MontgomeryTwoSquared = montgomeryTwoSquared;
 }
 
 internal static class MontgomeryDivisorDataCache
 {
     private static readonly ConcurrentDictionary<ulong, MontgomeryDivisorData> Cache = new();
 
-    public static MontgomeryDivisorData Get(ulong modulus)
-    {
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL);
-        }
-
-        return Cache.GetOrAdd(modulus, static m => Create(m));
-    }
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static MontgomeryDivisorData Get(ulong modulus) => Cache.GetOrAdd(modulus, static m => Create(m));
 
     private static MontgomeryDivisorData Create(ulong modulus)
     {
+        ulong nPrime = ComputeMontgomeryNPrime(modulus);
+        ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus);
+        ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus);
+        ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime);
+
         return new MontgomeryDivisorData(
             modulus,
-            ComputeMontgomeryNPrime(modulus),
-            ComputeMontgomeryResidue(1UL, modulus),
-            ComputeMontgomeryResidue(2UL, modulus));
+            nPrime,
+            montgomeryOne,
+            montgomeryTwo,
+            montgomeryTwoSquared);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus);
+    private static ulong ComputeMontgomeryResidue(UInt128 value, ulong modulus) => (ulong)(value % modulus);
 
     private static ulong ComputeMontgomeryNPrime(ulong modulus)
     {

PerfectNumbers.Core/UInt128Numbers.cs

@@ -6,6 +6,9 @@ namespace PerfectNumbers.Core
         public static readonly UInt128 Four = 4UL;
         public static readonly UInt128 SixtyFour = 64UL;
         public static readonly UInt128 OneHundredTwentySeven = 127UL;
+        public static readonly UInt128 OneShiftedLeft64 = UInt128.One << 64;
+        public static readonly UInt128 OneShiftedLeft64x2 = Two * OneShiftedLeft64;
+
         // TODO: Populate this table from the UInt128 residue benchmarks so GPU kernels can pull precomputed constants
         // without re-deriving them at runtime.
         // TODO: Extend the table with cached Mod3/Mod5 folding constants so both CPU and GPU residue helpers can reuse

PerfectNumbers.Core/ULongExtensions.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using PerfectNumbers.Core.Gpu;
@@ -6,6 +7,8 @@ namespace PerfectNumbers.Core;
 
 public static class ULongExtensions
 {
+    private const int Pow2WindowSize = 8;
+    private const ulong Pow2WindowFallbackThreshold = 32UL;
     public static ulong CalculateOrder(this ulong q)
     {
         if (q <= 2UL)
@@ -293,25 +296,85 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryMod(this ulong exponent, in MontgomeryDivisorData divisor)
+    public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
     {
         ulong modulus = divisor.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        if (exponent == 0UL)
         {
-            return 0UL;
+            return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus;
         }
 
+        if (exponent <= Pow2WindowFallbackThreshold)
+        {
+            return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery);
+        }
+
+        int bitLength = BitOperations.Log2(exponent) + 1;
+        int windowSize = GetWindowSize(bitLength);
+        int oddPowerCount = 1 << (windowSize - 1);
+
+        ulong result = divisor.MontgomeryOne;
+        ulong nPrime = divisor.NPrime;
+        Span<ulong> oddPowers = stackalloc ulong[oddPowerCount];
+        InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowers);
+
+        int index = bitLength - 1;
+        while (index >= 0)
+        {
+            if (((exponent >> index) & 1UL) == 0UL)
+            {
+                result = result.MontgomeryMultiply(result, modulus, nPrime);
+                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.MontgomeryMultiply(result, modulus, nPrime);
+            }
+
+            ulong mask = (1UL << windowLength) - 1UL;
+            ulong windowValue = (exponent >> windowStart) & mask;
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+            ulong multiplier = oddPowers[tableIndex];
+            result = result.MontgomeryMultiply(multiplier, modulus, nPrime);
+
+            index = windowStart - 1;
+        }
+
+        if (keepMontgomery)
+        {
+            return result;
+        }
+
+        return result.MontgomeryMultiply(1UL, modulus, nPrime);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery)
+    {
+        ulong modulus = divisor.Modulus;
+        ulong nPrime = divisor.NPrime;
         ulong result = divisor.MontgomeryOne;
         ulong baseVal = divisor.MontgomeryTwo;
-        ulong nPrime = divisor.NPrime;
         ulong remainingExponent = exponent;
 
-        while (remainingExponent > 0UL)
+        while (remainingExponent != 0UL)
         {
             if ((remainingExponent & 1UL) != 0UL)
             {
-                // TODO: Switch to the ProcessEightBitWindows helper once the scalar implementation ships so CPU callers
-                // inherit the same 2× speedup that the benchmarked windowed ladder delivered for large exponents.
                 result = result.MontgomeryMultiply(baseVal, modulus, nPrime);
             }
 
@@ -324,68 +387,68 @@ public static class ULongExtensions
             baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime);
         }
 
-        return result.MontgomeryMultiply(1UL, modulus, nPrime);
+        return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong Pow2MontgomeryModMontgomery(this ulong exponent, in MontgomeryDivisorData divisor)
+    private static int GetWindowSize(int bitLength)
     {
-        ulong modulus = divisor.Modulus;
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        if (bitLength <= 6)
         {
-            return 0UL;
+            return Math.Max(bitLength, 1);
         }
 
-        ulong result = divisor.MontgomeryOne;
-        ulong baseVal = divisor.MontgomeryTwo;
-        ulong remainingExponent = exponent;
-
-        while (remainingExponent > 0UL)
+        if (bitLength <= 23)
         {
-            if ((remainingExponent & 1UL) != 0UL)
-            {
-                // TODO: Route this Montgomery-domain variant through the same windowed pow2 helper once available so delta
-                // stepping matches the ProcessEightBitWindows gains highlighted by GpuPow2ModBenchmarks.
-                result = result.MontgomeryMultiply(baseVal, modulus, divisor.NPrime);
-            }
-
-            remainingExponent >>= 1;
-            if (remainingExponent == 0UL)
-            {
-                break;
-            }
-
-            baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, divisor.NPrime);
+            return 4;
         }
 
-        return result;
+        if (bitLength <= 79)
+        {
+            return 5;
+        }
+
+        if (bitLength <= 239)
+        {
+            return 6;
+        }
+
+        if (bitLength <= 671)
+        {
+            return 7;
+        }
+
+        return Pow2WindowSize;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span<ulong> oddPowers)
+    {
+        oddPowers[0] = divisor.MontgomeryTwo;
+        if (oddPowers.Length == 1)
+        {
+            return;
+        }
+
+        ulong square = divisor.MontgomeryTwoSquared;
+        for (int i = 1; i < oddPowers.Length; i++)
+        {
+            oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(square, modulus, nPrime);
+        }
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor)
     {
-        if (cycleLength == 0UL)
-        {
-            return exponent.Pow2MontgomeryMod(divisor);
-        }
-
-        // TODO: Swap this modulo with the divisor-cycle remainder helper from Pow2MontgomeryModCycleComputationBenchmarks so we
-        // advance the exponent using cached cycle deltas instead of paying for a fresh `%` on every query.
         ulong rotationCount = exponent % cycleLength;
-        return rotationCount.Pow2MontgomeryModFromCycleRemainder(divisor);
+        return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor)
     {
-        if (reducedExponent == 0UL)
-        {
-            return 1UL;
-        }
-
-        return reducedExponent.Pow2MontgomeryMod(divisor);
+        return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false);
     }
-
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
     {