From 557bc5d06f3cc08ca622471aca3a1fb8a0af34b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Fri, 24 Oct 2025 17:26:49 +0200 Subject: [PATCH] Use field-based exponent tracking for steppers --- .../ExponentRemainderStepperCpu.cs | 22 ++++---- .../Gpu/ExponentRemainderStepperGpu.cs | 51 ++++++++++--------- 2 files changed, 37 insertions(+), 36 deletions(-) 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;