diff --git a/PerfectNumbers.Core/BidirectionalCycleRemainderStepper.cs b/PerfectNumbers.Core/BidirectionalCycleRemainderStepper.cs new file mode 100644 index 0000000..1320e2f --- /dev/null +++ b/PerfectNumbers.Core/BidirectionalCycleRemainderStepper.cs @@ -0,0 +1,62 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +public struct BidirectionalCycleRemainderStepper(ulong modulus) +{ + private readonly ulong _modulus = modulus; + private ulong _previous = 0UL; + private ulong _remainder = 0UL; + + public readonly ulong Modulus + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _modulus; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Initialize(ulong value) + { + _remainder = value.ReduceCycleRemainder(_modulus); + _previous = value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong Step(ulong value) + { + long delta = unchecked((long)value - (long)_previous); + _previous = value; + if (delta == 0) + { + return _remainder; + } + + if (delta > 0) + { + ulong step = ((ulong)delta).ReduceCycleRemainder(_modulus); + ulong sum = _remainder + step; + if (sum >= _modulus) + { + sum -= _modulus; + } + + _remainder = sum; + return _remainder; + } + + ulong backward = ((ulong)-delta).ReduceCycleRemainder(_modulus); + if (backward > _remainder) + { + ulong diff = backward - _remainder; + ulong wrap = diff.ReduceCycleRemainder(_modulus); + _remainder = wrap == 0 ? 0UL : _modulus - wrap; + return _remainder; + } + + _remainder -= backward; + return _remainder; + } +} diff --git a/PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs b/PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs index 3238b3c..3fdb3d9 100644 --- a/PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs +++ b/PerfectNumbers.Core/BigIntegerCycleRemainderStepper.cs @@ -18,7 +18,7 @@ public struct BigIntegerCycleRemainderStepper(BigInteger cycleLength) public BigInteger Initialize(BigInteger value) { _previousValue = value; - _currentRemainder = value % _cycleLength; + _currentRemainder = value.ReduceCycleRemainder(_cycleLength); return _currentRemainder; } @@ -27,11 +27,7 @@ public struct BigIntegerCycleRemainderStepper(BigInteger cycleLength) { BigInteger delta = value - _previousValue; _previousValue = value; - _currentRemainder += delta; - if (_currentRemainder >= _cycleLength) - { - _currentRemainder %= _cycleLength; - } + _currentRemainder = (_currentRemainder + delta).ReduceCycleRemainder(_cycleLength); return _currentRemainder; } diff --git a/PerfectNumbers.Core/BigIntegerExtensions.cs b/PerfectNumbers.Core/BigIntegerExtensions.cs new file mode 100644 index 0000000..a97647c --- /dev/null +++ b/PerfectNumbers.Core/BigIntegerExtensions.cs @@ -0,0 +1,113 @@ +using System.Numerics; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +/// +/// BigInteger helpers including cycle-remainder reduction using the subtract ladder that benchmarked best. +/// Shared stats/warm set stay here so ulong reductions can reuse them through ULongExtensions. +/// +public static partial class BigIntegerExtensions +{ + // private const int UnrolledSubtractCount = 6; + // private static readonly HashSet WarmModuli; + // private static long _hits; + // private static long _misses; + + // static BigIntegerExtensions() + // { + // ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + // WarmModuli = new HashSet(snapshot.Length); + // for (int i = 0; i < snapshot.Length; i++) + // { + // ulong cycle = snapshot[i]; + // if (cycle > 1UL) + // { + // WarmModuli.Add(cycle); + // } + // } + // } + + // public static void DumpCycleRemainderStats() + // { + // long hits = Interlocked.Read(ref _hits); + // long misses = Interlocked.Read(ref _misses); + // Console.WriteLine($"CycleRemainderReducerCache: hits={hits}, misses={misses}"); + // } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static BigInteger ReduceCycleRemainder(this BigInteger value, BigInteger modulus) + { + // Track(modulus); + // DumpCycleRemainderStats(); + + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + return value % modulus; + } + + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // internal static void Track(in BigInteger modulus) + // { + // if (modulus <= ulong.MaxValue) + // { + // Track((ulong)modulus); + // } + // else + // { + // Interlocked.Increment(ref _misses); + // } + // } + + // [MethodImpl(MethodImplOptions.AggressiveInlining)] + // internal static void Track(ulong modulus) + // { + // if (WarmModuli.Contains(modulus)) + // { + // Interlocked.Increment(ref _hits); + // } + // else + // { + // Interlocked.Increment(ref _misses); + // } + // } +} diff --git a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs index 85db934..a04b62b 100644 --- a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs +++ b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs @@ -55,7 +55,7 @@ internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelera // 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 ExponentRemainderStepperCpu(divisorData); + var exponentStepper = new ExponentRemainderStepperCpu(divisorData, divisorCycle); var cycleStepper = new CycleRemainderStepper(divisorCycle); diff --git a/PerfectNumbers.Core/CycleRemainderStepper.cs b/PerfectNumbers.Core/CycleRemainderStepper.cs index a495bd3..b9297c0 100644 --- a/PerfectNumbers.Core/CycleRemainderStepper.cs +++ b/PerfectNumbers.Core/CycleRemainderStepper.cs @@ -19,14 +19,7 @@ public struct CycleRemainderStepper(ulong cycleLength) public ulong Initialize(ulong prime) { ulong cycleLength = _cycleLength; - ulong remainder = prime; - if (remainder >= cycleLength) - { - // TODO: Swap this `%` for the shared divisor-cycle remainder helper so initialization reuses the cached - // cycle deltas benchmarked faster than on-the-fly modulo work, matching the - // MersenneDivisorCycleLengthGpuBenchmarks winner for both CPU and GPU call sites. - remainder %= cycleLength; - } + ulong remainder = prime.ReduceCycleRemainder(cycleLength); _previousPrime = prime; _currentRemainder = remainder; @@ -43,9 +36,7 @@ public struct CycleRemainderStepper(ulong cycleLength) if (ulong.MaxValue - _currentRemainder < delta) { UInt128 extended = (UInt128)_currentRemainder + delta; - // TODO: Replace this `%` with the UInt128-aware cycle reducer so large deltas use the cached subtraction - // ladder instead of falling back to the slower modulo implementation. - _currentRemainder = (ulong)(extended % cycleLength); + _currentRemainder = ((ulong)(extended % cycleLength)).ReduceCycleRemainder(cycleLength); return _currentRemainder; } @@ -56,10 +47,7 @@ public struct CycleRemainderStepper(ulong cycleLength) _currentRemainder -= cycleLength; if (_currentRemainder >= cycleLength) { - // TODO: Route this `%` through the shared divisor-cycle helper so repeated wrap-arounds avoid - // modulo operations and match the benchmarked fast path highlighted in - // MersenneDivisorCycleLengthGpuBenchmarks. - _currentRemainder %= cycleLength; + _currentRemainder = _currentRemainder.ReduceCycleRemainder(cycleLength); } } return _currentRemainder; diff --git a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs index d486962..7e03a3b 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs @@ -1,3 +1,4 @@ +using System.Collections.Concurrent; using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; @@ -6,18 +7,20 @@ 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(in MontgomeryDivisorData divisor) +internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL) { -private readonly MontgomeryDivisorData _divisor = divisor; + private readonly MontgomeryDivisorData _divisor = divisor; private readonly ulong _modulus = divisor.Modulus; private readonly ulong _nPrime = divisor.NPrime; private readonly ulong _montgomeryOne = divisor.MontgomeryOne; + private readonly ulong _cycleLength = cycleLength; private ulong PreviousExponent = 0UL; private ulong _currentMontgomery = divisor.MontgomeryOne; + private static readonly ConcurrentDictionary<(ulong Cycle, ulong Delta), ulong> _cycleDeltaCache = new(); [MethodImpl(MethodImplOptions.AggressiveInlining)] - public readonly bool MatchesDivisor(in MontgomeryDivisorData divisor) - => _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne; + public readonly bool MatchesDivisor(in MontgomeryDivisorData divisor, ulong cycleLength) + => _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne && _cycleLength == cycleLength; /// /// Clears the cached Montgomery residue so the stepper can be reused for another divisor. @@ -61,14 +64,14 @@ private readonly MontgomeryDivisorData _divisor = divisor; /// /// 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)] + [MethodImpl(MethodImplOptions.AggressiveInlining)] 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 - // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks. - ulong multiplier = delta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); + ulong reducedDelta = ReduceDelta(delta); + ulong multiplier = reducedDelta == 0UL + ? _montgomeryOne + : reducedDelta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime); PreviousExponent = exponent; return ReduceCurrent(); @@ -85,9 +88,10 @@ private readonly MontgomeryDivisorData _divisor = divisor; public bool ComputeNextIsUnity(ulong exponent) { 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.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); + ulong reducedDelta = ReduceDelta(delta); + ulong multiplier = reducedDelta == 0UL + ? _montgomeryOne + : reducedDelta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime); PreviousExponent = exponent; return _currentMontgomery == _montgomeryOne; @@ -119,24 +123,62 @@ private readonly MontgomeryDivisorData _divisor = divisor; /// 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 - // 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.MontgomeryMultiplyCpu(montgomeryDelta, _modulus, _nPrime); + public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity) + { + ulong delta = exponent - PreviousExponent; + ulong reducedDelta = ReduceDelta(delta); + + // Prefer a cycle-reduced powmod so stepping stays aligned with the divisor-cycle cache. + ulong multiplier; + if (reducedDelta == 0UL) + { + multiplier = _montgomeryOne; + } + else if (_cycleLength != 0UL) + { + multiplier = reducedDelta.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); + } + else + { + multiplier = montgomeryDelta; + } + + _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime); PreviousExponent = exponent; isUnity = _currentMontgomery == _montgomeryOne; return true; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] + [MethodImpl(MethodImplOptions.AggressiveInlining)] private void InitializeCpuState(ulong exponent) { - _currentMontgomery = exponent.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); + ulong rotation = ReduceDelta(exponent); + _currentMontgomery = rotation == 0UL + ? _montgomeryOne + : rotation.Pow2MontgomeryModWindowedKeepMontgomeryCpu(_divisor); PreviousExponent = exponent; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private readonly ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private readonly ulong ReduceDelta(ulong delta) + { + ulong cycleLength = _cycleLength; + if (cycleLength == 0UL) + { + return delta; + } + + var key = (cycleLength, delta); + if (_cycleDeltaCache.TryGetValue(key, out ulong cached)) + { + return cached; + } + + ulong reduced = delta.ReduceCycleRemainder(cycleLength); + _cycleDeltaCache.TryAdd(key, reduced); + return reduced; + } } diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index b5f46e5..82977b1 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -30,6 +30,23 @@ public sealed class HeuristicPrimeTester private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + [ThreadStatic] + private static BidirectionalCycleRemainderStepper _mod10Stepper; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte GetMod10Remainder(ulong value) + { + if (_mod10Stepper.Modulus != 10UL) + { + _mod10Stepper = new BidirectionalCycleRemainderStepper(10UL); + _mod10Stepper.Initialize(value); + return (byte)_mod10Stepper.Step(value); + } + + // Caller must initialize once per thread before first use. + return (byte)_mod10Stepper.Step(value); + } + internal enum HeuristicDivisorGroup : byte { None = 0, @@ -70,7 +87,7 @@ public sealed class HeuristicPrimeTester [MethodImpl(MethodImplOptions.AggressiveInlining)] public static bool IsPrimeCpu(ulong n) { - byte nMod10 = (byte)n.Mod10(); + byte nMod10 = GetMod10Remainder(n); ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); if (maxDivisorSquare < 9UL) @@ -85,7 +102,7 @@ public sealed class HeuristicPrimeTester [MethodImpl(MethodImplOptions.AggressiveInlining)] public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n) { - byte nMod10 = (byte)n.Mod10(); + byte nMod10 = GetMod10Remainder(n); ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); return IsPrimeGpu(gpu, n, maxDivisorSquare, nMod10); } @@ -132,11 +149,11 @@ public sealed class HeuristicPrimeTester break; } - ulong divisor = groupADivisors[groupAIndex]; - if (n % divisor == 0UL) - { - return false; - } + ulong divisor = groupADivisors[groupAIndex]; + if (n.ReduceCycleRemainder(divisor) == 0UL) + { + return false; + } groupAIndex++; processed++; @@ -186,10 +203,10 @@ public sealed class HeuristicPrimeTester } ulong divisor = groupADivisors[groupAIndex]; - if (n % divisor == 0UL) - { - return false; - } + if (n.ReduceCycleRemainder(divisor) == 0UL) + { + return false; + } groupAIndex++; processed++; @@ -220,7 +237,7 @@ public sealed class HeuristicPrimeTester break; } - if (n % divisor == 0UL) + if (n.ReduceCycleRemainder(divisor) == 0UL) { return false; } @@ -241,7 +258,6 @@ public sealed class HeuristicPrimeTester private static bool HeuristicTrialDivisionGpuDetectsDivisor(PrimeOrderCalculatorAccelerator gpu, ulong n, ulong maxDivisorSquare, byte nMod10) { - // GpuPrimeWorkLimiter.Acquire(); int acceleratorIndex = gpu.AcceleratorIndex; var kernelLauncher = gpu.HeuristicCombinedTrialDivisionKernelLauncher; bool compositeDetected; @@ -252,13 +268,6 @@ public sealed class HeuristicPrimeTester var stream = AcceleratorStreamPool.Rent(acceleratorIndex); flagView1DView.CopyFromCPU(stream, ref compositeFlag, 1); - // Index1D index, - // ArrayView resultFlag, - // ulong n, - // ulong maxDivisorSquare, - // HeuristicGpuDivisorTableKind tableKind, - // HeuristicGpuDivisorTables tables) - kernelLauncher( stream, 1, @@ -271,48 +280,9 @@ public sealed class HeuristicPrimeTester flagView1DView.CopyToCPU(stream, ref compositeFlag, 1); stream.Synchronize(); - compositeDetected = compositeFlag != 0; - // compositeDetected = compositeFlag != 0; - - // if (!compositeDetected) - // { - // ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); - // for (int i = 0; i < endingOrder.Length && !compositeDetected; i++) - // { - // byte ending = endingOrder[i]; - // HeuristicGpuDivisorTableKind tableKind = ending switch - // { - // 1 => HeuristicGpuDivisorTableKind.GroupBEnding1, - // 7 => HeuristicGpuDivisorTableKind.GroupBEnding7, - // 9 => HeuristicGpuDivisorTableKind.GroupBEnding9, - // _ => HeuristicGpuDivisorTableKind.GroupA, - // }; - - // if (tableKind == HeuristicGpuDivisorTableKind.GroupA) - // { - // continue; - // } - - // int divisorLength = GetGroupBDivisors(ending).Length; - - // compositeFlag = 0; - // flagView1D.CopyFromCPU(stream, ref compositeFlag, 1); - // kernelLauncher( - // stream, - // 1, - // flagView1D, - // n, - // maxDivisorSquare, - // tableKind, - // gpu.HeuristicGpuTables); - // flagView1D.CopyToCPU(stream, ref compositeFlag, 1); - // stream.Synchronize(); - // compositeDetected = compositeFlag != 0; - // } - // } - AcceleratorStreamPool.Return(acceleratorIndex, stream); - // GpuPrimeWorkLimiter.Release(); + compositeDetected = compositeFlag != 0; + return compositeDetected; } @@ -556,6 +526,10 @@ public sealed class HeuristicPrimeTester private readonly ReadOnlySpan groupADivisorSquares; private int groupAIndex; private Span groupBStates; + private CycleRemainderStepper endingStepper; + private CycleRemainderStepper wheelStepper; + private bool hasRemainders; + private ulong lastDivisor; public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) { @@ -566,6 +540,10 @@ public sealed class HeuristicPrimeTester int count = InitializeGroupBStates(nMod10, groupBBuffer); groupBStates = count == 0 ? Span.Empty : groupBBuffer[..count]; + endingStepper = new CycleRemainderStepper(10UL); + wheelStepper = new CycleRemainderStepper(Wheel210); + hasRemainders = false; + lastDivisor = 0UL; } public bool TryGetNext(out HeuristicDivisorCandidate candidate) @@ -586,13 +564,14 @@ public sealed class HeuristicPrimeTester HeuristicDivisorGroup group = currentIndex < GroupAConstantCount ? HeuristicDivisorGroup.GroupAConstant : HeuristicDivisorGroup.GroupAWheel; + (byte ending, ushort wheelResidue) = UpdateRemainders(divisor); candidate = new HeuristicDivisorCandidate( divisor, group, - (byte)(divisor % 10UL), + ending, 0, - (ushort)(divisor % Wheel210)); + wheelResidue); return true; } @@ -623,14 +602,35 @@ public sealed class HeuristicPrimeTester ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex]; bestState.Advance(); + (byte candidateEnding, ushort wheelResidueBest) = UpdateRemainders(bestCandidate); candidate = new HeuristicDivisorCandidate( bestCandidate, HeuristicDivisorGroup.GroupB, - bestState.Ending, + candidateEnding, bestState.PriorityIndex, - (ushort)(bestCandidate % Wheel210)); + wheelResidueBest); return true; } + + private (byte ending, ushort wheelResidue) UpdateRemainders(ulong divisor) + { + byte ending; + ushort wheelResidue; + if (!hasRemainders) + { + ending = (byte)endingStepper.Initialize(divisor); + wheelResidue = (ushort)wheelStepper.Initialize(divisor); + hasRemainders = true; + } + else + { + ending = (byte)endingStepper.ComputeNext(divisor); + wheelResidue = (ushort)wheelStepper.ComputeNext(divisor); + } + + lastDivisor = divisor; + return (ending, wheelResidue); + } } @@ -684,6 +684,10 @@ public sealed class HeuristicPrimeTester private readonly GpuUInt128 step; private readonly GpuUInt128 limit; private GpuUInt128 current; + private CycleRemainderStepper endingStepper; + private CycleRemainderStepper wheelStepper; + private byte currentEnding; + private ushort currentWheelResidue; private MersenneDivisorResidueStepper residueStepper; private bool active; private ulong processedCount; @@ -705,6 +709,21 @@ public sealed class HeuristicPrimeTester current = hasCandidates ? firstDivisor : GpuUInt128.Zero; residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default; + if (hasCandidates) + { + endingStepper = new CycleRemainderStepper(10UL); + currentEnding = (byte)endingStepper.Initialize(firstDivisor.Low); + + wheelStepper = new CycleRemainderStepper(Wheel210); + currentWheelResidue = (ushort)wheelStepper.Initialize(firstDivisor.Low); + } + else + { + endingStepper = default; + wheelStepper = default; + currentEnding = 0; + currentWheelResidue = 0; + } active = hasCandidates; processedCount = 0UL; lastDivisor = 0UL; @@ -724,7 +743,7 @@ public sealed class HeuristicPrimeTester if (admissible) { - candidate = CreateCandidate(value); + candidate = CreateCandidate(value, currentEnding, currentWheelResidue); return true; } } @@ -746,6 +765,8 @@ public sealed class HeuristicPrimeTester current = next; residueStepper.Advance(); + currentEnding = (byte)endingStepper.ComputeNext(next.Low); + currentWheelResidue = (ushort)wheelStepper.ComputeNext(next.Low); } public readonly ulong ProcessedCount => processedCount; @@ -754,9 +775,8 @@ public sealed class HeuristicPrimeTester public readonly bool Exhausted => !active; - private static HeuristicDivisorCandidate CreateCandidate(ulong value) + private static HeuristicDivisorCandidate CreateCandidate(ulong value, byte ending, ushort wheelResidue) { - byte ending = (byte)(value % 10UL); HeuristicDivisorGroup group; if (value == 3UL || value == 7UL || value == 11UL || value == 13UL) @@ -772,8 +792,7 @@ public sealed class HeuristicPrimeTester group = HeuristicDivisorGroup.GroupB; } - ushort residue = (ushort)(value % Wheel210); - return new HeuristicDivisorCandidate(value, group, ending, 0, residue); + return new HeuristicDivisorCandidate(value, group, ending, 0, wheelResidue); } } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 0c5fb2e..f6047c0 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -266,7 +266,7 @@ internal static partial class PrimeOrderCalculator int exponent = entries[i].Value; for (int iteration = 0; iteration < exponent; iteration++) { - if ((order % primeFactor) != UInt128.Zero) + if (order.ReduceCycleRemainder(primeFactor) != UInt128.Zero) { UInt128 candidate = order / primeFactor; if (Pow2ModWideCpu(candidate, prime, divisorData) == UInt128.One) @@ -339,14 +339,13 @@ internal static partial class PrimeOrderCalculator buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value)); - // TODO: Investigating using stepper to replace modulo operation for (int i = 0; i < length; i++) { UInt128 primeFactor = buffer[i].Value; int exponent = buffer[i].Exponent; for (int iteration = 0; iteration < exponent; iteration++) { - if ((order % primeFactor) == UInt128.Zero) + if (order.ReduceCycleRemainder(primeFactor) == UInt128.Zero) { UInt128 reduced = order / primeFactor; if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) @@ -405,7 +404,7 @@ internal static partial class PrimeOrderCalculator UInt128 reduced = order; for (int iteration = 0; iteration < span[i].Exponent; iteration++) { - if ((reduced % primeFactor) != UInt128.Zero) + if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) { break; } @@ -2583,7 +2582,7 @@ internal static partial class PrimeOrderCalculator // continue; // } - if ((remainingLocal % primeValue) != 0UL) + if (remainingLocal.ReduceCycleRemainder(primeValue) != 0UL) { continue; } @@ -2745,16 +2744,16 @@ internal static partial class PrimeOrderCalculator private static ulong AdvancePolynomialCpu(ulong x, ulong c, ulong modulus) { UInt128 value = (UInt128)x * x + c; - return (ulong)(value % modulus); + return (ulong)value.ReduceCycleRemainder(modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong AdvancePolynomialTwiceCpu(ulong x, ulong c, ulong modulus) { UInt128 value = (UInt128)x * x + c; - var remainder = (ulong)(value % modulus); + var remainder = (ulong)value.ReduceCycleRemainder(modulus); value = (UInt128)remainder * remainder + c; - return (ulong)(value % modulus); + return (ulong)value.ReduceCycleRemainder(modulus); } private static ulong BinaryGcd(ulong a, ulong b) @@ -2897,7 +2896,7 @@ internal static partial class PrimeOrderCalculator int exponent = entries[i].Value; for (int iteration = 0; iteration < exponent; iteration++) { - ulong remainder = order % primeFactor; + ulong remainder = order.ReduceCycleRemainder(primeFactor); if (remainder != 0UL) { break; @@ -2951,7 +2950,7 @@ internal static partial class PrimeOrderCalculator { int exponent = 1; ulong remaining = quotient; - while ((remaining % factor) == 0UL) + while (remaining.ReduceCycleRemainder(factor) == 0UL) { remaining /= factor; exponent++; @@ -3062,7 +3061,7 @@ internal static partial class PrimeOrderCalculator UInt128 reduced = candidate; for (int iteration = 0; iteration < span[i].Exponent; iteration++) { - if ((reduced % primeFactor) != UInt128.Zero) + if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) { break; } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index e8de8c4..53c0d8c 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -307,19 +307,19 @@ internal static partial class PrimeOrderCalculator { UInt128 primeFactor = span[i].Value; UInt128 reduced = order; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if ((reduced % primeFactor) != UInt128.Zero) - { - break; - } + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) + { + break; + } - reduced /= primeFactor; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } + reduced /= primeFactor; + if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) + { + return false; + } + } } return true; diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 7a0a504..c1f44ea 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -1,7 +1,5 @@ -using System; using System.Numerics; using System.Runtime.CompilerServices; -using ILGPU; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu.Accelerators; @@ -72,7 +70,7 @@ public sealed class PrimeTester } divisor = smallPrimeDivisors[i]; - uint deltaMod = delta >= divisor ? (uint)(delta % divisor) : (uint)delta; + uint deltaMod = delta.ReduceCycleRemainder(divisor); remainder = remainders[i] + deltaMod; if (remainder >= divisor) { @@ -105,7 +103,7 @@ public sealed class PrimeTester } divisor = smallPrimeDivisors[i]; - remainder = (uint)(n % divisor); + remainder = n.ReduceCycleRemainder(divisor); remainders[i] = remainder; if (remainder == 0UL) diff --git a/PerfectNumbers.Core/PrimeTesterByLastDigit.cs b/PerfectNumbers.Core/PrimeTesterByLastDigit.cs index ce5f9f1..06f4cd1 100644 --- a/PerfectNumbers.Core/PrimeTesterByLastDigit.cs +++ b/PerfectNumbers.Core/PrimeTesterByLastDigit.cs @@ -1,11 +1,29 @@ using System.Runtime.CompilerServices; -using ILGPU.Runtime; -using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core; public sealed class PrimeTesterByLastDigit { + [ThreadStatic] + private static BidirectionalCycleRemainderStepper _mod10Stepper; + + [ThreadStatic] + private static bool _mod10Initialized; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte GetMod10(ulong value) + { + if (!_mod10Initialized || _mod10Stepper.Modulus != 10UL) + { + _mod10Stepper = new BidirectionalCycleRemainderStepper(10UL); + _mod10Stepper.Initialize(value); + _mod10Initialized = true; + return (byte)_mod10Stepper.Step(value); + } + + return (byte)_mod10Stepper.Step(value); + } + public static bool IsPrimeCpu(ulong n) { // The below IFs never trigger in production code of EvenPerfectBitScanner on --mersenne=bydivisor path. @@ -27,7 +45,7 @@ public sealed class PrimeTesterByLastDigit uint[] smallPrimeDivisors; ulong[] smallPrimeDivisorsMul; - byte nMod10 = n.Mod10(); + byte nMod10 = GetMod10(n); switch (nMod10) { diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index 39b9f2a..8e27bef 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -360,9 +360,9 @@ namespace PerfectNumbers.Core private static bool _hasExponentRemainderStepperCpu; [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ExponentRemainderStepperCpu RentExponentStepperCpu(in MontgomeryDivisorData divisorData) + internal static ExponentRemainderStepperCpu RentExponentStepperCpu(in MontgomeryDivisorData divisorData, ulong cycleLength = 0UL) { - if (_hasExponentRemainderStepperCpu && _exponentRemainderStepperCpu.MatchesDivisor(divisorData)) + if (_hasExponentRemainderStepperCpu && _exponentRemainderStepperCpu.MatchesDivisor(divisorData, cycleLength)) { ExponentRemainderStepperCpu stepper = _exponentRemainderStepperCpu; _hasExponentRemainderStepperCpu = false; @@ -370,7 +370,7 @@ namespace PerfectNumbers.Core return stepper; } - return new ExponentRemainderStepperCpu(divisorData); + return new ExponentRemainderStepperCpu(divisorData, cycleLength); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 521c462..2968436 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -306,6 +306,53 @@ public static class UInt128Extensions return remainder.Pow2MontgomeryModWindowed(modulus); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 ReduceCycleRemainder(this UInt128 value, UInt128 modulus) + { + if (modulus == UInt128.Zero || value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + return value % modulus; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool ShouldUseSingleBit(UInt128 exponent) => (exponent >> 64) == UInt128.Zero && (ulong)exponent <= Pow2WindowFallbackThreshold; diff --git a/PerfectNumbers.Core/UIntExtensions.cs b/PerfectNumbers.Core/UIntExtensions.cs index c94b317..3fa26b6 100644 --- a/PerfectNumbers.Core/UIntExtensions.cs +++ b/PerfectNumbers.Core/UIntExtensions.cs @@ -13,9 +13,59 @@ public static class UIntExtensions public static uint Mod3(uint value) => value - (uint)((value * Mod3Mask) >> 33) * 3U; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod5(this uint value) => value - (uint)((value * Mod5Mask) >> 34) * 5U; + public static uint Mod5(this uint value) => value - (uint)((value * Mod5Mask) >> 34) * 5U; - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static uint Mod8(this uint value) => value & 7U; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static uint Mod8(this uint value) => value & 7U; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static uint ReduceCycleRemainder(this uint value, uint modulus) + { + if (modulus == 0U || value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + return value % modulus; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static uint ReduceCycleRemainder(this ulong value, uint modulus) + => (uint)value.ReduceCycleRemainder((ulong)modulus); } - diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index a24f589..931056a 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -7,8 +7,6 @@ namespace PerfectNumbers.Core; public static partial class ULongExtensions { private const ulong Pow2WindowFallbackThreshold = 32UL; - private const int Pow2WindowOddPowerCount = 1 << (PerfectNumberConstants.Pow2WindowSize - 1); - private static readonly UInt128 Pow2WindowedModulusThreshold = (UInt128)PerfectNumberConstants.MaxQForDivisorCycles; private static readonly byte[] BitStatsZeroCount = new byte[256]; @@ -1135,6 +1133,56 @@ public static partial class ULongExtensions return x1; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong ReduceCycleRemainder(this ulong value, ulong modulus) + { + // BigIntegerExtensions.Track(modulus); + // BigIntegerExtensions.DumpCycleRemainderStats(); + + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + value -= modulus; + if (value < modulus) + { + return value; + } + + return value % modulus; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static bool SharesFactorWithExponentMinusOneCpu(this ulong exponent) {