From 6887111128ebc053c5384d08dfba5eb6baca06f1 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Sun, 7 Dec 2025 19:16:45 +0100 Subject: [PATCH] BigInteger-friendly steppers added. MersenneResidueStepper converted to struct. Manual residue tracking replaced with steppers. Modulo operations replaced with steppers. --- .../BigIntegerCycleRemainderStepper.cs | 38 +++ .../BigIntegerExponentStepper.cs | 60 ++++ ...MersenneNumberDivisorByDivisorCpuTester.cs | 269 +++++------------- PerfectNumbers.Core/CycleRemainderStepper.cs | 49 +--- .../ExponentRemainderStepperCpu.cs | 36 +-- .../Gpu/ExponentRemainderStepperGpu.cs | 2 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 120 +------- .../MersenneDivisorResidueStepper.cs | 67 ++--- PerfectNumbers.Core/MersenneResidueStepper.cs | 18 +- 9 files changed, 238 insertions(+), 421 deletions(-) create mode 100644 PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs create mode 100644 PerfectNumbers.Core/BigIntegerExponentStepper.cs diff --git a/PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs b/PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs new file mode 100644 index 0000000..3238b3c --- /dev/null +++ b/PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs @@ -0,0 +1,38 @@ +using System.Numerics; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +/// +/// Maintains rolling remainders of an increasing sequence modulo a BigInteger cycle length. +/// Mirrors for values that exceed 64 bits. +/// +[method: MethodImpl(MethodImplOptions.AggressiveInlining)] +public struct BigIntegerCycleRemainderStepper(BigInteger cycleLength) +{ + private readonly BigInteger _cycleLength = cycleLength; + private BigInteger _previousValue = BigInteger.Zero; + private BigInteger _currentRemainder = BigInteger.Zero; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public BigInteger Initialize(BigInteger value) + { + _previousValue = value; + _currentRemainder = value % _cycleLength; + return _currentRemainder; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public BigInteger ComputeNext(BigInteger value) + { + BigInteger delta = value - _previousValue; + _previousValue = value; + _currentRemainder += delta; + if (_currentRemainder >= _cycleLength) + { + _currentRemainder %= _cycleLength; + } + + return _currentRemainder; + } +} diff --git a/PerfectNumbers.Core/BigIntegerExponentStepper.cs b/PerfectNumbers.Core/BigIntegerExponentStepper.cs new file mode 100644 index 0000000..de4bd99 --- /dev/null +++ b/PerfectNumbers.Core/BigIntegerExponentStepper.cs @@ -0,0 +1,60 @@ +using System.Numerics; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +/// +/// Tracks successive residues of baseValue^exponent mod modulus for ascending exponents using inputs. +/// This avoids recomputing full modular exponentiation when stepping forward by a known delta. +/// +public struct BigIntegerExponentStepper(BigInteger modulus, BigInteger baseValue) +{ + private readonly BigInteger _modulus = modulus; + private readonly BigInteger _baseValue = baseValue % modulus; + private BigInteger PreviousExponent = BigInteger.Zero; + private BigInteger CurrentResidue = BigInteger.One % modulus; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public BigInteger Initialize(BigInteger exponent) + { + CurrentResidue = BigInteger.ModPow(_baseValue, exponent, _modulus); + PreviousExponent = exponent; + return CurrentResidue; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool InitializeIsUnity(BigInteger exponent) + { + Initialize(exponent); + return CurrentResidue.IsOne; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public BigInteger ComputeNext(BigInteger exponent) + { + if (exponent == PreviousExponent) + { + return CurrentResidue; + } + + BigInteger delta = exponent - PreviousExponent; + BigInteger deltaResidue = BigInteger.ModPow(_baseValue, delta, _modulus); + CurrentResidue = (CurrentResidue * deltaResidue) % _modulus; + PreviousExponent = exponent; + return CurrentResidue; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool ComputeNextIsUnity(BigInteger exponent) + { + ComputeNext(exponent); + return CurrentResidue.IsOne; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Reset() + { + PreviousExponent = BigInteger.Zero; + CurrentResidue = BigInteger.One % _modulus; + } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index f34b688..5d80c20 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -2,6 +2,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using System.Globalization; using System.IO; +using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; @@ -294,40 +295,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ulong limitHigh = limit.High; ulong divisorLow = divisor.Low; - // Keep the divisibility filters aligned with the divisor-cycle generator so the - // CPU path never requests cycles that were skipped during cache creation. - LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); - - byte step10; - byte step8; - byte step3; - byte step7; - byte step11; - byte step13; - byte step17; - byte step19; - - byte remainder10; - byte remainder8; - byte remainder3; - byte remainder7; - byte remainder11; - byte remainder13; - byte remainder17; - byte remainder19; - if (stepHigh == 0UL && limitHigh == 0UL) { - step10 = (byte)(stepLow % 10UL); - step8 = (byte)(stepLow & 7UL); - step3 = (byte)(stepLow % 3UL); - step7 = (byte)(stepLow % 7UL); - step11 = (byte)(stepLow % 11UL); - step13 = (byte)(stepLow % 13UL); - step17 = (byte)(stepLow % 17UL); - step19 = (byte)(stepLow % 19UL); - ulong maxK = allowedMax > 0UL ? (allowedMax - 1UL) / stepLow : 0UL; if (maxK == 0UL) { @@ -350,15 +319,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv allowedMax, startK, maxK, - decimalMask, - step10, - step8, - step3, - step7, - step11, - step13, - step17, - step19, ref currentK, out processedTop, out foundDivisor); @@ -381,15 +341,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv allowedMax, 1UL, lowerEnd, - decimalMask, - step10, - step8, - step3, - step7, - step11, - step13, - step17, - step19, ref currentK, out processedBottom, out foundDivisor); @@ -400,29 +351,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } ulong divisorHigh = divisor.High; - step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL); - step8 = (byte)(stepLow % 8UL); - step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL); - step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL); - step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL); - step13 = (byte)((((stepHigh % 13UL) * 3UL) + (stepLow % 13UL)) % 13UL); - step17 = (byte)(((stepHigh % 17UL) + (stepLow % 17UL)) % 17UL); - step19 = (byte)((((stepHigh % 19UL) * 17UL) + (stepLow % 19UL)) % 19UL); - - remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL); - remainder8 = (byte)(divisorLow % 8UL); - remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL); - remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL); - remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL); - remainder13 = (byte)((((divisorHigh % 13UL) * 3UL) + (divisorLow % 13UL)) % 13UL); - remainder17 = (byte)(((divisorHigh % 17UL) + (divisorLow % 17UL)) % 17UL); - remainder19 = (byte)((((divisorHigh % 19UL) * 17UL) + (divisorLow % 19UL)) % 19UL); + var residueStepper = new MersenneDivisorResidueStepper(prime, step, divisor); var divisorPool = MontgomeryDivisorDataPool.Shared; while (divisor.CompareTo(limit) <= 0) { - bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; - if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + if (residueStepper.IsAdmissible()) { ulong candidate = divisor.Low; MontgomeryDivisorData divisorData = divisorPool.FromModulus(candidate); @@ -468,14 +402,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv divisor.Add(step); currentK++; - remainder10 = AddMod10(remainder10, step10); - remainder8 = AddMod8(remainder8, step8); - remainder3 = AddMod3(remainder3, step3); - remainder7 = AddMod7(remainder7, step7); - remainder11 = AddMod11(remainder11, step11); - remainder13 = AddMod13(remainder13, step13); - remainder17 = AddMod17(remainder17, step17); - remainder19 = AddMod19(remainder19, step19); + residueStepper.Advance(); } processedAll = true; @@ -510,27 +437,29 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); - byte step10 = (byte)(step % 10); - byte step8 = (byte)(step % 8); - byte step3 = (byte)(step % 3); - byte step7 = (byte)(step % 7); - byte step11 = (byte)(step % 11); - byte step13 = (byte)(step % 13); - byte step17 = (byte)(step % 17); - byte step19 = (byte)(step % 19); + var rem10 = new BigIntegerCycleRemainderStepper(10); + var rem8 = new BigIntegerCycleRemainderStepper(8); + var rem5 = new BigIntegerCycleRemainderStepper(5); + var rem3 = new BigIntegerCycleRemainderStepper(3); + var rem7 = new BigIntegerCycleRemainderStepper(7); + var rem11 = new BigIntegerCycleRemainderStepper(11); + var rem13 = new BigIntegerCycleRemainderStepper(13); + var rem17 = new BigIntegerCycleRemainderStepper(17); + var rem19 = new BigIntegerCycleRemainderStepper(19); - byte remainder10 = (byte)(divisor % 10); - byte remainder8 = (byte)(divisor % 8); - byte remainder3 = (byte)(divisor % 3); - byte remainder7 = (byte)(divisor % 7); - byte remainder11 = (byte)(divisor % 11); - byte remainder13 = (byte)(divisor % 13); - byte remainder17 = (byte)(divisor % 17); - byte remainder19 = (byte)(divisor % 19); + byte remainder10 = (byte)rem10.Initialize(divisor); + byte remainder8 = (byte)rem8.Initialize(divisor); + byte remainder5 = (byte)rem5.Initialize(divisor); + byte remainder3 = (byte)rem3.Initialize(divisor); + byte remainder7 = (byte)rem7.Initialize(divisor); + byte remainder11 = (byte)rem11.Initialize(divisor); + byte remainder13 = (byte)rem13.Initialize(divisor); + byte remainder17 = (byte)rem17.Initialize(divisor); + byte remainder19 = (byte)rem19.Initialize(divisor); while (divisor <= allowedMax) { - bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; + bool passesSmallModuli = remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) { if (IsProbablePrimeBigInteger(divisor)) @@ -552,14 +481,15 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv currentK += BigInteger.One; divisor += step; - remainder10 = AddMod10(remainder10, step10); - remainder8 = AddMod8(remainder8, step8); - remainder3 = AddMod3(remainder3, step3); - remainder7 = AddMod7(remainder7, step7); - remainder11 = AddMod11(remainder11, step11); - remainder13 = AddMod13(remainder13, step13); - remainder17 = AddMod17(remainder17, step17); - remainder19 = AddMod19(remainder19, step19); + remainder10 = (byte)rem10.ComputeNext(divisor); + remainder8 = (byte)rem8.ComputeNext(divisor); + remainder5 = (byte)rem5.ComputeNext(divisor); + remainder3 = (byte)rem3.ComputeNext(divisor); + remainder7 = (byte)rem7.ComputeNext(divisor); + remainder11 = (byte)rem11.ComputeNext(divisor); + remainder13 = (byte)rem13.ComputeNext(divisor); + remainder17 = (byte)rem17.ComputeNext(divisor); + remainder19 = (byte)rem19.ComputeNext(divisor); if (divisor > allowedMax) { @@ -579,15 +509,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ulong allowedMax, ulong startK, ulong endK, - ushort decimalMask, - byte step10, - byte step8, - byte step3, - byte step7, - byte step11, - byte step13, - byte step17, - byte step19, ref ulong currentK, out bool processedAll, out ulong foundDivisor) @@ -619,14 +540,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ulong startDivisor = (ulong)startDivisor128; ulong rangeLimit = (ulong)rangeLimit128; - byte remainder10 = (byte)(startDivisor % 10UL); - byte remainder8 = (byte)(startDivisor & 7UL); - byte remainder3 = (byte)(startDivisor % 3UL); - byte remainder7 = (byte)(startDivisor % 7UL); - byte remainder11 = (byte)(startDivisor % 11UL); - byte remainder13 = (byte)(startDivisor % 13UL); - byte remainder17 = (byte)(startDivisor % 17UL); - byte remainder19 = (byte)(startDivisor % 19UL); + var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)step128, (GpuUInt128)startDivisor128); return CheckDivisors64( gpu, @@ -634,23 +548,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv step, rangeLimit, startDivisor, - decimalMask, - step10, - step8, - step3, - step7, - step11, - step13, - step17, - step19, - remainder10, - remainder8, - remainder3, - remainder7, - remainder11, - remainder13, - remainder17, - remainder19, + ref residueStepper, ref currentK, out processedAll, out foundDivisor); @@ -662,23 +560,7 @@ private bool CheckDivisors64( ulong step, ulong limit, ulong divisor, - ushort decimalMask, - byte step10, - byte step8, - byte step3, - byte step7, - byte step11, - byte step13, - byte step17, - byte step19, - byte remainder10, - byte remainder8, - byte remainder3, - byte remainder7, - byte remainder11, - byte remainder13, - byte remainder17, - byte remainder19, + ref MersenneDivisorResidueStepper residueStepper, ref ulong currentK, out bool processedAll, out ulong foundDivisor) @@ -687,44 +569,41 @@ private bool CheckDivisors64( Queue divisorPool = MontgomeryDivisorDataPool.Shared; if (!canAdvance) { - if (divisor <= limit) + if (divisor <= limit && residueStepper.IsAdmissible()) { - if (remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0 && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); + ulong divisorCycle; + if (!_tryCalculateCycleLengthForExponent( + divisor, + gpu, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) { - MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); - ulong divisorCycle; - if (!_tryCalculateCycleLengthForExponent( - divisor, - gpu, - prime, - divisorData, - out ulong computedCycle, - out bool primeOrderFailed) || computedCycle == 0UL) - { - divisorCycle = _calculateCycleLength( - divisor, - gpu, - divisorData, - skipPrimeOrderHeuristic: primeOrderFailed); - } - else - { - divisorCycle = computedCycle; - } + divisorCycle = _calculateCycleLength( + divisor, + gpu, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } - divisorPool.Return(divisorData); + divisorPool.Return(divisorData); - if (divisorCycle == prime) - { - foundDivisor = divisor; - processedAll = true; - return true; - } + if (divisorCycle == prime) + { + foundDivisor = divisor; + processedAll = true; + return true; + } - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); } } @@ -743,8 +622,7 @@ private bool CheckDivisors64( ulong remainingIterations = ((limit - divisor) / step) + 1UL; while (true) { - bool passesSmallModuli = remainder3 != 0 && remainder7 != 0 && remainder11 != 0 && remainder13 != 0 && remainder17 != 0 && remainder19 != 0; - if (passesSmallModuli && (remainder8 == 1 || remainder8 == 7) && ((decimalMask >> remainder10) & 1) != 0) + if (residueStepper.IsAdmissible()) { MontgomeryDivisorData divisorData = divisorPool.FromModulus(divisor); ulong divisorCycle; @@ -792,14 +670,7 @@ private bool CheckDivisors64( divisor += step; currentK++; - remainder10 = AddMod10(remainder10, step10); - remainder8 = AddMod8(remainder8, step8); - remainder3 = AddMod3(remainder3, step3); - remainder7 = AddMod7(remainder7, step7); - remainder11 = AddMod11(remainder11, step11); - remainder13 = AddMod13(remainder13, step13); - remainder17 = AddMod17(remainder17, step17); - remainder19 = AddMod19(remainder19, step19); + residueStepper.Advance(); } } @@ -870,7 +741,8 @@ private bool CheckDivisors64( continue; } - BigInteger x = BigInteger.ModPow(baseValue, d, value); + var expStepper = new BigIntegerExponentStepper(value, baseValue); + BigInteger x = expStepper.Initialize(d); if (x == 1 || x == value - 1) { continue; @@ -879,7 +751,8 @@ private bool CheckDivisors64( bool witnessFound = true; for (int r = 1; r < s; r++) { - x = BigInteger.ModPow(x, 2, value); + BigInteger targetExponent = d << r; + x = expStepper.ComputeNext(targetExponent); if (x == value - 1) { witnessFound = false; diff --git a/PerfectNumbers.Core/CycleRemainderStepper.cs b/PerfectNumbers.Core/CycleRemainderStepper.cs index a9e809a..a495bd3 100644 --- a/PerfectNumbers.Core/CycleRemainderStepper.cs +++ b/PerfectNumbers.Core/CycleRemainderStepper.cs @@ -2,37 +2,20 @@ using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; -public struct CycleRemainderStepper +public struct CycleRemainderStepper(ulong cycleLength) { - private readonly ulong _cycleLength; - private ulong _previousPrime; - private ulong _currentRemainder; - private bool _hasState; - - public CycleRemainderStepper(ulong cycleLength) - { - if (cycleLength == 0UL) - { - throw new ArgumentOutOfRangeException(nameof(cycleLength), "Cycle length must be non-zero."); - } - - _cycleLength = cycleLength; - _previousPrime = 0UL; - _currentRemainder = 0UL; - _hasState = false; - } + private readonly ulong _cycleLength = cycleLength; + private ulong _previousPrime = 0UL; + private ulong _currentRemainder = 0UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Reset() { - // TODO: Inline this reset at the call sites so the hot loops reuse struct reinitialization - // measured fastest in MersenneDivisorCycleLengthGpuBenchmarks, avoiding the extra - // method call when scanners need to restart stepping. _previousPrime = 0UL; _currentRemainder = 0UL; - _hasState = false; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong Initialize(ulong prime) { ulong cycleLength = _cycleLength; @@ -44,24 +27,15 @@ public struct CycleRemainderStepper // MersenneDivisorCycleLengthGpuBenchmarks winner for both CPU and GPU call sites. remainder %= cycleLength; } + _previousPrime = prime; _currentRemainder = remainder; - _hasState = true; return remainder; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong ComputeNext(ulong prime) { - if (!_hasState) - { - throw new InvalidOperationException("CycleRemainderStepper must be initialized before computing."); - } - - if (prime <= _previousPrime) - { - throw new ArgumentOutOfRangeException(nameof(prime), "Primes must be processed in strictly increasing order."); - } - ulong cycleLength = _cycleLength; ulong delta = prime - _previousPrime; _previousPrime = prime; @@ -90,13 +64,4 @@ public struct CycleRemainderStepper } return _currentRemainder; } - - public void CaptureState(out ulong previousPrime, out ulong previousRemainder, out bool hasState) - { - // TODO: Expose these fields directly once the residue scanners adopt the single-cycle helper - // so the hot path can read them without paying for an additional wrapper call per iteration. - previousPrime = _previousPrime; - previousRemainder = _currentRemainder; - hasState = _hasState; - } } diff --git a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs index 83e960f..d486962 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs @@ -1,28 +1,21 @@ +using System.Runtime.CompilerServices; + namespace PerfectNumbers.Core; /// /// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path. /// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order. /// -internal struct ExponentRemainderStepperCpu +internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor) { - private readonly MontgomeryDivisorData _divisor; - private readonly ulong _modulus; - private readonly ulong _nPrime; - private readonly ulong _montgomeryOne; - public ulong PreviousExponent; - private ulong _currentMontgomery; - - public ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor) - { - _divisor = divisor; - _modulus = divisor.Modulus; - _nPrime = divisor.NPrime; - _montgomeryOne = divisor.MontgomeryOne; - PreviousExponent = 0UL; - _currentMontgomery = divisor.MontgomeryOne; - } +private readonly MontgomeryDivisorData _divisor = divisor; + private readonly ulong _modulus = divisor.Modulus; + private readonly ulong _nPrime = divisor.NPrime; + private readonly ulong _montgomeryOne = divisor.MontgomeryOne; + private ulong PreviousExponent = 0UL; + private ulong _currentMontgomery = divisor.MontgomeryOne; + [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly bool MatchesDivisor(in MontgomeryDivisorData divisor) => _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne; @@ -30,6 +23,7 @@ internal struct ExponentRemainderStepperCpu /// 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. /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] public void Reset() { PreviousExponent = 0UL; @@ -42,6 +36,7 @@ internal struct ExponentRemainderStepperCpu /// /// First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata. /// The canonical residue produced by 2^exponent mod divisor. + [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong InitializeCpu(ulong exponent) { InitializeCpuState(exponent); @@ -53,6 +48,7 @@ internal struct ExponentRemainderStepperCpu /// /// 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. + [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool InitializeCpuIsUnity(ulong exponent) { InitializeCpuState(exponent); @@ -65,6 +61,7 @@ internal struct ExponentRemainderStepperCpu /// /// 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. + [MethodImpl(MethodImplOptions.AggressiveInlining)] public ulong ComputeNextCpu(ulong exponent) { ulong delta = exponent - PreviousExponent; @@ -84,6 +81,7 @@ internal struct ExponentRemainderStepperCpu /// /// Next exponent in the ascending sequence. /// when the canonical residue equals one. + [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool ComputeNextIsUnity(ulong exponent) { ulong delta = exponent - PreviousExponent; @@ -103,6 +101,7 @@ internal struct ExponentRemainderStepperCpu /// 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. + [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity) { _currentMontgomery = montgomeryResult; @@ -119,6 +118,7 @@ internal struct ExponentRemainderStepperCpu /// Montgomery residue delta corresponding to the exponent gap. /// Outputs whether the resulting residue equals Montgomery one. /// when the advance succeeds. + [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity) { // TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead @@ -130,11 +130,13 @@ internal struct ExponentRemainderStepperCpu return true; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private void InitializeCpuState(ulong exponent) { _currentMontgomery = exponent.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); PreviousExponent = exponent; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private readonly ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime); } diff --git a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs index 2cc7926..0c3f5a3 100644 --- a/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs +++ b/PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs @@ -11,7 +11,7 @@ internal struct ExponentRemainderStepperGpu { private readonly ulong _modulus; private readonly ReadOnlyGpuUInt128 _modulusWide; - public ulong PreviousExponent; + private ulong PreviousExponent; private GpuUInt128 _currentResidue; [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 458993c..5362d5a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -5,6 +5,7 @@ using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; +using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core.Gpu; @@ -253,16 +254,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; ulong startK = minK < 1UL ? 1UL : minK; - ulong stepValue = (ulong)twoP128; - byte step10 = (byte)(stepValue % 10UL); - byte step8 = (byte)(stepValue & 7UL); - byte step3 = (byte)(stepValue % 3UL); - byte step7 = (byte)(stepValue % 7UL); - byte step11 = (byte)(stepValue % 11UL); - - LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - ushort decimalMask = DivisorGenerator.GetDecimalMask(lastDigit); - Span divisorStorage = divisors.AsSpan(); Span exponentStorage = exponents.AsSpan(); Span divisorDataStorage = divisorData.AsSpan(); @@ -289,12 +280,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN ulong prime, UInt128 twoP128, UInt128 allowedMax128, - byte step10, - byte step8, - byte step3, - byte step7, - byte step11, - ushort decimalMask, ArrayView1D divisorDataView, ArrayView1D offsetView, ArrayView1D countView, @@ -329,6 +314,7 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN return true; } + var residueStepper = new MersenneDivisorResidueStepper(prime, (GpuUInt128)twoP128, (GpuUInt128)startDivisor128); UInt128 rangeLimit128 = (twoP128 * rangeEndK) + UInt128.One; if (rangeLimit128 > allowedMax128) { @@ -338,12 +324,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN UInt128 rangeCount128 = ((rangeLimit128 - startDivisor128) / twoP128) + UInt128.One; ulong rangeRemaining = rangeCount128 > ulong.MaxValue ? ulong.MaxValue : (ulong)rangeCount128; UInt128 currentRangeDivisor128 = startDivisor128; - ulong currentRangeDivisorValue = (ulong)currentRangeDivisor128; - byte rangeRemainder10 = (byte)(currentRangeDivisorValue % 10UL); - byte rangeRemainder8 = (byte)(currentRangeDivisorValue & 7UL); - byte rangeRemainder3 = (byte)(currentRangeDivisorValue % 3UL); - byte rangeRemainder7 = (byte)(currentRangeDivisorValue % 7UL); - byte rangeRemainder11 = (byte)(currentRangeDivisorValue % 11UL); while (rangeRemaining > 0UL && !composite) { @@ -356,16 +336,11 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN UInt128 nextDivisor128 = currentRangeDivisor128; int admissibleCount = 0; - byte localRemainder10 = rangeRemainder10; - byte localRemainder8 = rangeRemainder8; - byte localRemainder3 = rangeRemainder3; - byte localRemainder7 = rangeRemainder7; - byte localRemainder11 = rangeRemainder11; + var localStepper = residueStepper; for (int i = 0; i < chunkCount; i++) { - bool passesSmallModuli = localRemainder3 != 0 && localRemainder7 != 0 && localRemainder11 != 0; - if (passesSmallModuli && (localRemainder8 == 1 || localRemainder8 == 7) && ((decimalMask >> localRemainder10) & 1) != 0) + if (localStepper.IsAdmissible()) { ulong divisorValue = (ulong)nextDivisor128; MontgomeryDivisorData montgomeryData = divisorPool.FromModulus(divisorValue); @@ -389,18 +364,10 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN } nextDivisor128 += twoP128; - localRemainder10 = AddMod10(localRemainder10, step10); - localRemainder8 = AddMod8(localRemainder8, step8); - localRemainder3 = AddMod3(localRemainder3, step3); - localRemainder7 = AddMod7(localRemainder7, step7); - localRemainder11 = AddMod11(localRemainder11, step11); + localStepper.Advance(); } - rangeRemainder10 = localRemainder10; - rangeRemainder8 = localRemainder8; - rangeRemainder3 = localRemainder3; - rangeRemainder7 = localRemainder7; - rangeRemainder11 = localRemainder11; + residueStepper = localStepper; processedCountLocal += (ulong)chunkCount; @@ -456,12 +423,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN prime, twoP128, allowedMax128, - step10, - step8, - step3, - step7, - step11, - decimalMask, divisorDataView, offsetView, countView, @@ -502,12 +463,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN prime, twoP128, allowedMax128, - step10, - step8, - step3, - step7, - step11, - decimalMask, divisorDataView, offsetView, countView, @@ -550,69 +505,6 @@ public sealed partial class MersenneNumberDivisorByDivisorGpuTester : IMersenneN return computedCycle; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod3(byte value, byte delta) - { - const int Modulus = 3; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod7(byte value, byte delta) - { - const int Modulus = 7; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod10(byte value, byte delta) - { - const int Modulus = 10; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod11(byte value, byte delta) - { - const int Modulus = 11; - int sum = value + delta; - - if (sum >= Modulus) - { - sum -= Modulus; - } - - return (byte)sum; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte AddMod8(byte value, byte delta) - { - return (byte)((value + delta) & 7); - } - private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime) { diff --git a/PerfectNumbers.Core/MersenneDivisorResidueStepper.cs b/PerfectNumbers.Core/MersenneDivisorResidueStepper.cs index e4a1bf6..a7492a0 100644 --- a/PerfectNumbers.Core/MersenneDivisorResidueStepper.cs +++ b/PerfectNumbers.Core/MersenneDivisorResidueStepper.cs @@ -3,52 +3,34 @@ using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; -internal struct MersenneDivisorResidueStepper +[method: MethodImpl(MethodImplOptions.AggressiveInlining)] +internal struct MersenneDivisorResidueStepper(ulong prime, in GpuUInt128 step, in GpuUInt128 firstDivisor) { private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); - private readonly ushort decimalMask; - private readonly byte step10; - private readonly byte step8; - private readonly byte step5; - private readonly byte step3; - private readonly byte step7; - private readonly byte step11; + private readonly ushort decimalMask = (prime & 3UL) == 3UL ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + private readonly byte step10 = ComputeModulo(step, 10, multiplier: 6); + private readonly byte step8 = (byte)(step.Low % 8UL); + private readonly byte step5 = ComputeModulo(step, 5, multiplier: 1); + private readonly byte step3 = ComputeModulo(step, 3, multiplier: 1); + private readonly byte step7 = ComputeModulo(step, 7, multiplier: 2); + private readonly byte step11 = ComputeModulo(step, 11, multiplier: 5); + private readonly byte step13 = ComputeModulo(step, 13, multiplier: 3); + private readonly byte step17 = ComputeModulo(step, 17, multiplier: 1); + private readonly byte step19 = ComputeModulo(step, 19, multiplier: 17); - public byte Remainder10 { get; private set; } + private byte Remainder10 = ComputeModulo(firstDivisor, 10, multiplier: 6); + private byte Remainder8 = (byte)(firstDivisor.Low % 8UL); + private byte Remainder5 = ComputeModulo(firstDivisor, 5, multiplier: 1); + private byte Remainder3 = ComputeModulo(firstDivisor, 3, multiplier: 1); + private byte Remainder7 = ComputeModulo(firstDivisor, 7, multiplier: 2); + private byte Remainder11 = ComputeModulo(firstDivisor, 11, multiplier: 5); + private byte Remainder13 = ComputeModulo(firstDivisor, 13, multiplier: 3); + private byte Remainder17 = ComputeModulo(firstDivisor, 17, multiplier: 1); + private byte Remainder19 = ComputeModulo(firstDivisor, 19, multiplier: 17); - public byte Remainder8 { get; private set; } - - public byte Remainder5 { get; private set; } - - public byte Remainder3 { get; private set; } - - public byte Remainder7 { get; private set; } - - public byte Remainder11 { get; private set; } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public MersenneDivisorResidueStepper(ulong prime, in GpuUInt128 step, in GpuUInt128 firstDivisor) - { - decimalMask = (prime & 3UL) == 3UL ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; - - step10 = ComputeModulo(step, 10, multiplier: 6); - step8 = (byte)(step.Low % 8UL); - step5 = ComputeModulo(step, 5, multiplier: 1); - step3 = ComputeModulo(step, 3, multiplier: 1); - step7 = ComputeModulo(step, 7, multiplier: 2); - step11 = ComputeModulo(step, 11, multiplier: 5); - - Remainder10 = ComputeModulo(firstDivisor, 10, multiplier: 6); - Remainder8 = (byte)(firstDivisor.Low % 8UL); - Remainder5 = ComputeModulo(firstDivisor, 5, multiplier: 1); - Remainder3 = ComputeModulo(firstDivisor, 3, multiplier: 1); - Remainder7 = ComputeModulo(firstDivisor, 7, multiplier: 2); - Remainder11 = ComputeModulo(firstDivisor, 11, multiplier: 5); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] + [MethodImpl(MethodImplOptions.AggressiveInlining)] public readonly bool IsAdmissible() { if (((decimalMask >> Remainder10) & 1) == 0) @@ -61,7 +43,7 @@ internal struct MersenneDivisorResidueStepper return false; } - if (Remainder3 == 0 || Remainder5 == 0 || Remainder7 == 0 || Remainder11 == 0) + if (Remainder3 == 0 || Remainder5 == 0 || Remainder7 == 0 || Remainder11 == 0 || Remainder13 == 0 || Remainder17 == 0 || Remainder19 == 0) { return false; } @@ -78,6 +60,9 @@ internal struct MersenneDivisorResidueStepper Remainder3 = AddMod(Remainder3, step3, 3); Remainder7 = AddMod(Remainder7, step7, 7); Remainder11 = AddMod(Remainder11, step11, 11); + Remainder13 = AddMod(Remainder13, step13, 13); + Remainder17 = AddMod(Remainder17, step17, 17); + Remainder19 = AddMod(Remainder19, step19, 19); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/MersenneResidueStepper.cs b/PerfectNumbers.Core/MersenneResidueStepper.cs index 22aa270..5ebdd63 100644 --- a/PerfectNumbers.Core/MersenneResidueStepper.cs +++ b/PerfectNumbers.Core/MersenneResidueStepper.cs @@ -1,19 +1,17 @@ +using System.Runtime.CompilerServices; + namespace PerfectNumbers.Core; -public class MersenneResidueStepper +public struct MersenneResidueStepper { - private UInt128 m; + private readonly UInt128 m; private ulong p; private UInt128 pow2_mod; private UInt128 M_mod; + [MethodImpl(MethodImplOptions.AggressiveInlining)] public MersenneResidueStepper(UInt128 modulus, ulong p0, UInt128 M_mod_m_at_p0) { - if (modulus == 0) - { - throw new ArgumentException("Modulus must be positive."); - } - m = modulus; p = p0; @@ -32,8 +30,10 @@ public class MersenneResidueStepper // TODO: Inline this accessor into callers so the hot residue path reads the backing fields directly // instead of jumping through a wrapper method every iteration. - public (ulong, UInt128) State() => (p, M_mod); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public readonly (ulong, UInt128) State() => (p, M_mod); + [MethodImpl(MethodImplOptions.AggressiveInlining)] public (ulong, UInt128) Step(long delta) { if (delta == 0) @@ -67,6 +67,7 @@ public class MersenneResidueStepper return State(); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public (ulong, UInt128) To(ulong p_target) { // TODO: Collapse this wrapper once callers can invoke Step directly; the extra indirection shows up @@ -76,6 +77,7 @@ public class MersenneResidueStepper return Step(delta); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static UInt128 ModInverse(UInt128 a, UInt128 m) { UInt128 m0 = m;