Merge pull request #166 from mlemanczyk/codex/locate-montgomery-arithmetic-usage-in-gpu-kernels-ir9oyy

Align prime-order GPU pow2 kernel with plain modular helper

EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs

@@ -158,36 +158,16 @@ public class MulMod64Benchmarks
     }
 
     /// <summary>
-    /// Uses the GPU-friendly <see cref="ULongExtensions.MulMod64GpuCompatible(ulong, ulong, ulong)"/> helper; dense inputs now
-    /// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at
-    /// 3.59 ns.
+    /// Uses the production <see cref="ULongExtensions.MulMod64(ulong, ulong, ulong)"/> helper shared by CPU and GPU callers.
     /// </summary>
     /// <remarks>
-    /// Observed means: CrossWordBlend 24.422 ns (6.86×), MixedBitPattern 61.850 ns, NearFullRange 7.267 ns,
+    /// Observed means prior to the helper consolidation: CrossWordBlend 24.422 ns, MixedBitPattern 61.850 ns, NearFullRange 7.267 ns,
     /// PrimeSizedModulus 55.114 ns, SparseOperands 7.261 ns, ZeroOperands 3.587 ns.
     /// </remarks>
     [Benchmark]
-    public ulong GpuCompatibleBaseline()
+    public ulong UnifiedMulMod()
     {
-        // TODO: Callers in production should migrate to ULongExtensions.MulMod64 where GPU parity
-        // is not required; that baseline stayed 6–17× faster for dense 64-bit operands.
-        return Input.Left.MulMod64GpuCompatible(Input.Right, Input.Modulus);
-    }
-
-    /// <summary>
-    /// Uses the deferred native-modulo helper <see cref="ULongExtensions.MulMod64GpuCompatibleDeferred(ulong, ulong, ulong)"/>;
-    /// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns.
-    /// </summary>
-    /// <remarks>
-    /// Observed means: CrossWordBlend 146.559 ns (40.63×), MixedBitPattern 255.555 ns, NearFullRange 21.390 ns,
-    /// PrimeSizedModulus 294.486 ns, SparseOperands 18.182 ns, ZeroOperands 2.009 ns.
-    /// </remarks>
-    [Benchmark]
-    public ulong GpuCompatibleDeferred()
-    {
-        // TODO: Retire the deferred GPU shim from runtime paths and keep it only for benchmarks;
-        // ULongExtensions.MulMod64 avoids the 6×–82× slowdown on real-world operand mixes.
-        return Input.Left.MulMod64GpuCompatibleDeferred(Input.Right, Input.Modulus);
+        return Input.Left.MulMod64(Input.Right, Input.Modulus);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]

EvenPerfectBitScanner/Program.cs

@@ -1030,7 +1030,7 @@ internal static class Program
         Console.WriteLine("  --use-order            test primality via q order");
         Console.WriteLine("  --workaround-mod       avoid '%' operator on the GPU");
         // mod-automaton removed
-        Console.WriteLine("  --filter-p=<path>      process only p from previous run results (required for --mersenne=bydivisor)");
+        Console.WriteLine("  --filter-p=<path>      process only p from previous run results (required for --mersenne=bydivisor; primes must be strictly increasing)");
         Console.WriteLine("  --write-batch-size=<value> overwrite frequency of disk writes (default 100 lines)");
         Console.WriteLine("  --gcd-filter           enable early sieve based on GCD");
         Console.WriteLine("  --help, -help, --?, -?, /?   show this help message");

PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs

@@ -1,5 +1,6 @@
 using System;
 using FluentAssertions;
+using PerfectNumbers.Core;
 using PerfectNumbers.Core.Cpu;
 using Xunit;
 
@@ -7,20 +8,23 @@ namespace PerfectNumbers.Core.Tests;
 
 public class MersenneNumberDivisorCpuTesterTests
 {
+    private const ulong LargeDivisor1 = 4_000_169UL;
+    private const ulong LargeDivisor2 = 4_000_561UL;
+
     [Fact]
     [Trait("Category", "Fast")]
     public void ByDivisor_tester_tracks_divisors_across_primes()
     {
         var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(11UL);
+        tester.ConfigureFromMaxPrime(83UL);
 
-        tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue();
+        tester.IsPrime(59UL, out bool divisorsExhausted).Should().BeTrue();
         divisorsExhausted.Should().BeTrue();
 
-        tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue();
+        tester.IsPrime(61UL, out divisorsExhausted).Should().BeTrue();
         divisorsExhausted.Should().BeTrue();
 
-        tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse();
+        tester.IsPrime(67UL, out divisorsExhausted).Should().BeTrue();
         divisorsExhausted.Should().BeTrue();
     }
 
@@ -29,66 +33,94 @@ public class MersenneNumberDivisorCpuTesterTests
     public void ByDivisor_prime_test_limit_disables_divisor_exhaustion()
     {
         var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(11UL);
+        tester.ConfigureFromMaxPrime(83UL);
 
-        tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
+        tester.IsPrime(59UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue();
         divisorsExhausted.Should().BeFalse();
     }
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_session_checks_divisors_across_primes()
+    public void ByDivisor_session_checks_divisors_across_exponents()
     {
         var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(13UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] primes = { 5UL, 7UL, 11UL, 13UL };
-        byte[] hits = new byte[primes.Length];
+        ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL, 333_380UL };
+        byte[] hits = new byte[exponents.Length];
 
-        ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL, MontgomeryDivisorData.FromModulus(23UL));
-        session.CheckDivisor(23UL, MontgomeryDivisorData.FromModulus(23UL), cycle23, primes, hits);
-        hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
+        MontgomeryDivisorData divisor1 = MontgomeryDivisorData.FromModulus(LargeDivisor1);
+        ulong cycle1 = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor1, divisor1);
+        session.CheckDivisor(LargeDivisor1, divisor1, cycle1, exponents, hits);
+
+        byte[] expected1 = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected1[i] = exponents[i].Pow2ModWindowedCpu(divisor1.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        hits.Should().Equal(expected1);
 
         Array.Fill(hits, (byte)0);
-        ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL, MontgomeryDivisorData.FromModulus(31UL));
-        session.CheckDivisor(31UL, MontgomeryDivisorData.FromModulus(31UL), cycle31, primes, hits);
-        hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 });
+        MontgomeryDivisorData divisor2 = MontgomeryDivisorData.FromModulus(LargeDivisor2);
+        ulong cycle2 = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor2, divisor2);
+        session.CheckDivisor(LargeDivisor2, divisor2, cycle2, exponents, hits);
+
+        byte[] expected2 = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected2[i] = exponents[i].Pow2ModWindowedCpu(divisor2.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        hits.Should().Equal(expected2);
     }
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
+    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor1_as_composite()
     {
         var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(43UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
+        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(LargeDivisor1);
+        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor1, divisorData);
+
+        ulong[] exponents = { 59UL, 67UL, cycle, cycle * 2UL, cycle + 1UL };
         byte[] hits = new byte[exponents.Length];
+        session.CheckDivisor(LargeDivisor1, divisorData, cycle, exponents, hits);
 
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
-        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
-        session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
+        byte[] expected = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
 
-        hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+        hits.Should().Equal(expected);
     }
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
+    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor2_as_composite()
     {
         var tester = new MersenneNumberDivisorByDivisorCpuTester();
-        tester.ConfigureFromMaxPrime(61UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
+        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(LargeDivisor2);
+        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(LargeDivisor2, divisorData);
+
+        ulong[] exponents = { 61UL, 73UL, cycle, cycle * 2UL, cycle - 1UL };
         byte[] hits = new byte[exponents.Length];
+        session.CheckDivisor(LargeDivisor2, divisorData, cycle, exponents, hits);
 
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
-        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
-        session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
+        byte[] expected = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
 
-        hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+        hits.Should().Equal(expected);
     }
 }

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -12,12 +12,12 @@ namespace PerfectNumbers.Core.Tests;
 public class MersenneNumberDivisorGpuTesterTests
 {
     [Theory]
-    [InlineData(33UL, 7UL, true)]
-    [InlineData(35UL, 31UL, true)]
-    [InlineData(37UL, 223UL, true)]
-    [InlineData(33UL, 13UL, false)]
-    [InlineData(35UL, 73UL, false)]
-    [InlineData(37UL, 227UL, false)]
+    [InlineData(117_652UL, 4_000_169UL, true)]
+    [InlineData(333_380UL, 4_000_561UL, true)]
+    [InlineData(59UL, 4_000_169UL, false)]
+    [InlineData(61UL, 4_000_561UL, false)]
+    [InlineData(333_379UL, 4_000_561UL, false)]
+    [InlineData(117_653UL, 4_000_169UL, false)]
     [Trait("Category", "Fast")]
     public void IsDivisible_returns_expected(ulong exponent, ulong divisor, bool expected)
     {
@@ -152,23 +152,39 @@ public class MersenneNumberDivisorGpuTesterTests
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_session_checks_divisors_across_primes()
+    public void ByDivisor_session_checks_divisors_across_exponents()
     {
         var tester = new MersenneNumberDivisorByDivisorGpuTester();
-        tester.ConfigureFromMaxPrime(43UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] primes = { 31UL, 37UL, 41UL, 43UL };
-        byte[] hits = new byte[primes.Length];
+        ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL, 333_380UL };
+        byte[] hits = new byte[exponents.Length];
 
-        ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
-        session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
-        hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0 });
+        MontgomeryDivisorData divisor1 = MontgomeryDivisorData.FromModulus(4_000_169UL);
+        ulong cycle1 = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor1);
+        session.CheckDivisor(4_000_169UL, divisor1, cycle1, exponents, hits);
+
+        byte[] expected1 = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected1[i] = exponents[i].Pow2ModWindowedCpu(divisor1.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        hits.Should().Equal(expected1);
 
         Array.Fill(hits, (byte)0);
-        ulong cycle13367 = MersenneDivisorCycles.CalculateCycleLength(13367UL, MontgomeryDivisorData.FromModulus(13367UL));
-        session.CheckDivisor(13367UL, MontgomeryDivisorData.FromModulus(13367UL), cycle13367, primes, hits);
-        hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 });
+        MontgomeryDivisorData divisor2 = MontgomeryDivisorData.FromModulus(4_000_561UL);
+        ulong cycle2 = MersenneDivisorCycles.CalculateCycleLength(4_000_561UL, divisor2);
+        session.CheckDivisor(4_000_561UL, divisor2, cycle2, exponents, hits);
+
+        byte[] expected2 = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected2[i] = exponents[i].Pow2ModWindowedCpu(divisor2.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        hits.Should().Equal(expected2);
     }
 
     [Fact]
@@ -180,122 +196,126 @@ public class MersenneNumberDivisorGpuTesterTests
             GpuBatchSize = 2,
         };
 
-        tester.ConfigureFromMaxPrime(47UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
-        byte[] hits = new byte[primes.Length];
+        ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL };
+        byte[] hits = new byte[exponents.Length];
 
-        ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
-        session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
+        MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(4_000_169UL);
+        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor);
+        session.CheckDivisor(4_000_169UL, divisor, cycle, exponents, hits);
 
-        hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
+        byte[] expected = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected[i] = exponents[i].Pow2ModWindowedCpu(divisor.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        hits.Should().Equal(expected);
     }
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_seven_as_composite()
+    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor1_as_composite()
     {
         var tester = new MersenneNumberDivisorByDivisorGpuTester
         {
             GpuBatchSize = 5,
         };
-        tester.ConfigureFromMaxPrime(43UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] exponents = { 6UL, 7UL, 9UL, 10UL, 12UL };
+        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_169UL);
+        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisorData);
+
+        ulong[] exponents = { 59UL, 67UL, cycle, cycle * 2UL, cycle + 1UL };
         byte[] hits = new byte[exponents.Length];
 
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(7UL);
-        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(7UL, divisorData);
-
-        session.CheckDivisor(7UL, divisorData, cycle, exponents, hits);
+        session.CheckDivisor(4_000_169UL, divisorData, cycle, exponents, hits);
 
         var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession)
             .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!;
         ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!;
         Span<ulong> residues = hostBuffer.AsSpan(0, exponents.Length);
-        byte[] computedHits = new byte[exponents.Length];
+        byte[] expected = new byte[exponents.Length];
         for (int i = 0; i < residues.Length; i++)
         {
-            computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0;
+            expected[i] = residues[i] == 1UL ? (byte)1 : (byte)0;
         }
 
-        computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+        byte[] direct = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            direct[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        expected.Should().Equal(direct);
+        hits.Should().Equal(direct);
     }
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_eleven_as_composite()
+    public void ByDivisor_session_marks_mersenne_numbers_divisible_by_large_divisor2_as_composite()
     {
         var tester = new MersenneNumberDivisorByDivisorGpuTester
         {
             GpuBatchSize = 5,
         };
-        tester.ConfigureFromMaxPrime(61UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] exponents = { 10UL, 11UL, 20UL, 21UL, 30UL };
+        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_561UL);
+        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_561UL, divisorData);
+
+        ulong[] exponents = { 61UL, 73UL, cycle, cycle * 2UL, cycle - 1UL };
         byte[] hits = new byte[exponents.Length];
 
-        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(11UL);
-        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(11UL, divisorData);
-
-        session.CheckDivisor(11UL, divisorData, cycle, exponents, hits);
+        session.CheckDivisor(4_000_561UL, divisorData, cycle, exponents, hits);
 
         var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession)
             .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!;
         ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!;
         Span<ulong> residues = hostBuffer.AsSpan(0, exponents.Length);
-        byte[] computedHits = new byte[exponents.Length];
+        byte[] expected = new byte[exponents.Length];
         for (int i = 0; i < residues.Length; i++)
         {
-            computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0;
+            expected[i] = residues[i] == 1UL ? (byte)1 : (byte)0;
         }
 
-        computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 });
+        byte[] direct = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            direct[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        expected.Should().Equal(direct);
+        hits.Should().Equal(direct);
     }
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_gpu_tester_skips_divisors_excluded_by_small_cycle_generation()
+    public void ByDivisor_gpu_tester_generates_cycles_for_large_divisors()
     {
-        var cycles = MersenneDivisorCycles.Shared;
-        var tableField = typeof(MersenneDivisorCycles).GetField("_table", BindingFlags.NonPublic | BindingFlags.Instance)!;
-        var smallCyclesField = typeof(MersenneDivisorCycles).GetField("_smallCycles", BindingFlags.NonPublic | BindingFlags.Instance)!;
+        var tester = new MersenneNumberDivisorByDivisorGpuTester();
+        tester.ConfigureFromMaxPrime(197UL);
 
-        var originalTable = (List<(ulong Divisor, ulong Cycle)>)tableField.GetValue(cycles)!;
-        var originalSmall = (ulong[]?)smallCyclesField.GetValue(cycles);
+        using var session = tester.CreateDivisorSession();
+        ulong[] exponents = { 59UL, 67UL, 117_652UL };
+        byte[] hits = new byte[exponents.Length];
 
-        try
+        MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(4_000_169UL);
+        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisorData);
+
+        session.CheckDivisor(4_000_169UL, divisorData, cycle, exponents, hits);
+
+        byte[] expected = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
         {
-            var patchedTable = new List<(ulong Divisor, ulong Cycle)>();
-            ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1];
-
-            MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(191UL);
-            patchedSmall[191] = MersenneDivisorCycles.CalculateCycleLength(191UL, divisorData);
-
-            tableField.SetValue(cycles, patchedTable);
-            smallCyclesField.SetValue(cycles, patchedSmall);
-            DivisorCycleCache.Shared.RefreshSnapshot();
-
-            var tester = new MersenneNumberDivisorByDivisorGpuTester();
-            tester.ConfigureFromMaxPrime(19UL);
-
-            typeof(MersenneNumberDivisorByDivisorGpuTester)
-                .GetField("_divisorLimit", BindingFlags.NonPublic | BindingFlags.Instance)!
-                .SetValue(tester, 200UL);
-
-            tester.IsPrime(19UL, out bool divisorsExhausted).Should().BeTrue();
-            divisorsExhausted.Should().BeTrue();
-        }
-        finally
-        {
-            tableField.SetValue(cycles, originalTable);
-            smallCyclesField.SetValue(cycles, originalSmall);
-            DivisorCycleCache.Shared.RefreshSnapshot();
-            GpuContextPool.DisposeAll();
+            expected[i] = exponents[i].Pow2ModWindowedCpu(divisorData.Modulus) == 1UL ? (byte)1 : (byte)0;
         }
+
+        hits.Should().Equal(expected);
     }
 
     [Fact]
@@ -307,16 +327,23 @@ public class MersenneNumberDivisorGpuTesterTests
             BatchSize = 2,
         };
 
-        tester.ConfigureFromMaxPrime(47UL);
+        tester.ConfigureFromMaxPrime(197UL);
 
         using var session = tester.CreateDivisorSession();
-        ulong[] primes = { 31UL, 37UL, 41UL, 43UL, 47UL };
-        byte[] hits = new byte[primes.Length];
+        ulong[] exponents = { 59UL, 61UL, 67UL, 71UL, 73UL, 79UL, 117_652UL };
+        byte[] hits = new byte[exponents.Length];
 
-        ulong cycle223 = MersenneDivisorCycles.CalculateCycleLength(223UL, MontgomeryDivisorData.FromModulus(223UL));
-        session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), cycle223, primes, hits);
+        MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(4_000_169UL);
+        ulong cycle = MersenneDivisorCycles.CalculateCycleLength(4_000_169UL, divisor);
+        session.CheckDivisor(4_000_169UL, divisor, cycle, exponents, hits);
 
-        hits.Should().ContainInOrder(new byte[] { 0, 1, 0, 0, 0 });
+        byte[] expected = new byte[exponents.Length];
+        for (int i = 0; i < exponents.Length; i++)
+        {
+            expected[i] = exponents[i].Pow2ModWindowedCpu(divisor.Modulus) == 1UL ? (byte)1 : (byte)0;
+        }
+
+        hits.Should().Equal(expected);
     }
 }
 

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -633,17 +633,23 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
     private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit)
     {
-        if (prime <= 1UL)
-        {
-            return 0UL;
-        }
+        // The by-divisor scanner never emits exponents below the first odd prime, so the defensive guard stays commented out to
+        // avoid branching on production workloads.
+        // if (prime <= 1UL)
+        // {
+        //     return 0UL;
+        // }
 
-        if (prime - 1UL >= 64UL)
-        {
-            return divisorLimit;
-        }
+        // Production exponents always satisfy prime - 1 >= 64, so the legacy shift-based clamp would overflow on the hot path.
+        // Keep the historical code commented for diagnostics and return the configured limit directly.
+        // if (prime - 1UL >= 64UL)
+        // {
+        //     return divisorLimit;
+        // }
+        //
+        // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
 
-        return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
+        return divisorLimit;
     }
 
     private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -652,7 +658,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         private bool _disposed;
         private ulong[]? _primeDeltas;
         private ulong[]? _cycleRemainders;
-        private ulong[]? _montgomeryDeltas;
+        private ulong[]? _residues;
         private int _bufferCapacity;
         private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator;
 
@@ -674,10 +680,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             }
 
             int length = primes.Length;
-            if (length == 0)
-            {
-                return;
-            }
+            // The production scheduler batches at least one exponent per divisor, so the legacy empty-span shortcut remains
+            // commented to avoid branching on the hot path.
+            // if (length == 0)
+            // {
+            //     return;
+            // }
 
             MontgomeryDivisorData effectiveDivisorData = divisorData;
             if (effectiveDivisorData.Modulus != divisor)
@@ -765,13 +773,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan);
 
             EnsureGpuEvaluator();
-            Span<ulong> residueSpan = _montgomeryDeltas!.AsSpan(0, length);
+            Span<ulong> residueSpan = _residues!.AsSpan(0, length);
             _gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan);
 
-            ulong montgomeryOne = divisorData.MontgomeryOne;
             for (int i = 0; i < length; i++)
             {
-                hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0;
+                hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == 1UL ? (byte)1 : (byte)0;
             }
         }
 
@@ -801,7 +808,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 return;
             }
 
-            if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null)
+            // Sessions release pooled buffers between scans; a single null check covers every array because buffers are rented and returned together.
+            if (requiredLength <= _bufferCapacity && _primeDeltas is not null)
             {
                 return;
             }
@@ -809,6 +817,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             int newCapacity = _bufferCapacity;
             if (newCapacity == 0)
             {
+                // Fresh sessions start without buffers; seed the first rent with a single-slot array.
                 newCapacity = 1;
             }
 
@@ -821,51 +830,42 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
             _primeDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
             _cycleRemainders = ArrayPool<ulong>.Shared.Rent(newCapacity);
-            _montgomeryDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
+            _residues = ArrayPool<ulong>.Shared.Rent(newCapacity);
             _bufferCapacity = newCapacity;
         }
 
         private void ReturnBuffers()
         {
-            if (_primeDeltas is not null)
+            // Buffers are lazily allocated per session; the _primeDeltas guard applies to every pooled array because sessions rent and release them together.
+            if (_primeDeltas is null)
             {
-                ArrayPool<ulong>.Shared.Return(_primeDeltas, clearArray: false);
-                _primeDeltas = null;
-            }
-
-            if (_cycleRemainders is not null)
-            {
-                ArrayPool<ulong>.Shared.Return(_cycleRemainders, clearArray: false);
-                _cycleRemainders = null;
-            }
-
-            if (_montgomeryDeltas is not null)
-            {
-                ArrayPool<ulong>.Shared.Return(_montgomeryDeltas, clearArray: false);
-                _montgomeryDeltas = null;
+                return;
             }
 
+            ArrayPool<ulong>.Shared.Return(_primeDeltas, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(_cycleRemainders!, clearArray: false);
+            ArrayPool<ulong>.Shared.Return(_residues!, clearArray: false);
+            _primeDeltas = null;
+            _cycleRemainders = null;
+            _residues = null;
             _bufferCapacity = 0;
         }
 
         private static void ComputePrimeDeltas(ReadOnlySpan<ulong> primes, Span<ulong> deltas)
         {
-            if (primes.Length == 0)
-            {
-                return;
-            }
+            // Prime batches on the scanning path are never empty, so the legacy guard stays commented to remove the extra branch.
+            // if (primes.Length == 0)
+            // {
+            //     return;
+            // }
 
             ulong previous = primes[0];
             deltas[0] = previous;
 
+            // Production batches feed strictly increasing primes; callers outside the scanning path must uphold the contract.
             for (int i = 1; i < primes.Length; i++)
             {
                 ulong current = primes[i];
-                if (current <= previous)
-                {
-                    throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing.");
-                }
-
                 deltas[i] = current - previous;
                 previous = current;
             }

PerfectNumbers.Core/DivisorCycleCache.cs

@@ -86,11 +86,10 @@ public sealed class DivisorCycleCache
                 }
 
                 cycles[i] = cached;
+                continue;
             }
-            else
-            {
-                missingBuffer[missingCount++] = i;
-            }
+
+            missingBuffer[missingCount++] = i;
         }
 
         if (missingCount != 0)
@@ -104,6 +103,7 @@ public sealed class DivisorCycleCache
             {
                 ComputeCyclesCpu(divisors, cycles, missingIndices);
             }
+
         }
 
         if (rentedMissing is not null)

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -35,7 +35,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
 
     private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeExponent(Accelerator accelerator) =>
-        _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeMontgomeryExponentKernel));
+        _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeExponentKernel));
 
     private Action<Index1D, ArrayView<ulong>, ArrayView<byte>, byte> GetRemainderDeltaKernel(Accelerator accelerator) =>
         _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<byte>, byte>(ComputeRemainderDeltasKernel));
@@ -530,46 +530,60 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
     private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
     {
-        if (maxPrime <= 1UL)
-        {
-            return 0UL;
-        }
+        // The production configuration seeds the scan with exponents above the first odd prime, so the defensive guard remains
+        // commented to document the invariant without branching.
+        // if (maxPrime <= 1UL)
+        // {
+        //     return 0UL;
+        // }
 
-        if (maxPrime - 1UL >= 64UL)
-        {
-            return ulong.MaxValue;
-        }
+        // Real workloads push maxPrime well past the 64-bit shift boundary, so keep the historical implementation commented out
+        // and return ulong.MaxValue directly.
+        // if (maxPrime - 1UL >= 64UL)
+        // {
+        //     return ulong.MaxValue;
+        // }
+        //
+        // return (1UL << (int)(maxPrime - 1UL)) - 1UL;
 
-        return (1UL << (int)(maxPrime - 1UL)) - 1UL;
+        return ulong.MaxValue;
     }
 
     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 the defensive guard remains
+        // commented to document that invariant.
+        // if (prime <= 1UL)
+        // {
+        //     return 0UL;
+        // }
 
-        if (prime - 1UL >= 64UL)
-        {
-            return divisorLimit;
-        }
+        // Production exponents are always large enough that the shift-based clamp would overflow, so skip the legacy branch and
+        // return the configured divisor limit directly.
+        // if (prime - 1UL >= 64UL)
+        // {
+        //     return divisorLimit;
+        // }
+        //
+        // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
 
-        return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit);
+        return divisorLimit;
     }
 
     private static void CheckKernel(Index1D index, ArrayView<MontgomeryDivisorData> divisors, ArrayView<ulong> exponents, ArrayView<byte> hits)
     {
         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.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = exponent.Pow2ModWindowedGpu(modulus) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -609,16 +623,20 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             }
 
             int length = primes.Length;
-            if (length == 0)
-            {
-                return;
-            }
+            // Production batches always contain at least one exponent, so the empty-span guard remains commented out to remove
+            // the redundant branch.
+            // if (length == 0)
+            // {
+            //     return;
+            // }
 
-            if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL)
-            {
-                hits.Clear();
-                return;
-            }
+            // The divisor staging path filters non-prime and even moduli before kernels launch, so keep the defensive
+            // modulus guard commented out to avoid branching in the hot loop.
+            // if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL)
+            // {
+            //     hits.Clear();
+            //     return;
+            // }
 
             ArrayView1D<ulong, Stride1D.Dense> exponentsView = _exponentsBuffer.View;
             ArrayView1D<ulong, Stride1D.Dense> resultsView = _resultsBuffer.View;
@@ -777,33 +795,34 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         mask[globalIndex] = 1;
     }
 
-    private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
+    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.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
+        results[index] = exponent.Pow2ModWindowedGpu(modulus);
     }
 
     public ulong DivisorLimit
     {
         get
         {
-            lock (_sync)
-            {
-                if (!_isConfigured)
-                {
-                    throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
-                }
+            // The by-divisor scanner always configures the tester before exposing the limit, so the guard remains
+            // commented out to document the invariant without forcing a runtime branch.
+            // if (!_isConfigured)
+            // {
+            //     throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
+            // }
 
-                return _divisorLimit;
-            }
+            return _divisorLimit;
         }
     }
 
@@ -825,10 +844,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         internal BatchResources(MersenneNumberDivisorByDivisorGpuTester owner, Accelerator accelerator, int capacity)
         {
             int cycleCapacity = DivisorCycleCache.Shared.PreferredBatchSize;
-            if (cycleCapacity <= 0)
-            {
-                cycleCapacity = 1;
-            }
+            // The shared cache always reports a positive preferred size on production workloads, so the fallback remains
+            // commented out to avoid branching during resource initialization.
+            // if (cycleCapacity <= 0)
+            // {
+            //     cycleCapacity = 1;
+            // }
 
             int actualCapacity = Math.Max(capacity, cycleCapacity);
 

PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs

@@ -104,21 +104,25 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable
     private static void BuildKernel(Index1D index, ArrayView<ulong> moduli, ArrayView<MontgomeryDivisorData> destination)
     {
         int idx = index;
-        if (idx >= moduli.Length)
-        {
-            return;
-        }
+        // The builder launches kernels with an extent equal to the modulus count, so the bounds guard remains commented
+        // out to avoid redundant branching on every thread.
+        // if (idx >= moduli.Length)
+        // {
+        //     return;
+        // }
 
         ulong modulus = moduli[idx];
-        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
-        {
-            destination[idx] = default;
-            return;
-        }
+        // All staged divisors are odd primes on the scanner path, so keep the defensive modulus guard commented out to
+        // remove the branch from the kernel.
+        // if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        // {
+        //     destination[idx] = default;
+        //     return;
+        // }
 
         ulong nPrime = ComputeMontgomeryNPrime(modulus);
         ulong montgomeryOne = ComputeMontgomeryResidue64(modulus);
-        ulong montgomeryTwo = montgomeryOne.MulMod64GpuCompatible(2UL, modulus);
+        ulong montgomeryTwo = montgomeryOne.MulMod64(2UL, modulus);
         ulong montgomeryTwoSquared = MontgomeryMultiplyGpu(montgomeryTwo, montgomeryTwo, modulus, nPrime);
 
         destination[idx] = new MontgomeryDivisorData(

PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs

@@ -103,12 +103,23 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
     {
         int idx = index;
         int length = (int)exponents.Length;
-        if (idx >= length)
-        {
-            return;
-        }
+        // Batched launches match the buffer length on production workloads, so the bounds guard stays commented out to remove
+        // the redundant comparison.
+        // if (idx >= length)
+        // {
+        //     return;
+        // }
 
         ulong exponent = exponents[idx];
-        results[idx] = exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: true);
+        ulong modulus = divisorData.Modulus;
+        // 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.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.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 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;
             }
@@ -848,10 +850,11 @@ internal static partial class PrimeOrderGpuHeuristics
         {
             ulong primeFactor = factors[i];
             int exponent = exponents[i];
-            if (primeFactor <= 1UL)
-            {
-                continue;
-            }
+            // Factorization never yields unity on the scanning path, so skip the redundant guard.
+            // if (primeFactor <= 1UL)
+            // {
+            //     continue;
+            // }
 
             for (int iteration = 0; iteration < exponent; iteration++)
             {
@@ -861,7 +864,7 @@ internal static partial class PrimeOrderGpuHeuristics
                 }
 
                 ulong reduced = order / primeFactor;
-                if (Pow2EqualsOneKernel(reduced, divisor))
+                if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
                 {
                     order = reduced;
                     continue;
@@ -887,12 +890,14 @@ internal static partial class PrimeOrderGpuHeuristics
         ArrayView1D<ulong, Stride1D.Dense> compositeStack,
         ArrayView1D<byte, Stride1D.Dense> statusOut)
     {
-        if (order == 0UL)
-        {
-            return false;
-        }
+        // Order values handed to the kernel are strictly positive on the production path, so the zero guard stays commented
+        // out to remove the redundant branch.
+        // if (order == 0UL)
+        // {
+        //     return false;
+        // }
 
-        if (!Pow2EqualsOneKernel(order, divisor))
+        if ((order).Pow2ModWindowedGpu(divisor.Modulus) != 1UL)
         {
             return false;
         }
@@ -922,10 +927,11 @@ internal static partial class PrimeOrderGpuHeuristics
         {
             ulong primeFactor = factors[i];
             int exponent = exponents[i];
-            if (primeFactor <= 1UL)
-            {
-                continue;
-            }
+            // Factorization never yields one on the by-divisor path, so the defensive continue stays commented out.
+            // if (primeFactor <= 1UL)
+            // {
+            //     continue;
+            // }
 
             ulong reduced = order;
             for (int iteration = 0; iteration < exponent; iteration++)
@@ -936,7 +942,7 @@ internal static partial class PrimeOrderGpuHeuristics
                 }
 
                 reduced /= primeFactor;
-                if (Pow2EqualsOneKernel(reduced, divisor))
+                if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
                 {
                     return false;
                 }
@@ -998,28 +1004,32 @@ internal static partial class PrimeOrderGpuHeuristics
         long candidateCapacity = candidates.Length;
         int candidateLimit = candidateCapacity < HeuristicCandidateLimit ? (int)candidateCapacity : HeuristicCandidateLimit;
         int candidateCount = BuildCandidatesKernel(order, workFactors, workExponents, factorCount, candidates, stackIndex, stackExponent, stackProduct, candidateLimit);
-        if (candidateCount == 0)
-        {
-            return false;
-        }
+        // Candidate enumeration always yields at least one order on the configured workloads, so the zero-count guard remains
+        // commented out to eliminate the redundant branch.
+        // if (candidateCount == 0)
+        // {
+        //     return false;
+        // }
 
         SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount);
 
         int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks;
-        if (powBudget <= 0)
-        {
-            powBudget = candidateCount;
-        }
+        // The heuristics always reserve at least one check, so the zero-budget fallback remains commented out.
+        // if (powBudget <= 0)
+        // {
+        //     powBudget = candidateCount;
+        // }
 
         int powUsed = 0;
 
         for (int i = 0; i < candidateCount && powUsed < powBudget; i++)
         {
             ulong candidate = candidates[i];
-            if (candidate <= 1UL)
-            {
-                continue;
-            }
+            // Candidate generation skips 1 by construction, so the defensive continue stays commented out.
+            // if (candidate <= 1UL)
+            // {
+            //     continue;
+            // }
 
             if (powUsed >= powBudget)
             {
@@ -1027,7 +1037,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             powUsed++;
-            if (!Pow2EqualsOneKernel(candidate, divisor))
+            if ((candidate).Pow2ModWindowedGpu(divisor.Modulus) != 1UL)
             {
                 continue;
             }
@@ -1307,10 +1317,11 @@ internal static partial class PrimeOrderGpuHeuristics
         {
             ulong primeFactor = factors[i];
             int exponent = exponents[i];
-            if (primeFactor <= 1UL)
-            {
-                continue;
-            }
+            // Factorization never yields one for candidate orders, so keep the guard commented out to avoid re-checking inputs.
+            // if (primeFactor <= 1UL)
+            // {
+            //     continue;
+            // }
 
             ulong reduced = candidate;
             for (int iteration = 0; iteration < exponent; iteration++)
@@ -1327,7 +1338,7 @@ internal static partial class PrimeOrderGpuHeuristics
                 }
 
                 powUsed++;
-                if (Pow2EqualsOneKernel(reduced, divisor))
+                if ((reduced).Pow2ModWindowedGpu(divisor.Modulus) == 1UL)
                 {
                     return false;
                 }
@@ -1349,17 +1360,21 @@ internal static partial class PrimeOrderGpuHeuristics
         ArrayView1D<ulong, Stride1D.Dense> compositeStack,
         ArrayView1D<byte, Stride1D.Dense> statusOut)
     {
-        if (initial <= 1UL)
-        {
-            return true;
-        }
+        // Composite candidates sent here are always greater than one during production scans, so the trivial-return branch stays
+        // commented out.
+        // if (initial <= 1UL)
+        // {
+        //     return true;
+        // }
 
         int stackCapacity = (int)compositeStack.Length;
-        if (stackCapacity <= 0)
-        {
-            statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow;
-            return false;
-        }
+        // Pollard-Rho always reserves a positive stack on production workloads, so the overflow guard remains commented out
+        // to avoid branching here.
+        // if (stackCapacity <= 0)
+        // {
+        //     statusOut[0] = (byte)PrimeOrderKernelStatus.PollardOverflow;
+        //     return false;
+        // }
 
         int stackTop = 0;
         compositeStack[stackTop] = initial;
@@ -1369,10 +1384,12 @@ internal static partial class PrimeOrderGpuHeuristics
         {
             stackTop--;
             ulong composite = compositeStack[stackTop];
-            if (composite <= 1UL)
-            {
-                continue;
-            }
+            // Composite values enqueued for factoring are always greater than one on production workloads, so the guard stays
+            // commented out to avoid extra branching.
+            // if (composite <= 1UL)
+            // {
+            //     continue;
+            // }
 
             if (!PeelSmallPrimesKernel(
                     composite,
@@ -1580,9 +1597,7 @@ internal static partial class PrimeOrderGpuHeuristics
             bLocal >>= BitOperations.TrailingZeroCount(bLocal);
             if (aLocal > bLocal)
             {
-                ulong temp = aLocal;
-                aLocal = bLocal;
-                bLocal = temp;
+                (aLocal, bLocal) = (bLocal, aLocal);
             }
 
             bLocal -= aLocal;
@@ -1593,11 +1608,6 @@ internal static partial class PrimeOrderGpuHeuristics
         }
     }
 
-    private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
-    {
-        return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
-    }
-
     private static ulong CalculateByDoublingKernel(ulong prime)
     {
         ulong order = 1UL;
@@ -1619,10 +1629,11 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> remainders)
     {
-        if (exponents.Length == 0)
-        {
-            return GpuPow2ModStatus.Success;
-        }
+        // Residue batches are never empty on production workloads, so skip the zero-length success shortcut.
+        // if (exponents.Length == 0)
+        // {
+        //     return GpuPow2ModStatus.Success;
+        // }
 
         if (remainders.Length < exponents.Length)
         {
@@ -1698,19 +1709,22 @@ 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);
         }
     }
 
     private static bool TryComputeOnGpuWide(ReadOnlySpan<UInt128> exponents, UInt128 prime, Span<UInt128> results)
     {
         int length = exponents.Length;
-        if (length == 0)
-        {
-            return true;
-        }
+        // Production launches always carry work, so keep the empty-batch shortcut commented out.
+        // if (length == 0)
+        // {
+        //     return true;
+        // }
 
         GpuUInt128[]? rentedExponents = null;
         GpuUInt128[]? rentedResults = null;
@@ -1764,6 +1778,8 @@ internal static partial class PrimeOrderGpuHeuristics
         }
         finally
         {
+            // Small batches use the stack-allocated fast path, leaving the rented arrays null. Retain the guards so pooled
+            // buffers are only returned when they were actually leased.
             if (rentedExponents is not null)
             {
                 ArrayPool<GpuUInt128>.Shared.Return(rentedExponents, clearArray: false);
@@ -1788,16 +1804,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 =>
@@ -1810,58 +1816,17 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
-        ulong exponent = exponents[index];
         ulong modulus = divisor.Modulus;
-		// ulong nPrime = divisor.NPrime;
-		// ulong result = divisor.MontgomeryOne;
-		// ulong baseValue = divisor.MontgomeryTwo;
-		// ulong remainingExponent = exponent;
+        // 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;
+        // }
 
-		// if (remainingExponent == 0UL)
-		// {
-		//     remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
-		//     return;
-		// }
-
-		// while (remainingExponent != 0UL)
-		// {
-		//     if ((remainingExponent & 1UL) != 0UL)
-		//     {
-		//         result = result.MontgomeryMultiply(baseValue, modulus, nPrime);
-		//     }
-
-		//     remainingExponent >>= 1;
-		//     if (remainingExponent == 0UL)
-		//     {
-		//         break;
-		//     }
-
-		//     // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term.
-		//     baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime);
-		// }
-
-		// remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime);
-
-		// Plain modular exponentiation without Montgomery reduction kept for reference.
-		ulong plainResult = 1UL;// % modulus;
-		ulong plainBase = 2UL;// % modulus;
-		while (exponent != 0UL)
-		{
-			if ((exponent & 1UL) != 0UL)
-			{
-				plainResult = (plainResult * plainBase) % modulus;
-			}
-
-			exponent >>= 1;
-			if (exponent == 0UL)
-			{
-				break;
-			}
-
-			plainBase = (plainBase * plainBase) % modulus;
-		}
-		
-        remainders[index] = plainResult;
+        ulong exponent = exponents[index];
+        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

@@ -74,23 +74,27 @@ public class MersenneDivisorCycles
             return true;
         }
 
-        if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
-        {
-            return false;
-        }
+        // Exponent and divisor validation happens during scheduling, so the defensive guard remains commented out to keep the
+        // hot path branch-free.
+        // if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
+        // {
+        //     return false;
+        // }
 
-        if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
-        {
-            var small = Shared._smallCycles;
-            if (small is not null)
-            {
-                ulong cached = small[(int)divisor];
-                if (cached != 0UL)
-                {
-                    return cached == exponent;
-                }
-            }
-        }
+        // Production scans never revisit divisors in the small-cycle snapshot, so skip the legacy lookup here. Factoring and
+        // diagnostics can query TryMatchSmallCycleForFactoring when they need those cached values.
+        // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+        // {
+        //     var small = Shared._smallCycles;
+        //     if (small is not null)
+        //     {
+        //         ulong cached = small[(int)divisor];
+        //         if (cached != 0UL)
+        //         {
+        //             return cached == exponent;
+        //         }
+        //     }
+        // }
 
         if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL)
         {
@@ -98,6 +102,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);
@@ -109,7 +114,7 @@ public class MersenneDivisorCycles
             return cycleLength == exponent;
         }
 
-        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+        return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL;
     }
 
     public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
@@ -119,35 +124,70 @@ public class MersenneDivisorCycles
             return false;
         }
 
-        if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+        // Candidate divisors on the production path sit well above the precomputed snapshot, so skip the small-cycle lookup
+        // here and leave factoring scenarios to call TryMatchSmallCycleForFactoring when necessary.
+        // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
+        // {
+        //     MersenneDivisorCycles shared = Shared;
+        //     ulong[]? small = shared._smallCycles;
+        //     if (small is not null)
+        //     {
+        //         ulong cached = small[(int)divisor];
+        //         if (cached != 0UL)
+        //         {
+        //             return cached == exponent;
+        //         }
+        //     }
+        // }
+
+        return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL;
+    }
+
+    public static bool TryMatchSmallCycleForFactoring(ulong divisor, ulong exponent, out ulong cycleLength)
+    {
+        if (divisor > PerfectNumberConstants.MaxQForDivisorCycles)
         {
-            MersenneDivisorCycles shared = Shared;
-            ulong[]? small = shared._smallCycles;
-            if (small is not null)
-            {
-                ulong cached = small[(int)divisor];
-                if (cached != 0UL)
-                {
-                    return cached == exponent;
-                }
-            }
+            // Factoring helpers probe arbitrary divisors, so anything above the snapshot limit legitimately misses here.
+            cycleLength = 0UL;
+            return false;
         }
 
-        return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL;
+        ulong[]? small = Shared._smallCycles;
+        if (small is null)
+        {
+            // Snapshot recycling during diagnostics can temporarily clear the table; surface the miss so callers recompute.
+            cycleLength = 0UL;
+            return false;
+        }
+
+        ulong cached = small[(int)divisor];
+        if (cached == 0UL)
+        {
+            // Only a subset of small divisors carry cached cycles; signal the miss so factoring falls back to the GPU helper.
+            cycleLength = 0UL;
+            return false;
+        }
+
+        cycleLength = cached;
+        return cached == exponent;
     }
 
     private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent)
     {
-        if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
-        {
-            return false;
-        }
+        // Production workloads pre-filter invalid inputs, so skip the defensive guards here. Factoring and diagnostics should
+        // validate inputs before calling into the candidate path.
+        // if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL)
+        // {
+        //     return false;
+        // }
 
         UInt128 step = (UInt128)exponent << 1;
-        if (step == UInt128.Zero)
-        {
-            return false;
-        }
+        // The overflow bail-out is handled by factoring-specific call sites; keep the historical branch commented to document the
+        // legacy behavior.
+        // if (step == UInt128.Zero)
+        // {
+        //     return false;
+        // }
 
         UInt128 adjusted = (UInt128)divisor - UInt128.One;
         return adjusted % step == UInt128.Zero;

PerfectNumbers.Core/ULongExtensions.cs

@@ -12,13 +12,17 @@ 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)
     {
-        if (q <= 2UL)
-        {
-            return 0UL;
-        }
+        // Order computations only run on admissible divisors, so the small-prime guard remains commented out to keep the
+        // path branch-free.
+        // if (q <= 2UL)
+        // {
+        //     return 0UL;
+        // }
 
         ulong order = q - 1UL, prime, temp;
         uint[] smallPrimes = PrimesGenerator.SmallPrimes;
@@ -77,7 +81,7 @@ public static class ULongExtensions
     public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift)
     {
         ulong high = value >> shift;
-		ulong mask = (1UL << shift) - 1UL;
+        ulong mask = (1UL << shift) - 1UL;
         ulong low = value & mask;
 
         UInt128 highContribution = (UInt128)high * multiplier;
@@ -288,26 +292,220 @@ public static class ULongExtensions
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    // TODO: Replace this fallback with the UInt128 Montgomery helper measured fastest in
-    // MulMod64Benchmarks so CPU callers stop paying for triple modulo operations.
-    public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus);
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus)
+    public static ulong MulMod64(this ulong a, ulong b, ulong modulus)
     {
-        // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64,
-        // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs.
-        GpuUInt128 state = new(a % modulus);
-        return state.MulMod(b, modulus);
+        // Scanner inputs always provide odd prime moduli, so the legacy modulus guard stays commented out to avoid extra
+        // branching.
+        // if (modulus <= 1UL)
+        // {
+        //     return 0UL;
+        // }
+
+        ulong reducedA = a >= modulus ? a % modulus : a;
+        ulong reducedB = b >= modulus ? b % modulus : b;
+        GpuUInt128 state = new(reducedA);
+        return state.MulMod(reducedB, modulus);
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus)
+    public static ulong Pow2ModWindowedCpu(this ulong exponent, ulong modulus)
     {
-        // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the
-        // 5-40× slowdown seen across real-world operand distributions.
-        GpuUInt128 state = new(a);
-        return state.MulModWithNativeModulo(b, modulus);
+        // 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;
+        // }
+
+        // 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 diagnostic callers
+        // still hit this small-exponent fast path, so keep the fallback in place.
+        if (exponent <= Pow2WindowFallbackThreshold)
+        {
+            return Pow2ModBinaryFallback(exponent, modulus);
+        }
+
+        int bitLength = GetPortableBitLength(exponent);
+        int windowSize = GetWindowSize(bitLength);
+        int oddPowerCount = 1 << (windowSize - 1);
+
+        Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
+            ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount]
+            : new ulong[oddPowerCount];
+        oddPowerStorage = oddPowerStorage[..oddPowerCount];
+        InitializePlainOddPowers(modulus, oddPowerStorage);
+
+        // 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)
+        {
+            if (((exponent >> index) & 1UL) == 0UL)
+            {
+                result = result.MulMod64(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.MulMod64(result, modulus);
+            }
+
+            ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL);
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+            result = result.MulMod64(oddPowerStorage[tableIndex], modulus);
+
+            index = windowStart - 1;
+        }
+
+        return result;
+    }
+
+    [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;
+        // Production moduli are always odd primes greater than two, so normalize the seed values directly.
+        ulong result = 1UL;
+        ulong baseValue = 2UL;
+
+        while (remainingExponent != 0UL)
+        {
+            if ((remainingExponent & 1UL) != 0UL)
+            {
+                result = result.MulMod64(baseValue, modulus);
+            }
+
+            remainingExponent >>= 1;
+            if (remainingExponent == 0UL)
+            {
+                break;
+            }
+
+            baseValue = baseValue.MulMod64(baseValue, modulus);
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus)
+    {
+        // The GPU pipeline feeds only odd prime moduli, so the removed guard would never fire for
+        // EvenPerfectBitScanner workloads.
+        // if (modulus <= 1UL)
+        // {
+        //     return 0UL;
+        // }
+
+        // The GPU pipeline only emits positive exponents on production workloads, so skip the zero guard to avoid an extra branch.
+        // Diagnostic harnesses must normalize their inputs before invoking this helper.
+        // if (exponent == 0UL)
+        // {
+        //     return 1UL;
+        // }
+
+        // The GPU kernels mostly process large exponents as well, but the binary fallback still serves
+        // diagnostic harnesses that 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.MulMod64(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.MulMod64(result, modulus);
+            }
+
+            ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL);
+            int tableIndex = (int)((windowValue - 1UL) >> 1);
+            result = result.MulMod64(oddPowers.Get(tableIndex), modulus);
+
+            index = windowStart - 1;
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Pow2ModBinaryGpuFallback(ulong exponent, ulong modulus)
+    {
+        ulong remainingExponent = exponent;
+        // Production moduli are always odd primes greater than two, so normalize the seed values directly.
+        ulong result = 1UL;
+        ulong baseValue = 2UL;
+
+        while (remainingExponent != 0UL)
+        {
+            if ((remainingExponent & 1UL) != 0UL)
+            {
+                result = result.MulMod64(baseValue, modulus);
+            }
+
+            remainingExponent >>= 1;
+            if (remainingExponent == 0UL)
+            {
+                break;
+            }
+
+            baseValue = baseValue.MulMod64(baseValue, modulus);
+        }
+
+        return result;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -347,10 +545,16 @@ public static class ULongExtensions
         ulong result = divisor.MontgomeryOne;
         ulong nPrime = divisor.NPrime;
 
-        Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount] : default;
+        Span<ulong> oddPowerStorage = oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount
+            ? stackalloc ulong[PerfectNumberConstants.MaxOddPowersCount]
+            : default;
         ulong[]? pooledOddPowers = null;
 
-		if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold)
+        if (oddPowerCount <= PerfectNumberConstants.MaxOddPowersCount)
+        {
+            oddPowerStorage = oddPowerStorage[..oddPowerCount];
+        }
+        else if (oddPowerCount < PerfectNumberConstants.PooledArrayThreshold)
         {
             oddPowerStorage = new ulong[oddPowerCount];
         }
@@ -548,6 +752,187 @@ public static class ULongExtensions
     }
 
     [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;
+        }
+
+        ulong baseValue = 2UL % modulus;
+        destination[0] = baseValue;
+        if (oddPowerCount == 1)
+        {
+            return;
+        }
+
+        ulong square = baseValue.MulMod64(baseValue, modulus);
+        for (int i = 1; i < oddPowerCount; i++)
+        {
+            // Reuse baseValue to accumulate successive odd powers and avoid allocating another local.
+            baseValue = baseValue.MulMod64(square, modulus);
+            destination[i] = baseValue;
+        }
+    }
+
+    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;
+
+            // Window planning always requests at least one odd power, so keep the zero-count check commented out to avoid
+            // branching here.
+            // if (_count == 0)
+            // {
+            //     return;
+            // }
+
+            ulong baseValue = 2UL;
+            _value0 = baseValue;
+            if (_count == 1)
+            {
+                return;
+            }
+
+            ulong square = baseValue.MulMod64(baseValue, modulus);
+            ulong current = baseValue;
+            for (int i = 1; i < _count; i++)
+            {
+                current = current.MulMod64(square, modulus);
+                SetValue(i, current);
+            }
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public ulong Get(int index)
+        {
+            // Windowed lookups never exceed the populated table bounds on production workloads, so leave the check commented
+            // out to remove the branch.
+            // 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;

todo_list.txt

@@ -122,7 +122,6 @@ PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:203:
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207:                // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214:                // TODO: Remove this no-cycle branch and require on-demand cycle computation (without
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240:            // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the
-PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262:                : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283:        // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356:                // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips
 PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396:            // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here
@@ -266,7 +265,6 @@ todo_list.txt:121:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTeste
 todo_list.txt:122:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207:                // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device
 todo_list.txt:123:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214:                // TODO: Remove this no-cycle branch and require on-demand cycle computation (without
 todo_list.txt:124:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240:            // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the
-todo_list.txt:125:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262:                : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the
 todo_list.txt:126:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283:        // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI
 todo_list.txt:127:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356:                // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips
 todo_list.txt:128:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396:            // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here
@@ -407,17 +405,6 @@ PerfectNumbers.Core\ModResidueTracker.cs:257:                        // TODO: Us
 PerfectNumbers.Core\ModResidueTracker.cs:263:                // TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths
 PerfectNumbers.Core\ModResidueTracker.cs:295:                                // TODO: Route this through the planned UInt128 windowed powmod so the inner MulMod adopts the
 PerfectNumbers.Core\ModResidueTracker.cs:314:                // TODO: Replace this double-and-add fallback with the UInt128 intrinsic-backed multiplier once the
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:49:        set => _useDivisorCycles = value; // TODO: Remove this flag once GPU runs unconditionally fetch cycle lengths so every divisor check benefits from the measured speedups.
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:298:                // TODO: Stage MontgomeryDivisorData in a GPU-ready cache keyed by divisor so we can copy
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:372:        // TODO: Replace the block lease with a direct single-cycle lookup so the GPU path
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:395:                        // TODO: Drop this block hand-off once the cache exposes a direct
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:409:                            // TODO: Replace this `%` with the cycle-cache helper that subtracts precomputed
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:425:                        // TODO: Trigger an immediate GPU computation for the single missing cycle without storing
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:476:                    // TODO: Preload ProcessEightBitWindows-friendly divisor metadata for GPU batches so this
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:563:        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.
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:632:                // TODO: Replace this guard with an on-demand cycle computation that runs on the device selected
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:787:                        // TODO: Replace this `%` with the divisor-cycle remainder helper so we reuse cached
-PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:858:        results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); // TODO: Swap this Montgomery-domain pow2 with the windowed helper slated to replace the slower ladder in Pow2Minus1Mod benchmarks.
 PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:43:                // TODO: Swap this per-run enumerator creation for the shared PrimeIterator so by-divisor scans reuse the
 PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:219:                        // TODO: Swap Task.Factory.StartNew for the shared low-overhead work queue once the scanner's
 PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:240:                                bool useDivisorCycles = tester.UseDivisorCycles; // TODO: Remove the conditional path once all testers always enable divisor cycles to keep every divisor scan on the faster cached-length track.