From 6bd8945565d8d403af83a81bca6116bca4e715c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 9 Oct 2025 04:31:03 +0200 Subject: [PATCH] Use analytic delta updates in CPU divisor sessions --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 240 +++++++++++++++++- docs/bydivisor_cpu_path.md | 119 +++++++++ 2 files changed, 348 insertions(+), 11 deletions(-) create mode 100644 docs/bydivisor_cpu_path.md diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 194a642..0d2c874 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/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 _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 primes, Span 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. - // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors. - var exponentStepper = new ExponentRemainderStepper(cachedData); - if (!exponentStepper.IsValidModulus) + var cycleStepper = new CycleRemainderStepper(divisorCycle); + ExponentRemainderStepper exponentStepper = default; + bool exponentStepperInitialized = false; + byte hitValue = 0; // Reused for every prime processed in this batch. + + bool EnsureExponentStepper() { - hits.Clear(); - return; + if (exponentStepperInitialized) + { + 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 - ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0) - : (byte)0; + if (remainder == 0UL) + { + 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) diff --git a/docs/bydivisor_cpu_path.md b/docs/bydivisor_cpu_path.md new file mode 100644 index 0000000..4650e4d --- /dev/null +++ b/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.