diff --git a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
index df66c67..dc22238 100644
--- a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
+++ b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
@@ -10,7 +10,7 @@ internal struct ExponentRemainderStepperCpu
private readonly ulong _modulus;
private readonly ulong _nPrime;
private readonly ulong _montgomeryOne;
- private ulong _previousExponent;
+ public ulong PreviousExponent;
private ulong _currentMontgomery;
public ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor)
@@ -19,12 +19,10 @@ internal struct ExponentRemainderStepperCpu
_modulus = divisor.Modulus;
_nPrime = divisor.NPrime;
_montgomeryOne = divisor.MontgomeryOne;
- _previousExponent = 0UL;
+ PreviousExponent = 0UL;
_currentMontgomery = divisor.MontgomeryOne;
}
- public ulong PreviousExponent => _previousExponent;
-
public bool MatchesDivisor(in MontgomeryDivisorData divisor)
{
return _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne;
@@ -36,7 +34,7 @@ internal struct ExponentRemainderStepperCpu
///
public void Reset()
{
- _previousExponent = 0UL;
+ PreviousExponent = 0UL;
_currentMontgomery = _montgomeryOne;
}
@@ -71,13 +69,13 @@ internal struct ExponentRemainderStepperCpu
/// The canonical residue produced by 2^exponent mod divisor.
public ulong ComputeNextCpu(ulong exponent)
{
- ulong delta = exponent - _previousExponent;
+ 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
// cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
- _previousExponent = exponent;
+ PreviousExponent = exponent;
return ReduceCurrent();
}
@@ -90,12 +88,12 @@ internal struct ExponentRemainderStepperCpu
/// when the canonical residue equals one.
public bool ComputeNextIsUnity(ulong exponent)
{
- ulong delta = exponent - _previousExponent;
+ 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.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
- _previousExponent = exponent;
+ PreviousExponent = exponent;
return _currentMontgomery == _montgomeryOne;
}
@@ -110,7 +108,7 @@ internal struct ExponentRemainderStepperCpu
public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity)
{
_currentMontgomery = montgomeryResult;
- _previousExponent = exponent;
+ PreviousExponent = exponent;
isUnity = _currentMontgomery == _montgomeryOne;
return true;
}
@@ -129,7 +127,7 @@ internal struct ExponentRemainderStepperCpu
// of relying on the caller-provided delta so incremental scans remain in sync with the snapshot
// without computing additional cycles or mutating cache state.
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime);
- _previousExponent = exponent;
+ PreviousExponent = exponent;
isUnity = _currentMontgomery == _montgomeryOne;
return true;
}
@@ -137,7 +135,7 @@ internal struct ExponentRemainderStepperCpu
private void InitializeCpuState(ulong exponent)
{
_currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
- _previousExponent = exponent;
+ PreviousExponent = exponent;
}
private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime);
diff --git a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs
index 0d74bc5..c4c979f 100644
--- a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs
+++ b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs
@@ -13,23 +13,17 @@ internal struct ExponentRemainderStepperGpu
private const int WindowSizeMax = 8;
private readonly ulong _modulus;
- private ulong _previousExponent;
+ public ulong PreviousExponent;
private ulong _currentResidue;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public ExponentRemainderStepperGpu(in MontgomeryDivisorData divisor)
{
_modulus = divisor.Modulus;
- _previousExponent = 0UL;
+ PreviousExponent = 0UL;
_currentResidue = 1UL;
}
- public ulong PreviousExponent
- {
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- get => _previousExponent;
- }
-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public bool MatchesDivisor(in MontgomeryDivisorData divisor)
{
@@ -43,7 +37,7 @@ internal struct ExponentRemainderStepperGpu
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void Reset()
{
- _previousExponent = 0UL;
+ PreviousExponent = 0UL;
_currentResidue = 1UL;
}
@@ -103,26 +97,26 @@ internal struct ExponentRemainderStepperGpu
private void InitializeStateGpu(ulong exponent)
{
_currentResidue = Pow2ModWindowed(exponent, _modulus);
- _previousExponent = exponent;
+ PreviousExponent = exponent;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void AdvanceStateGpu(ulong exponent)
{
- ulong delta = exponent - _previousExponent;
+ ulong delta = exponent - PreviousExponent;
ulong multiplier = Pow2ModWindowed(delta, _modulus);
- _currentResidue = _currentResidue.MulMod64Gpu(multiplier, _modulus);
- _previousExponent = exponent;
+ _currentResidue = MultiplyMod(_currentResidue, multiplier, _modulus);
+ PreviousExponent = exponent;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static ulong Pow2ModWindowed(ulong exponent, ulong modulus)
{
- if (exponent == 0UL)
- {
- // The GPU by-divisor flow never routes modulus ≤ 1 here because divisors follow q = 2kp + 1 ≥ 3.
- return 1UL;
- }
+ // EvenPerfectBitScanner never provides zero exponents on the GPU path, so the guard stays commented out.
+ // if (exponent == 0UL)
+ // {
+ // return 1UL;
+ // }
int bitLength = GetBitLength(exponent);
int windowSize = GetWindowSize(bitLength);
@@ -132,7 +126,7 @@ internal struct ExponentRemainderStepperGpu
while (index >= 0)
{
ulong currentBit = (exponent >> index) & 1UL;
- ulong squared = result.MulMod64Gpu(result, modulus);
+ ulong squared = MultiplyMod(result, result, modulus);
bool processWindow = currentBit != 0UL;
result = processWindow ? result : squared;
@@ -147,7 +141,7 @@ internal struct ExponentRemainderStepperGpu
for (int square = 0; square < windowLength; square++)
{
- result = result.MulMod64Gpu(result, modulus);
+ result = MultiplyMod(result, result, modulus);
}
if (!processWindow)
@@ -158,13 +152,21 @@ internal struct ExponentRemainderStepperGpu
ulong mask = (1UL << windowLength) - 1UL;
ulong windowValue = (exponent >> windowStart) & mask;
ulong multiplier = ComputeWindowedOddPower(windowValue, modulus);
- result = result.MulMod64Gpu(multiplier, modulus);
+ result = MultiplyMod(result, multiplier, modulus);
index = windowStart - 1;
}
return result;
}
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static ulong MultiplyMod(ulong left, ulong right, ulong modulus)
+ {
+ // GpuUInt128.MulMod is the fastest GPU-compatible multiply-reduce helper available to this kernel path.
+ var product = new GpuUInt128(left);
+ return product.MulMod(right, modulus);
+ }
+
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong ComputeWindowedOddPower(ulong windowValue, ulong modulus)
{
@@ -172,18 +174,19 @@ internal struct ExponentRemainderStepperGpu
ulong baseValue = 2UL;
if (windowValue == 1UL)
{
+ // Single-bit windows occur frequently in production primes, so keep the fast return for window value 1.
return baseValue;
}
ulong result = baseValue;
ulong remaining = (windowValue - 1UL) >> 1;
- ulong squareBase = baseValue.MulMod64Gpu(baseValue, modulus);
+ ulong squareBase = MultiplyMod(baseValue, baseValue, modulus);
while (remaining != 0UL)
{
- result = ((remaining & 1UL) != 0UL) ? result.MulMod64Gpu(squareBase, modulus) : result;
+ result = ((remaining & 1UL) != 0UL) ? MultiplyMod(result, squareBase, modulus) : result;
remaining >>= 1;
- squareBase = (remaining != 0UL) ? squareBase.MulMod64Gpu(squareBase, modulus) : squareBase;
+ squareBase = (remaining != 0UL) ? MultiplyMod(squareBase, squareBase, modulus) : squareBase;
}
return result;