Merge pull request #133 from mlemanczyk/codex/analyze-evenperfectbitscanner-performance-with-mersenne

Implement analytic delta stepping for CPU divisor scans

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -343,6 +343,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
     private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
     {
         private readonly MersenneNumberDivisorByDivisorCpuTester _owner;
+        private readonly Dictionary<ulong, PrimeDivisorState> _primeStates = new();
         private bool _disposed;
 
         internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner)
@@ -353,6 +354,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         internal void Reset()
         {
             _disposed = false;
+            _primeStates.Clear();
         }
 
         public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> hits)
@@ -384,21 +386,50 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 }
             }
 
-            // Keep these remainder steppers in place so future updates continue reusing the previously computed residues.
+            var cycleStepper = new CycleRemainderStepper(divisorCycle);
-            // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors.
+            ExponentRemainderStepper exponentStepper = default;
-            var exponentStepper = new ExponentRemainderStepper(cachedData);
+            bool exponentStepperInitialized = false;
-            if (!exponentStepper.IsValidModulus)
+            byte hitValue = 0; // Reused for every prime processed in this batch.
+
+            bool EnsureExponentStepper()
             {
-                hits.Clear();
+                if (exponentStepperInitialized)
-                return;
+                {
+                    return true;
+                }
+
+                exponentStepper = new ExponentRemainderStepper(cachedData);
+                if (!exponentStepper.IsValidModulus)
+                {
+                    hits.Clear();
+                    return false;
+                }
+
+                exponentStepperInitialized = true;
+                return true;
             }
 
-            var cycleStepper = new CycleRemainderStepper(divisorCycle);
-
             ulong remainder = cycleStepper.Initialize(primes[0]);
-            hits[0] = remainder == 0UL
+            if (remainder == 0UL)
-                ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0)
+            {
-                : (byte)0;
+                if (TryComputeAnalyticHit(primes[0], divisor, cachedData, out hitValue))
+                {
+                    hits[0] = hitValue;
+                }
+                else
+                {
+                    if (!EnsureExponentStepper())
+                    {
+                        return;
+                    }
+
+                    hits[0] = exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0;
+                }
+            }
+            else
+            {
+                hits[0] = 0;
+            }
 
             for (int i = 1; i < length; i++)
             {
@@ -409,10 +440,197 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                     continue;
                 }
 
+                if (TryComputeAnalyticHit(primes[i], divisor, cachedData, out hitValue))
+                {
+                    hits[i] = hitValue;
+                    continue;
+                }
+
+                if (!EnsureExponentStepper())
+                {
+                    return;
+                }
+
                 hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0;
             }
         }
 
+        private bool TryComputeAnalyticHit(ulong prime, ulong divisor, in MontgomeryDivisorData divisorData, out byte hit)
+        {
+            hit = 0;
+
+            if (!TryGetMersenneMinusOne(prime, out ulong mersenne))
+            {
+                return false;
+            }
+
+            if (divisor <= 1UL)
+            {
+                return false;
+            }
+
+            ulong step = prime << 1;
+            if (step == 0UL)
+            {
+                return false;
+            }
+
+            if (!_primeStates.TryGetValue(prime, out PrimeDivisorState state) || !state.HasState)
+            {
+                state = CreatePrimeState(step, divisor, mersenne, divisorData);
+                _primeStates[prime] = state;
+                hit = state.Remainder == 0UL ? (byte)1 : (byte)0;
+                return true;
+            }
+
+            if (state.Step != step || divisor <= state.LastDivisor)
+            {
+                state = CreatePrimeState(step, divisor, mersenne, divisorData);
+                _primeStates[prime] = state;
+                hit = state.Remainder == 0UL ? (byte)1 : (byte)0;
+                return true;
+            }
+
+            ulong difference = divisor - state.LastDivisor;
+            if (difference == 0UL)
+            {
+                hit = state.Remainder == 0UL ? (byte)1 : (byte)0;
+                return true;
+            }
+
+            if (difference % step != 0UL)
+            {
+                state = CreatePrimeState(step, divisor, mersenne, divisorData);
+                _primeStates[prime] = state;
+                hit = state.Remainder == 0UL ? (byte)1 : (byte)0;
+                return true;
+            }
+
+            ulong increments = difference / step;
+            ulong quotient = state.Quotient;
+            ulong remainder = state.Remainder;
+            ulong currentDivisor = state.LastDivisor;
+
+            for (ulong index = 0UL; index < increments; index++)
+            {
+                ulong nextDivisor = currentDivisor + step;
+                Int128 stepTimesQuotient = (Int128)step * (Int128)quotient;
+                Int128 numerator = stepTimesQuotient - (Int128)remainder;
+                Int128 divisorInt = (Int128)nextDivisor;
+                Int128 deltaInt = numerator >= 0
+                    ? (numerator + divisorInt - Int128.One) / divisorInt
+                    : numerator / divisorInt;
+
+                if (deltaInt < Int128.Zero)
+                {
+                    deltaInt = Int128.Zero;
+                }
+
+                ulong delta = (ulong)deltaInt;
+                if (delta > quotient)
+                {
+                    state = CreatePrimeState(step, divisor, mersenne, divisorData);
+                    _primeStates[prime] = state;
+                    hit = state.Remainder == 0UL ? (byte)1 : (byte)0;
+                    return true;
+                }
+
+                Int128 t = (Int128)remainder - stepTimesQuotient;
+                Int128 nextRemainder = t + deltaInt * divisorInt;
+                if (nextRemainder < Int128.Zero || nextRemainder >= divisorInt)
+                {
+                    state = CreatePrimeState(step, divisor, mersenne, divisorData);
+                    _primeStates[prime] = state;
+                    hit = state.Remainder == 0UL ? (byte)1 : (byte)0;
+                    return true;
+                }
+
+                quotient -= delta;
+                remainder = (ulong)nextRemainder;
+                currentDivisor = nextDivisor;
+            }
+
+            ulong residue = ComputeResidueFromRemainder(remainder, divisor);
+            state.LastDivisor = divisor;
+            state.Quotient = quotient;
+            state.Remainder = remainder;
+            state.MontgomeryResidue = ToMontgomeryResidue(residue, divisorData);
+            // Reuse state to persist the updated quotient and remainder tuple for the next divisor.
+            _primeStates[prime] = state;
+
+            hit = remainder == 0UL ? (byte)1 : (byte)0;
+            return true;
+        }
+
+        private static PrimeDivisorState CreatePrimeState(ulong step, ulong divisor, ulong mersenne, in MontgomeryDivisorData divisorData)
+        {
+            PrimeDivisorState state = new()
+            {
+                Step = step,
+                LastDivisor = divisor,
+                HasState = true,
+            };
+
+            ulong quotient = Math.DivRem(mersenne, divisor, out ulong remainder);
+            state.Quotient = quotient;
+            state.Remainder = remainder;
+            ulong residue = ComputeResidueFromRemainder(remainder, divisor);
+            state.MontgomeryResidue = ToMontgomeryResidue(residue, divisorData);
+            return state;
+        }
+
+        private static bool TryGetMersenneMinusOne(ulong exponent, out ulong mersenne)
+        {
+            if (exponent <= 63UL)
+            {
+                mersenne = (1UL << (int)exponent) - 1UL;
+                return true;
+            }
+
+            if (exponent == 64UL)
+            {
+                mersenne = ulong.MaxValue;
+                return true;
+            }
+
+            mersenne = 0UL;
+            return false;
+        }
+
+        private static ulong ComputeResidueFromRemainder(ulong remainder, ulong divisor)
+        {
+            ulong residue = remainder + 1UL;
+            if (residue >= divisor)
+            {
+                residue -= divisor;
+            }
+
+            return residue;
+        }
+
+        private static ulong ToMontgomeryResidue(ulong value, in MontgomeryDivisorData divisorData)
+        {
+            ulong modulus = divisorData.Modulus;
+            if (modulus <= 1UL)
+            {
+                return 0UL;
+            }
+
+            value %= modulus;
+            UInt128 shifted = (UInt128)value << 64;
+            return (ulong)(shifted % modulus);
+        }
+
+        private struct PrimeDivisorState
+        {
+            public ulong Step;
+            public ulong LastDivisor;
+            public ulong Quotient;
+            public ulong Remainder;
+            public ulong MontgomeryResidue;
+            public bool HasState;
+        }
+
         public void Dispose()
         {
             if (_disposed)

docs/bydivisor_cpu_path.md

@@ -0,0 +1,119 @@
+# CPU `--mersenne=bydivisor` pipeline notes
+
+This document summarizes how the `EvenPerfectBitScanner` executable executes the
+`--mersenne=bydivisor` mode on the CPU, and how the analytic observations we
+derived for the quotients
+
+\[
+Q_k = \frac{N}{d_k}, \qquad d_k = 2pk + 1, \qquad N = 2^p - 1,
+\]
+
+can help us identify concrete optimization opportunities.
+
+## High-level control flow
+
+* `Program.Main` wires `_byDivisorTester` to a CPU implementation whenever the
+  CLI switches resolve to `--mersenne=bydivisor` and
+  `--mersenne-device=cpu`. The CLI also forces a block size of 1 so that the
+  pipeline pushes exponents to the tester one at a time.
+  (`Program.cs`, lines 459-520.)
+* When the mode starts, `RunByDivisorMode` delegates to
+  `MersenneNumberDivisorByDivisorTester.Run`, passing the candidate list from
+  `--filter-p`, the CPU tester instance and the optional replay of previous
+  outcomes. (`Program.cs`, lines 809-826.)
+* `MersenneNumberDivisorByDivisorTester.Run` sorts and deduplicates the
+  candidate exponents, drops the ones already confirmed in prior runs, validates
+  primality with the `Open.Numeric.Primes` enumerator and configures the tester
+  with the largest remaining exponent. (`MersenneNumberDivisorByDivisorTester.cs`,
+  lines 14-162.)
+* For each confirmed prime exponent, the CPU tester computes the largest divisor
+  it needs to probe through `ComputeAllowedMaxDivisor`. If the range is non-empty
+  it enters the `CheckDivisors` loop. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
+  lines 83-110 and 161-266.)
+
+## Inner loop behaviour (`CheckDivisors`)
+
+The hot path is the `while (divisor <= limit)` loop in
+`MersenneNumberDivisorByDivisorCpuTester.CheckDivisors`.
+
+1. The loop enumerates divisor candidates of the form `d_k = 2pk + 1` by
+   incrementing `divisor` by `step = 2p` every iteration. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
+   lines 177-256.)
+2. Small modulus filters track `d_k (mod 3·5·7·8·10·11)` to discard impossible
+   candidates without touching the divisor-cycle cache. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
+   lines 195-263.)
+3. For every admissible divisor the tester fetches its Montgomery data and the
+   cached cycle length. When the cache does not contain an entry the
+   implementation recomputes the cycle by factoring the divisor. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
+   lines 223-254.)
+4. A cycle equal to the tested exponent proves divisibility, so the loop aborts
+   as soon as one hit is found. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
+   lines 241-247.)
+
+`DivisorScanSession.CheckDivisor` follows the same `d_k` progression when a
+batch API is used. It keeps a running Montgomery residue via
+`ExponentRemainderStepper`, but the stepper currently recomputes the delta by
+calling a full `Pow2MontgomeryModWindowed` whenever the exponent jumps.
+(`MersenneNumberDivisorByDivisorCpuTester.cs`, lines 358-405 and
+`ExponentRemainderStepper.cs`, lines 26-95.)
+
+## Where the analytic model fits
+
+The formulas we established for the quotient sequence
+
+\[
+\frac{Q_{k+1}}{Q_k} = \frac{2pk + 1}{2p(k+1) + 1}, \qquad
+Q_k - Q_{k+1} = \frac{N·2p}{(2pk + 1)(2p(k+1) + 1)}
+\]
+
+tell us that the successive quotients and remainders shrink smoothly and that the
+per-step update depends only on the previous `(Q_k, r_k)` state and the next
+`d_{k+1}`. This lines up with several TODO comments that highlight missing
+incremental updates:
+
+* `ExponentRemainderStepper` already caches the previous exponent and a running
+  Montgomery residue, but it falls back to a full modular exponentiation to
+  compute the delta between consecutive divisors. (`ExponentRemainderStepper.cs`,
+  lines 49-72.)
+* `CheckDivisors` rebuilds the Montgomery data for every hit candidate even
+  though the sequence of `d_k` values is strictly monotone.
+
+The quotient identity gives us a way to compute `Δ_k = Q_k - Q_{k+1}` by reading
+`Q_k` and the inexpensive division `(2pQ_k - r_k) / d_{k+1}`. Once we know `Δ_k`
+we can update
+
+\[
+Q_{k+1} = Q_k - Δ_k, \qquad r_{k+1} = r_k - 2p·Q_k + Δ_k · d_{k+1},
+\]
+
+without ever evaluating `2^p mod d_{k+1}` from scratch. Incorporating this into
+`ExponentRemainderStepper` would let the CPU path reuse a single Montgomery state
+across the entire scan and swap the current delta exponentiation (which costs a
+full modular ladder) for one multiply-high with a precomputed reciprocal of
+`d_{k+1}`. The loop already maintains `d_k` and `step`, so the extra bookkeeping
+fits naturally.
+
+## Optimization ideas grounded in the code
+
+1. **Delta-driven Montgomery updates.** Extend `DivisorScanSession` so it tracks
+   `(Q_k, r_k)` alongside the Montgomery residue and uses the analytic delta
+   formula to derive the next residue. This replaces the per-step powmod in
+   `ExponentRemainderStepper.ComputeNext` with a handful of 64-bit operations.
+   (`MersenneNumberDivisorByDivisorCpuTester.cs`, lines 358-405;
+   `ExponentRemainderStepper.cs`, lines 49-72.)
+2. **Bit-range skipping.** Because `Q_{k+1} / Q_k` is exactly the rational
+   `(2pk + 1)/(2p(k+1) + 1)`, we can predict when `⌊log₂ Q_k⌋` loses a bit and
+   jump the loop counter to the first `k` that crosses a target quotient bound
+   instead of checking each divisor. This is compatible with the existing
+   `allowedMax` guard because it only modifies how we arrive at the stopping
+   divisor. (`MersenneNumberDivisorByDivisorCpuTester.cs`, lines 161-266.)
+3. **Reuse Montgomery metadata.** The loop already calls `MontgomeryDivisorDataCache.Get`
+   for each divisor, but the cache lookup is repeated even though the next
+   divisor is a small increment away. Persisting the previous `MontgomeryDivisorData`
+   together with the `(Q_k, r_k)` tuple lets us skip redundant cache reads when
+   the divisor stays within the hot range. (`MersenneNumberDivisorByDivisorCpuTester.cs`,
+   lines 223-256.)
+
+These adjustments directly leverage the quotient progression we analysed and
+align with the existing TODO list, giving us a clear path to reduce the number of
+modular exponentiations and cache lookups in the CPU by-divisor path.