From 7feeb089437785c2f8c7afb829f2f4e5017e58e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 23 Oct 2025 19:50:39 +0200 Subject: [PATCH] Clarify Pollard-Rho parity guard and seed CPU batches --- .../Cpu/MersenneCpuDivisorScanSession.cs | 7 +- .../ExponentRemainderStepper.cs | 152 +++++++++++------- PerfectNumbers.Core/MersenneDivisorCycles.cs | 14 +- .../PrimeOrderCalculator.Cpu.cs | 61 ++++++- 4 files changed, 160 insertions(+), 74 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs index 1279369..10a92a9 100644 --- a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs +++ b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs @@ -43,16 +43,13 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi // 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) - { - throw new InvalidOperationException($"Invalid modulus {cachedData.Modulus} for divisor scans."); - } var cycleStepper = new CycleRemainderStepper(divisorCycle); + bool initialUnity = exponentStepper.InitializeCpuIsUnity(primes[0]); ulong remainder = cycleStepper.Initialize(primes[0]); hits[0] = remainder == 0UL - ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0) + ? (initialUnity ? (byte)1 : (byte)0) : (byte)0; for (int i = 1; i < length; i++) diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepper.cs index 5cbef65..0722880 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepper.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepper.cs @@ -1,5 +1,9 @@ namespace PerfectNumbers.Core; +/// +/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates. +/// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order. +/// internal struct ExponentRemainderStepper { private readonly MontgomeryDivisorData _divisor; @@ -8,7 +12,6 @@ internal struct ExponentRemainderStepper private readonly ulong _montgomeryOne; private ulong _previousExponent; private ulong _currentMontgomery; - private bool _hasState; public ExponentRemainderStepper(in MontgomeryDivisorData divisor) { @@ -18,13 +21,8 @@ internal struct ExponentRemainderStepper _montgomeryOne = divisor.MontgomeryOne; _previousExponent = 0UL; _currentMontgomery = divisor.MontgomeryOne; - _hasState = false; } - public bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL; - - public bool HasState => _hasState; - public ulong PreviousExponent => _previousExponent; public bool MatchesDivisor(in MontgomeryDivisorData divisor) @@ -32,29 +30,47 @@ internal struct ExponentRemainderStepper return _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne; } + /// + /// Clears the cached Montgomery residue so the stepper can be reused for another divisor. + /// Callers must reinitialize the stepper before consuming any residues after a reset. + /// public void Reset() { _previousExponent = 0UL; _currentMontgomery = _montgomeryOne; - _hasState = false; } - public ulong ComputeNext(ulong exponent) + /// + /// Initializes the stepper for CPU calculations using the first exponent in a strictly ascending sequence. + /// The returned value is the canonical (non-Montgomery) residue corresponding to the provided exponent. + /// + /// First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata. + /// The canonical residue produced by 2^exponent mod divisor. + public ulong InitializeCpu(ulong exponent) { - if (!IsValidModulus) - { - _hasState = false; - return 0UL; - } + InitializeCpuState(exponent); + return ReduceCurrent(); + } - if (!_hasState || exponent <= _previousExponent) - { - _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); - _previousExponent = exponent; - _hasState = true; - return ReduceCurrent(); - } + /// + /// Initializes the CPU stepping state and reports whether the seed exponent already yields unity. + /// + /// First exponent evaluated for the current divisor. Must match the divisor passed to the constructor. + /// when the seed exponent produces a Montgomery residue equal to one. + public bool InitializeCpuIsUnity(ulong exponent) + { + InitializeCpuState(exponent); + return _currentMontgomery == _montgomeryOne; + } + /// + /// Advances the CPU stepping state to the next exponent in a strictly ascending sequence that shares the same divisor. + /// Callers must invoke beforehand. + /// + /// Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent. + /// The canonical residue produced by 2^exponent mod divisor. + public ulong ComputeNextCpu(ulong exponent) + { ulong delta = exponent - _previousExponent; // TODO: Once divisor cycle lengths are mandatory, pull the delta multiplier from the // single-block divisor-cycle snapshot so we can skip the powmod entirely and reuse the @@ -65,23 +81,37 @@ internal struct ExponentRemainderStepper return ReduceCurrent(); } + /// + /// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence. + /// + /// First exponent evaluated on the GPU path. The same divisor metadata must be reused for future calls. + /// The canonical residue produced by 2^exponent mod divisor. + public ulong InitializeGpu(ulong exponent) + { + InitializeGpuState(exponent); + return ReduceCurrent(); + } + + /// + /// Initializes the GPU stepping state and reports whether the seed exponent produces unity. + /// + /// First exponent evaluated on the GPU path. Must match the divisor supplied during construction. + /// when the seed exponent produces a Montgomery residue equal to one. + public bool InitializeGpuIsUnity(ulong exponent) + { + InitializeGpuState(exponent); + return _currentMontgomery == _montgomeryOne; + } + + /// + /// Advances the CPU stepping state to the next exponent and reports whether it evaluates to unity. + /// Callers must seed the stepper through or and then + /// pass exponents in strictly ascending order. + /// + /// Next exponent in the ascending sequence. + /// when the canonical residue equals one. public bool ComputeNextIsUnity(ulong exponent) { - if (!IsValidModulus) - { - throw new InvalidOperationException($"Modulus is invalid for exponent {exponent}"); - // _hasState = false; - // return false; - } - - if (!_hasState || exponent <= _previousExponent) - { - _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); - _previousExponent = exponent; - _hasState = true; - return _currentMontgomery == _montgomeryOne; - } - ulong delta = exponent - _previousExponent; // TODO: Reuse the divisor-cycle derived Montgomery delta once the cache exposes single-cycle // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor. @@ -91,40 +121,32 @@ internal struct ExponentRemainderStepper return _currentMontgomery == _montgomeryOne; } + /// + /// Seeds the stepper with a precomputed Montgomery residue produced by the caller. + /// The provided exponent becomes the new baseline and all subsequent calls must continue with greater exponents. + /// + /// Exponent associated with . + /// Residue in Montgomery form produced externally for the configured divisor. + /// Outputs whether the supplied residue equals Montgomery one. + /// when the residue is accepted and cached. public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity) { - if (!IsValidModulus) - { - _hasState = false; - isUnity = false; - return false; - } - _currentMontgomery = montgomeryResult; _previousExponent = exponent; - _hasState = true; isUnity = _currentMontgomery == _montgomeryOne; return true; } + /// + /// Advances the stepper using a caller-supplied Montgomery delta that represents the gap to the next exponent. + /// Callers must seed the stepper first, reuse the same divisor metadata, and continue feeding exponents in ascending order so the delta remains valid. + /// + /// Next exponent in the ascending sequence. + /// Montgomery residue delta corresponding to the exponent gap. + /// Outputs whether the resulting residue equals Montgomery one. + /// when the advance succeeds. public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity) { - if (!IsValidModulus) - { - _hasState = false; - isUnity = false; - return false; - } - - if (!_hasState || exponent <= _previousExponent) - { - _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); - _previousExponent = exponent; - _hasState = true; - isUnity = _currentMontgomery == _montgomeryOne; - return true; - } - // TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead // of relying on the caller-provided delta so incremental scans remain in sync with the snapshot // without computing additional cycles or mutating cache state. @@ -134,5 +156,17 @@ internal struct ExponentRemainderStepper return true; } + private void InitializeCpuState(ulong exponent) + { + _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); + _previousExponent = exponent; + } + + private void InitializeGpuState(ulong exponent) + { + _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); + _previousExponent = exponent; + } + private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime); } diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index b0e677e..c3fa2b1 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -564,7 +564,8 @@ public class MersenneDivisorCycles stepper.Reset(); - for (int j = 0; j < actual; j++) + evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]); + for (int j = 1; j < actual; j++) { evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); } @@ -582,10 +583,13 @@ public class MersenneDivisorCycles private static ulong PollardRho64(ulong n) { - if ((n & 1UL) == 0UL) - { - return 2UL; - } + // The EvenPerfectBitScanner path removes all factors of two before invoking Pollard-Rho, so this branch never executes + // outside synthetic tests and benchmarks. Keep the guard documented to prevent future regressions without paying the hot- + // path cost. + // if ((n & 1UL) == 0UL) + // { + // return 2UL; + // } while (true) { diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 3545015..5372881 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -165,7 +165,14 @@ internal static partial class PrimeOrderCalculator ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData); int candidateCount = candidates.Length; - for (int i = 0; i < candidateCount; i++) + + if (stepper.InitializeCpuIsUnity(candidates[0])) + { + ThreadStaticPools.ReturnExponentStepper(stepper); + return false; + } + + for (int i = 1; i < candidateCount; i++) { if (stepper.ComputeNextIsUnity(candidates[i])) { @@ -300,7 +307,8 @@ internal static partial class PrimeOrderCalculator stepper.Reset(); - for (int j = 0; j < actual; j++) + evaluations[0] = stepper.InitializeCpuIsUnity(candidates[0]); + for (int j = 1; j < actual; j++) { evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); } @@ -459,7 +467,12 @@ internal static partial class PrimeOrderCalculator stepper.Reset(); - for (int j = 0; j < actual; j++) + if (stepper.InitializeCpuIsUnity(candidates[0])) + { + return true; + } + + for (int j = 1; j < actual; j++) { if (stepper.ComputeNextIsUnity(candidates[j])) { @@ -535,6 +548,7 @@ internal static partial class PrimeOrderCalculator bool allowGpuBatch = true; Span candidateSpan = CollectionsMarshal.AsSpan(candidates); ExponentRemainderStepper powStepper = ThreadStaticPools.RentExponentStepper(divisorData); + bool powStepperInitialized = false; // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); int index = 0; @@ -594,7 +608,33 @@ internal static partial class PrimeOrderCalculator } } - for (int i = 0; i < batchSize && powUsed < powBudget; i++) + int startIndex = 0; + + if (!gpuSuccess && !powStepperInitialized) + { + ulong candidate = batch[0]; + powUsed++; + + bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); + powStepperInitialized = true; + if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepper(powStepper); + result = candidate; + return true; + } + + startIndex = 1; + } + + for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) { ulong candidate = batch[i]; powUsed++; @@ -936,7 +976,18 @@ internal static partial class PrimeOrderCalculator stepper.Reset(); - for (int j = 0; j < actual; j++) + if (powUsed >= powBudget && powBudget > 0) + { + return true; + } + + powUsed++; + if (stepper.InitializeCpuIsUnity(candidates[0])) + { + return true; + } + + for (int j = 1; j < actual; j++) { if (powUsed >= powBudget && powBudget > 0) {