Merge pull request #227 from mlemanczyk/codex/refactor-exponentremainderstepper-structure

Clarify Pollard-Rho guard and seed CPU stepping upfront
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-10-23 20:15:59 +02:00
commit 6af0ae900e
4 files changed
+160 -74

No files matched your search

@@ -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++)
+93 -59
View File
@@ -1,5 +1,9 @@
namespace PerfectNumbers.Core;
/// <summary>
/// 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.
/// </summary>
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;
}
/// <summary>
/// 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.
/// </summary>
public void Reset()
{
_previousExponent = 0UL;
_currentMontgomery = _montgomeryOne;
_hasState = false;
}
public ulong ComputeNext(ulong exponent)
/// <summary>
/// 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.
/// </summary>
/// <param name="exponent">First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata.</param>
/// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
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();
}
/// <summary>
/// Initializes the CPU stepping state and reports whether the seed exponent already yields unity.
/// </summary>
/// <param name="exponent">First exponent evaluated for the current divisor. Must match the divisor passed to the constructor.</param>
/// <returns><see langword="true"/> when the seed exponent produces a Montgomery residue equal to one.</returns>
public bool InitializeCpuIsUnity(ulong exponent)
{
InitializeCpuState(exponent);
return _currentMontgomery == _montgomeryOne;
}
/// <summary>
/// Advances the CPU stepping state to the next exponent in a strictly ascending sequence that shares the same divisor.
/// Callers must invoke <see cref="InitializeCpu(ulong)"/> beforehand.
/// </summary>
/// <param name="exponent">Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent.</param>
/// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
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();
}
/// <summary>
/// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence.
/// </summary>
/// <param name="exponent">First exponent evaluated on the GPU path. The same divisor metadata must be reused for future calls.</param>
/// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
public ulong InitializeGpu(ulong exponent)
{
InitializeGpuState(exponent);
return ReduceCurrent();
}
/// <summary>
/// Initializes the GPU stepping state and reports whether the seed exponent produces unity.
/// </summary>
/// <param name="exponent">First exponent evaluated on the GPU path. Must match the divisor supplied during construction.</param>
/// <returns><see langword="true"/> when the seed exponent produces a Montgomery residue equal to one.</returns>
public bool InitializeGpuIsUnity(ulong exponent)
{
InitializeGpuState(exponent);
return _currentMontgomery == _montgomeryOne;
}
/// <summary>
/// Advances the CPU stepping state to the next exponent and reports whether it evaluates to unity.
/// Callers must seed the stepper through <see cref="InitializeCpu(ulong)"/> or <see cref="InitializeCpuIsUnity(ulong)"/> and then
/// pass exponents in strictly ascending order.
/// </summary>
/// <param name="exponent">Next exponent in the ascending sequence.</param>
/// <returns><see langword="true"/> when the canonical residue equals one.</returns>
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;
}
/// <summary>
/// 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.
/// </summary>
/// <param name="exponent">Exponent associated with <paramref name="montgomeryResult"/>.</param>
/// <param name="montgomeryResult">Residue in Montgomery form produced externally for the configured divisor.</param>
/// <param name="isUnity">Outputs whether the supplied residue equals Montgomery one.</param>
/// <returns><see langword="true"/> when the residue is accepted and cached.</returns>
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;
}
/// <summary>
/// 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.
/// </summary>
/// <param name="exponent">Next exponent in the ascending sequence.</param>
/// <param name="montgomeryDelta">Montgomery residue delta corresponding to the exponent gap.</param>
/// <param name="isUnity">Outputs whether the resulting residue equals Montgomery one.</param>
/// <returns><see langword="true"/> when the advance succeeds.</returns>
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);
}
+9 -5
View File
@@ -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)
{
@@ -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<ulong> 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)
{