diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs index b794d5b..2d6731e 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs @@ -59,7 +59,7 @@ public class Pow2MontgomeryModBenchmarks for (int i = 0; i < SampleCount; i++) { - checksum ^= exponents[i].Pow2MontgomeryMod(divisors[i]); + checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false); } return checksum; @@ -221,14 +221,20 @@ public class Pow2MontgomeryModBenchmarks { if (modulus <= 1UL || (modulus & 1UL) == 0UL) { - return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); + return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL); } + ulong nPrime = ComputeMontgomeryNPrime(modulus); + ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus); + ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus); + ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); + return new MontgomeryDivisorData( modulus, - ComputeMontgomeryNPrime(modulus), - ComputeMontgomeryResidue(1UL, modulus), - ComputeMontgomeryResidue(2UL, modulus)); + nPrime, + montgomeryOne, + montgomeryTwo, + montgomeryTwoSquared); } private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); @@ -472,3 +478,4 @@ public class Pow2MontgomeryModBenchmarks return (ulong)(product % modulus); } } + diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs index f6ab991..de6810d 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModCycleComputationBenchmarks.cs @@ -1,4 +1,5 @@ using System; +using System.Buffers.Binary; using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; using PerfectNumbers.Core; @@ -21,20 +22,25 @@ namespace EvenPerfectBitScanner.Benchmarks; public class Pow2MontgomeryModCycleComputationBenchmarks { private const int SampleCount = 256; - + private const ulong OneShiftedLeft63 = 1UL << 63; + private const long OneShiftedLeft60 = 1L << 60; private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount]; private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount]; + // private readonly MontgomeryDivisorData[] _veryLargeDivisors = new MontgomeryDivisorData[SampleCount]; private readonly ulong[] _smallExponents = new ulong[SampleCount]; private readonly ulong[] _largeExponents = new ulong[SampleCount]; + // private readonly ulong[] _veryLargeExponents = new ulong[SampleCount]; private readonly ulong[] _smallCycles = new ulong[SampleCount]; private readonly ulong[] _largeCycles = new ulong[SampleCount]; + // private readonly ulong[] _veryLargeCycles = new ulong[SampleCount]; private readonly Random _random = new(113); public enum InputScale { Small, - Large + Large, + // VeryLarge } [ParamsAllValues] @@ -45,6 +51,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks { for (int i = 0; i < SampleCount; i++) { + Console.WriteLine($"Generating sample {i}"); ulong smallModulus = NextSmallOddModulus(); _smallDivisors[i] = CreateMontgomeryDivisorData(smallModulus); _smallExponents[i] = NextSmallExponent(); @@ -54,6 +61,12 @@ public class Pow2MontgomeryModCycleComputationBenchmarks _largeDivisors[i] = CreateMontgomeryDivisorData(largeModulus); _largeExponents[i] = NextLargeExponent(); _largeCycles[i] = largeCycle; + + // ulong veryLargeModulus = NextVeryLargeOddModulus(); + // _veryLargeDivisors[i] = CreateMontgomeryDivisorData(veryLargeModulus); + // _veryLargeExponents[i] = NextVeryLargeExponent(); + // Console.WriteLine($"Calculating cycle length {veryLargeModulus}"); + // _veryLargeCycles[i] = MersenneDivisorCycles.CalculateCycleLengthGpu(veryLargeModulus); } } @@ -69,7 +82,7 @@ public class Pow2MontgomeryModCycleComputationBenchmarks for (int i = 0; i < SampleCount; i++) { - checksum ^= exponents[i].Pow2MontgomeryMod(divisors[i]); + checksum ^= exponents[i].Pow2MontgomeryModWindowed(divisors[i], keepMontgomery: false); } return checksum; @@ -114,23 +127,39 @@ public class Pow2MontgomeryModCycleComputationBenchmarks private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles) { - if (Scale == InputScale.Small) + switch (Scale) { - exponents = _smallExponents; - divisors = _smallDivisors; - cycles = _smallCycles; - } - else - { - exponents = _largeExponents; - divisors = _largeDivisors; - cycles = _largeCycles; + case InputScale.Small: + exponents = _smallExponents; + divisors = _smallDivisors; + cycles = _smallCycles; + break; + case InputScale.Large: + exponents = _largeExponents; + divisors = _largeDivisors; + cycles = _largeCycles; + break; + // case InputScale.VeryLarge: + // exponents = _veryLargeExponents; + // divisors = _veryLargeDivisors; + // cycles = _veryLargeCycles; + // break; + default: + throw new ArgumentOutOfRangeException(nameof(Scale), Scale, null); } } private ulong NextSmallExponent() => (ulong)_random.NextInt64(1L, 1_000_000L); - private ulong NextLargeExponent() => (ulong)_random.NextInt64(1L << 60, long.MaxValue); + private ulong NextLargeExponent() => (ulong)_random.NextInt64(OneShiftedLeft60, long.MaxValue); + + private ulong NextVeryLargeExponent() + { + Span buffer = stackalloc byte[8]; + _random.NextBytes(buffer); + ulong value = checked(BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63); + return value; + } private ulong NextSmallOddModulus() { @@ -145,18 +174,32 @@ public class Pow2MontgomeryModCycleComputationBenchmarks return (modulus, (ulong)bitLength); } + private ulong NextVeryLargeOddModulus() + { + Span buffer = stackalloc byte[8]; + _random.NextBytes(buffer); + ulong value = BinaryPrimitives.ReadUInt64LittleEndian(buffer) | OneShiftedLeft63 | 1UL; + return value; + } + private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) { if (modulus <= 1UL || (modulus & 1UL) == 0UL) { - return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); + return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL, 0UL); } + ulong nPrime = ComputeMontgomeryNPrime(modulus); + ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus); + ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus); + ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); + return new MontgomeryDivisorData( modulus, - ComputeMontgomeryNPrime(modulus), - ComputeMontgomeryResidue(1UL, modulus), - ComputeMontgomeryResidue(2UL, modulus)); + nPrime, + montgomeryOne, + montgomeryTwo, + montgomeryTwoSquared); } private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index c411c60..c8da4b7 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -128,7 +128,7 @@ internal static class Program mersenneOption = arg.AsSpan(arg.IndexOf('=') + 1); if (mersenneOption.Equals("pow2mod", StringComparison.OrdinalIgnoreCase)) { - kernelType = GpuKernelType.Pow2Mod; // TODO: Switch this mode to the ProcessEightBitWindows kernel once Pow2Minus1Mod adopts the benchmarked faster windowed ladder. + kernelType = GpuKernelType.Pow2Mod; // Dispatch to the ProcessEightBitWindows pow2mod kernel. } else if (mersenneOption.Equals("lucas", StringComparison.OrdinalIgnoreCase)) { @@ -485,7 +485,7 @@ internal static class Program { MersenneTesters = new ThreadLocal(() => { - // TODO: Swap the underlying pow2mod kernels to ProcessEightBitWindows once Pow2Minus1Mod migrates to the windowed helper highlighted in GpuPow2ModBenchmarks. + // ProcessEightBitWindows windowed pow2 ladder is the default kernel. var tester = new MersenneNumberTester( useIncremental: !useLucas, useOrderCache: false, @@ -1362,7 +1362,7 @@ internal static class Program searchedMersenne = true; if (_useByDivisorMode) { - // TODO: Route by-divisor scans through the ProcessEightBitWindows pow2mod kernel once the windowed Pow2Minus1Mod helper replaces the slower single-bit ladder. + // Windowed pow2mod kernel handles by-divisor scans. return _byDivisorTester!.IsPrime(p, out detailedCheck); } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 9370406..476fb4a 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,5 +1,6 @@ using System; using System.Collections.Concurrent; +using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; @@ -173,6 +174,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } + Dictionary? factorCache = null; ulong divisor = 3UL; while (divisor <= allowedMax) { @@ -183,7 +185,19 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv // TODO: Hoist MontgomeryDivisorData acquisition into the divisor-cycle cache so we reuse the // staged ProcessEightBitWindows-ready metadata instead of reloading the slower Montgomery // structs for every divisor. - ulong divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + ulong divisorCycle; + if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + { + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + } + else + { + factorCache ??= new Dictionary(8); + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, prime, factorCache, out divisorCycle) || divisorCycle == 0UL) + { + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + } + } if (divisorCycle != 0UL) { byte hit = CheckDivisor(prime, divisorCycle, divisorData); @@ -220,12 +234,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) { - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL || divisorCycle == 0UL) - { - return 0; - } - ulong residue = prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData); return residue == 1UL ? (byte)1 : (byte)0; } @@ -308,12 +316,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } - if (divisor <= 1UL || (divisor & 1UL) == 0UL) - { - hits.Clear(); - return; - } - if (divisorCycle == 0UL) { divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepper.cs index 7d2f38c..39fb669 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepper.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepper.cs @@ -47,7 +47,7 @@ internal struct ExponentRemainderStepper // TODO: Route these state resets through the ProcessEightBitWindows helper once the scalar // Pow2MontgomeryMod implementation adopts it so delta stepping inherits the benchmarked // 2× gains for large exponents. - _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); + _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); _previousExponent = exponent; _hasState = true; return ReduceCurrent(); @@ -59,7 +59,7 @@ internal struct ExponentRemainderStepper // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks. // TODO: Replace this per-delta powmod with the upcoming windowed ladder so incremental // updates stop paying the single-bit cost highlighted in GpuPow2ModBenchmarks. - ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor); + ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _previousExponent = exponent; return ReduceCurrent(); @@ -77,7 +77,7 @@ internal struct ExponentRemainderStepper { // TODO: Switch this reload to the shared windowed pow2 helper once available so CPU // residue checks align with the optimized ProcessEightBitWindows timings. - _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); + _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); _previousExponent = exponent; _hasState = true; return _currentMontgomery == _montgomeryOne; @@ -88,7 +88,7 @@ internal struct ExponentRemainderStepper // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor. // TODO: Use the windowed delta pow2 helper here as well to avoid the single-bit ladder that // currently lags behind the benchmarked implementation for huge divisor cycles. - ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor); + ulong multiplier = delta.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _previousExponent = exponent; return _currentMontgomery == _montgomeryOne; @@ -123,7 +123,7 @@ internal struct ExponentRemainderStepper { // TODO: Replace this fallback path with the upcoming ProcessEightBitWindows helper so // fresh Montgomery states also benefit from the faster pow2 ladder measured on CPUs. - _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); + _currentMontgomery = exponent.Pow2MontgomeryModWindowed(_divisor, keepMontgomery: true); _previousExponent = exponent; _hasState = true; isUnity = _currentMontgomery == _montgomeryOne; diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index 4d54771..e213e93 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -1,3 +1,4 @@ +using System; using System.Numerics; using System.Runtime.CompilerServices; @@ -15,6 +16,10 @@ public struct GpuUInt128 : IComparable, IEquatable private const ulong NativeModuloChunkMask = (1UL << NativeModuloChunkBits) - 1UL; private static readonly ulong[] NativeModuloBitMasks = CreateNativeModuloBitMasks(); + private const int Pow2WindowSizeBits = 8; + private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1); + private const int Pow2WindowMaxExponent = (1 << Pow2WindowSizeBits) - 1; + public static readonly GpuUInt128 Zero = new(); public static readonly GpuUInt128 One = new(1UL); @@ -256,122 +261,24 @@ public struct GpuUInt128 : IComparable, IEquatable return Zero; } - // TODO: Replace this single-bit ladder with the ProcessEightBitWindows strategy measured - // in GpuPow2ModBenchmarks once the GPU kernels can share the windowed helper. The windowed - // version held 21.5 µs for 33-bit moduli versus 51.0 µs here when scanning large divisors. - GpuUInt128 result = Zero; - ulong i = 0UL; - while (i < exponent) + GpuUInt128 pow = Pow2Mod(exponent, modulus); + if (pow.IsZero) { - // result = (result * 2) % modulus - result += result; - if (result.CompareTo(modulus) >= 0) - { - result -= modulus; - } - - // result = (result + 1) % modulus - result += One; - if (result.CompareTo(modulus) >= 0) - { - result -= modulus; - } - - i++; + pow = modulus; } - return result; + pow.Sub(One); + return pow; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void Pow2Minus1ModBatch(GpuUInt128 modulus, ReadOnlySpan exponents, Span results) { - GpuUInt128 r0 = Zero, r1 = Zero, r2 = Zero, r3 = Zero, r4 = Zero, r5 = Zero, r6 = Zero, r7 = Zero; - ulong e0 = 0UL, e1 = 0UL, e2 = 0UL, e3 = 0UL, e4 = 0UL, e5 = 0UL, e6 = 0UL, e7 = 0UL; - bool a0 = false, a1 = false, a2 = false, a3 = false, a4 = false, a5 = false, a6 = false, a7 = false; - - // TODO: Mirror the eight-bit window batching once the scalar helper above switches; the - // current bit-serial loop lags by ~3× on 8k–33-bit moduli compared to the windowed path. - int len = exponents.Length; - if (len > 0) { e0 = exponents[0]; a0 = true; } - if (len > 1) { e1 = exponents[1]; a1 = true; } - if (len > 2) { e2 = exponents[2]; a2 = true; } - if (len > 3) { e3 = exponents[3]; a3 = true; } - if (len > 4) { e4 = exponents[4]; a4 = true; } - if (len > 5) { e5 = exponents[5]; a5 = true; } - if (len > 6) { e6 = exponents[6]; a6 = true; } - if (len > 7) { e7 = exponents[7]; a7 = true; } - - ulong max = 0UL; - if (a0 && e0 > max) max = e0; - if (a1 && e1 > max) max = e1; - if (a2 && e2 > max) max = e2; - if (a3 && e3 > max) max = e3; - if (a4 && e4 > max) max = e4; - if (a5 && e5 > max) max = e5; - if (a6 && e6 > max) max = e6; - if (a7 && e7 > max) max = e7; - - for (ulong i = 0UL; i < max; i++) + int count = exponents.Length; + for (int i = 0; i < count; i++) { - if (a0) - { - r0 += r0; if (r0.CompareTo(modulus) >= 0) r0 -= modulus; - r0 += One; if (r0.CompareTo(modulus) >= 0) r0 -= modulus; - if (i + 1UL == e0) a0 = false; - } - if (a1) - { - r1 += r1; if (r1.CompareTo(modulus) >= 0) r1 -= modulus; - r1 += One; if (r1.CompareTo(modulus) >= 0) r1 -= modulus; - if (i + 1UL == e1) a1 = false; - } - if (a2) - { - r2 += r2; if (r2.CompareTo(modulus) >= 0) r2 -= modulus; - r2 += One; if (r2.CompareTo(modulus) >= 0) r2 -= modulus; - if (i + 1UL == e2) a2 = false; - } - if (a3) - { - r3 += r3; if (r3.CompareTo(modulus) >= 0) r3 -= modulus; - r3 += One; if (r3.CompareTo(modulus) >= 0) r3 -= modulus; - if (i + 1UL == e3) a3 = false; - } - if (a4) - { - r4 += r4; if (r4.CompareTo(modulus) >= 0) r4 -= modulus; - r4 += One; if (r4.CompareTo(modulus) >= 0) r4 -= modulus; - if (i + 1UL == e4) a4 = false; - } - if (a5) - { - r5 += r5; if (r5.CompareTo(modulus) >= 0) r5 -= modulus; - r5 += One; if (r5.CompareTo(modulus) >= 0) r5 -= modulus; - if (i + 1UL == e5) a5 = false; - } - if (a6) - { - r6 += r6; if (r6.CompareTo(modulus) >= 0) r6 -= modulus; - r6 += One; if (r6.CompareTo(modulus) >= 0) r6 -= modulus; - if (i + 1UL == e6) a6 = false; - } - if (a7) - { - r7 += r7; if (r7.CompareTo(modulus) >= 0) r7 -= modulus; - r7 += One; if (r7.CompareTo(modulus) >= 0) r7 -= modulus; - if (i + 1UL == e7) a7 = false; - } + results[i] = Pow2Minus1Mod(exponents[i], modulus); } - - if (len > 0) results[0] = r0; - if (len > 1) results[1] = r1; - if (len > 2) results[2] = r2; - if (len > 3) results[3] = r3; - if (len > 4) results[4] = r4; - if (len > 5) results[5] = r5; - if (len > 6) results[6] = r6; - if (len > 7) results[7] = r7; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -483,7 +390,6 @@ public struct GpuUInt128 : IComparable, IEquatable } } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static GpuUInt128 Pow2Mod(ulong exponent, GpuUInt128 modulus) { @@ -492,36 +398,71 @@ public struct GpuUInt128 : IComparable, IEquatable return new(); } - const int WindowSizeBits = 8; - const ulong WindowMask = (1UL << WindowSizeBits) - 1UL; - - GpuUInt128 result = new(1UL); - GpuUInt128 baseVal = new(2UL); - - ulong e = exponent; - while (e != 0UL) + if (exponent == 0UL) { - int remainingBits = BitOperations.Log2(e) + 1; - int bitsToProcess = remainingBits >= WindowSizeBits ? WindowSizeBits : remainingBits; - ulong chunk = e & (WindowMask >> (WindowSizeBits - bitsToProcess)); + return new GpuUInt128(1UL); + } - for (int bitIndex = 0; bitIndex < bitsToProcess; bitIndex++) + Span oddPowers = stackalloc GpuUInt128[Pow2WindowOddPowerCount]; + GpuUInt128 power = new GpuUInt128(1UL); + int oddIndex = 0; + + for (int bit = 1; bit <= Pow2WindowMaxExponent; bit++) + { + power += power; + if (power.CompareTo(modulus) >= 0) { - if ((chunk & 1UL) != 0UL) - { - result.MulMod(baseVal, modulus); - } - - baseVal.MulMod(baseVal, modulus); - chunk >>= 1; + power.Sub(modulus); } - e >>= bitsToProcess; + if ((bit & 1) != 0) + { + oddPowers[oddIndex++] = power; + } + } + + GpuUInt128 result = new GpuUInt128(1UL); + int bitLength = BitOperations.Log2(exponent) + 1; + int index = bitLength - 1; + + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result.SquareMod(modulus); + index--; + continue; + } + + int windowStart = index - Pow2WindowSizeBits + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + int windowBitCount = index - windowStart + 1; + for (int square = 0; square < windowBitCount; square++) + { + result.SquareMod(modulus); + } + + ulong mask = (1UL << windowBitCount) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + + GpuUInt128 factor = oddPowers[tableIndex]; + result.MulMod(factor, modulus); + + index = windowStart - 1; } return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static int TrailingZeroCount(GpuUInt128 value) { @@ -1258,7 +1199,6 @@ public struct GpuUInt128 : IComparable, IEquatable p3 = hHH + c2; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong MulMod64(ulong a, ulong b, ulong modulus) { diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 6d451df..10690e6 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -1,6 +1,7 @@ using System; using System.Buffers; using System.Collections.Concurrent; +using System.Collections.Generic; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Numerics; @@ -371,7 +372,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulong exponent = exponents[index]; - hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; // TODO: Replace Pow2MontgomeryMod with the ProcessEightBitWindows helper once available so GPU by-divisor scans reuse the benchmarked fast pow2 ladder. + hits[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -384,6 +385,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private MemoryBuffer1D _resultsBuffer; private ulong[] _hostBuffer; private int[]? _positionBuffer; + private Dictionary? _factorCache; private int _capacity; private bool _disposed; @@ -403,6 +405,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal void Reset() { _disposed = false; + _factorCache?.Clear(); } public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) @@ -428,11 +431,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv if (divisorCycle == 0UL) { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); - if (divisorCycle == 0UL) + Dictionary cache = _factorCache ??= new Dictionary(16); + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, primes[0], cache, out divisorCycle) || divisorCycle == 0UL) { - hits.Clear(); - return; + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + if (divisorCycle == 0UL) + { + hits.Clear(); + return; + } } } @@ -569,6 +576,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } _disposed = true; + _factorCache?.Clear(); _owner._sessionPool.Add(this); } } @@ -584,7 +592,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong exponent = exponents[index]; - results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); // TODO: Swap this Montgomery-domain pow2 with the windowed helper slated to replace the slower ladder in Pow2Minus1Mod benchmarks. + results[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: true); } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index c2e0d27..f8cbd4c 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -1,9 +1,14 @@ using System.Buffers; +using System.Buffers.Binary; +using System.Collections.Generic; +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Security.Cryptography; using System.Text; +using System.Threading; using PerfectNumbers.Core.Gpu; using ILGPU; using ILGPU.Runtime; -using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; @@ -146,6 +151,342 @@ public class MersenneDivisorCycles return order; } + public static bool TryCalculateCycleLengthForExponent( + ulong divisor, + ulong exponent, + Dictionary? factorCache, + out ulong cycleLength) + { + cycleLength = 0UL; + + if ((divisor & (divisor - 1UL)) == 0UL) + { + cycleLength = 1UL; + return true; + } + + if (divisor <= 3UL || exponent <= 1UL) + { + cycleLength = CalculateCycleLength(divisor); + return true; + } + + ulong phi = divisor - 1UL; + if (phi == 0UL) + { + return false; + } + + int twoCount = BitOperations.TrailingZeroCount(phi); + ulong reducedPhi = phi >> twoCount; + if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL) + { + return false; + } + + ulong k = reducedPhi / exponent; + if (k == 0UL) + { + return false; + } + + Dictionary factorCounts = new(8); + if (twoCount > 0) + { + factorCounts[2UL] = twoCount; + } + + if (!AccumulateFactors(exponent, factorCache, factorCounts) || !AccumulateFactors(k, factorCache, factorCounts)) + { + return false; + } + + cycleLength = ReduceOrder(divisor, phi, factorCounts); + return true; + } + + private static bool AccumulateFactors( + ulong value, + Dictionary? cache, + Dictionary counts) + { + if (value <= 1UL) + { + return true; + } + + if (!TryGetFactorization(value, cache, out FactorCacheEntry factorization)) + { + return false; + } + + ulong[] primes = factorization.Primes; + byte[] exponents = factorization.Exponents; + int length = factorization.Count; + + for (int i = 0; i < length; i++) + { + AddFactor(counts, primes[i], exponents[i]); + } + + return true; + } + + private static bool TryGetFactorization( + ulong value, + Dictionary? cache, + out FactorCacheEntry factorization) + { + if (cache is not null && cache.TryGetValue(value, out factorization)) + { + return true; + } + + Dictionary scratch = new(8); + if (!TryFactorIntoCountsInternal(value, scratch)) + { + factorization = default!; + return false; + } + + int count = scratch.Count; + ulong[] primes = new ulong[count]; + byte[] exponents = new byte[count]; + int index = 0; + foreach (KeyValuePair entry in scratch) + { + primes[index] = entry.Key; + exponents[index] = checked((byte)entry.Value); + index++; + } + + Array.Sort(primes, exponents); + factorization = new FactorCacheEntry(primes, exponents, count); + + if (cache is not null) + { + cache[value] = factorization; + } + + return true; + } + + private static bool TryFactorIntoCountsInternal(ulong value, Dictionary counts) + { + if (value <= 1UL) + { + return true; + } + + ulong remaining = value; + uint[] smallPrimes = PrimesGenerator.SmallPrimes; + ulong[] smallPrimesSquared = PrimesGenerator.SmallPrimesPow2; + int smallLength = smallPrimes.Length; + + for (int i = 0; i < smallLength; i++) + { + ulong prime = smallPrimes[i]; + if (smallPrimesSquared[i] > remaining) + { + break; + } + + while ((remaining % prime) == 0UL) + { + AddFactor(counts, prime); + remaining /= prime; + } + } + + if (remaining == 1UL) + { + return true; + } + + if (PrimeTester.IsPrimeInternal(remaining, CancellationToken.None)) + { + AddFactor(counts, remaining); + return true; + } + + ulong factor = PollardRho64(remaining); + if (factor == 0UL || factor == remaining) + { + return false; + } + + ulong quotient = remaining / factor; + return TryFactorIntoCountsInternal(factor, counts) && TryFactorIntoCountsInternal(quotient, counts); + } + + private static void AddFactor(Dictionary counts, ulong factor) + { + if (counts.TryGetValue(factor, out int existing)) + { + counts[factor] = existing + 1; + } + else + { + counts[factor] = 1; + } + } + + private static void AddFactor(Dictionary counts, ulong factor, int multiplicity) + { + if (counts.TryGetValue(factor, out int existing)) + { + counts[factor] = existing + multiplicity; + } + else + { + counts[factor] = multiplicity; + } + } + + private static ulong ReduceOrder(ulong divisor, ulong initialOrder, Dictionary factorCounts) + { + if (factorCounts.Count == 0) + { + return initialOrder; + } + + MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); + int count = factorCounts.Count; + ulong[] primes = ArrayPool.Shared.Rent(count); + try + { + int index = 0; + foreach (KeyValuePair entry in factorCounts) + { + primes[index] = entry.Key; + index++; + } + + Array.Sort(primes, 0, count); + + ulong order = initialOrder; + for (int i = 0; i < count; i++) + { + ulong prime = primes[i]; + int multiplicity = factorCounts[prime]; + + for (int iteration = 0; iteration < multiplicity; iteration++) + { + if (order % prime != 0UL) + { + break; + } + + ulong candidate = order / prime; + if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + { + order = candidate; + continue; + } + + break; + } + } + + return order; + } + finally + { + ArrayPool.Shared.Return(primes, clearArray: false); + } + } + + private static ulong PollardRho64(ulong n) + { + if ((n & 1UL) == 0UL) + { + return 2UL; + } + + while (true) + { + ulong c = (NextRandomUInt64() % (n - 1UL)) + 1UL; + ulong x = (NextRandomUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + x = AdvancePolynomial(x, c, n); + y = AdvancePolynomial(y, c, n); + y = AdvancePolynomial(y, c, n); + + ulong diff = x > y ? x - y : y - x; + d = BinaryGcd(diff, n); + } + + if (d != n) + { + return d; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) + { + UInt128 value = (UInt128)x * x + c; + return (ulong)(value % modulus); + } + + private static ulong BinaryGcd(ulong a, ulong b) + { + if (a == 0UL) + { + return b; + } + + if (b == 0UL) + { + return a; + } + + int shift = BitOperations.TrailingZeroCount(a | b); + a >>= BitOperations.TrailingZeroCount(a); + + while (true) + { + b >>= BitOperations.TrailingZeroCount(b); + if (a > b) + { + (a, b) = (b, a); + } + + b -= a; + if (b == 0UL) + { + return a << shift; + } + } + } + + private static ulong NextRandomUInt64() + { + Span buffer = stackalloc byte[8]; + RandomNumberGenerator.Fill(buffer); + return BinaryPrimitives.ReadUInt64LittleEndian(buffer); + } + + public readonly struct FactorCacheEntry + { + public FactorCacheEntry(ulong[] primes, byte[] exponents, int count) + { + Primes = primes; + Exponents = exponents; + Count = count; + } + public ulong[] Primes { get; } + + public byte[] Exponents { get; } + + public int Count { get; } + } + public static (long nextPosition, long completeCount) FindLast(string path) { using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false); diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index 1a2291f..806ac5b 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -237,6 +237,7 @@ public static class MersenneNumberDivisorByDivisorTester Span hitsSpan = default; int hitIndex = 0; int index; + Dictionary? factorCache = null; ulong divisorCycle = 0UL; try @@ -254,7 +255,6 @@ public static class MersenneNumberDivisorByDivisorTester divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted, ref localDivisorCursor, ref localDivisorsRemaining); if (divisor != 0UL) { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); activeCount = BuildPrimeBuffer( divisor, @@ -281,6 +281,19 @@ public static class MersenneNumberDivisorByDivisorTester { continue; } + divisorCycle = 0UL; + if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + { + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + } + else + { + factorCache ??= new Dictionary(16); + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, primeBuffer[0], factorCache, out divisorCycle) || divisorCycle == 0UL) + { + divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); + } + } hitsSpan = hitsBuffer.AsSpan(0, activeCount); hitsSpan.Clear(); diff --git a/PerfectNumbers.Core/MontgomeryDivisorData.cs b/PerfectNumbers.Core/MontgomeryDivisorData.cs index fa73db7..0027afb 100644 --- a/PerfectNumbers.Core/MontgomeryDivisorData.cs +++ b/PerfectNumbers.Core/MontgomeryDivisorData.cs @@ -4,42 +4,42 @@ using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; -public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo) +public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo, ulong montgomeryTwoSquared) { - public ulong Modulus { get; } = modulus; + public readonly ulong Modulus = modulus; - public ulong NPrime { get; } = nPrime; + public readonly ulong NPrime = nPrime; - public ulong MontgomeryOne { get; } = montgomeryOne; + public readonly ulong MontgomeryOne = montgomeryOne; - public ulong MontgomeryTwo { get; } = montgomeryTwo; + public readonly ulong MontgomeryTwo = montgomeryTwo; + public readonly ulong MontgomeryTwoSquared = montgomeryTwoSquared; } internal static class MontgomeryDivisorDataCache { private static readonly ConcurrentDictionary Cache = new(); - public static MontgomeryDivisorData Get(ulong modulus) - { - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); - } - - return Cache.GetOrAdd(modulus, static m => Create(m)); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static MontgomeryDivisorData Get(ulong modulus) => Cache.GetOrAdd(modulus, static m => Create(m)); private static MontgomeryDivisorData Create(ulong modulus) { + ulong nPrime = ComputeMontgomeryNPrime(modulus); + ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus); + ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus); + ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); + return new MontgomeryDivisorData( modulus, - ComputeMontgomeryNPrime(modulus), - ComputeMontgomeryResidue(1UL, modulus), - ComputeMontgomeryResidue(2UL, modulus)); + nPrime, + montgomeryOne, + montgomeryTwo, + montgomeryTwoSquared); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); + private static ulong ComputeMontgomeryResidue(UInt128 value, ulong modulus) => (ulong)(value % modulus); private static ulong ComputeMontgomeryNPrime(ulong modulus) { diff --git a/PerfectNumbers.Core/UInt128Numbers.cs b/PerfectNumbers.Core/UInt128Numbers.cs index b2bde3d..bccaa5b 100644 --- a/PerfectNumbers.Core/UInt128Numbers.cs +++ b/PerfectNumbers.Core/UInt128Numbers.cs @@ -6,6 +6,9 @@ namespace PerfectNumbers.Core public static readonly UInt128 Four = 4UL; public static readonly UInt128 SixtyFour = 64UL; public static readonly UInt128 OneHundredTwentySeven = 127UL; + public static readonly UInt128 OneShiftedLeft64 = UInt128.One << 64; + public static readonly UInt128 OneShiftedLeft64x2 = Two * OneShiftedLeft64; + // TODO: Populate this table from the UInt128 residue benchmarks so GPU kernels can pull precomputed constants // without re-deriving them at runtime. // TODO: Extend the table with cached Mod3/Mod5 folding constants so both CPU and GPU residue helpers can reuse diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 64dc3d4..e3bac57 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -1,3 +1,4 @@ +using System; using System.Numerics; using System.Runtime.CompilerServices; using PerfectNumbers.Core.Gpu; @@ -6,6 +7,8 @@ namespace PerfectNumbers.Core; public static class ULongExtensions { + private const int Pow2WindowSize = 8; + private const ulong Pow2WindowFallbackThreshold = 32UL; public static ulong CalculateOrder(this ulong q) { if (q <= 2UL) @@ -293,25 +296,85 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryMod(this ulong exponent, in MontgomeryDivisorData divisor) + public static ulong Pow2MontgomeryModWindowed(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) { ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) + if (exponent == 0UL) { - return 0UL; + return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; } + if (exponent <= Pow2WindowFallbackThreshold) + { + return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + } + + int bitLength = BitOperations.Log2(exponent) + 1; + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + ulong result = divisor.MontgomeryOne; + ulong nPrime = divisor.NPrime; + Span oddPowers = stackalloc ulong[oddPowerCount]; + InitializeMontgomeryOddPowers(divisor, modulus, nPrime, oddPowers); + + int index = bitLength - 1; + while (index >= 0) + { + if (((exponent >> index) & 1UL) == 0UL) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (((exponent >> windowStart) & 1UL) == 0UL) + { + windowStart++; + } + + int windowLength = index - windowStart + 1; + for (int square = 0; square < windowLength; square++) + { + result = result.MontgomeryMultiply(result, modulus, nPrime); + } + + ulong mask = (1UL << windowLength) - 1UL; + ulong windowValue = (exponent >> windowStart) & mask; + int tableIndex = (int)((windowValue - 1UL) >> 1); + ulong multiplier = oddPowers[tableIndex]; + result = result.MontgomeryMultiply(multiplier, modulus, nPrime); + + index = windowStart - 1; + } + + if (keepMontgomery) + { + return result; + } + + return result.MontgomeryMultiply(1UL, modulus, nPrime); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong Pow2MontgomeryModSingleBit(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { + ulong modulus = divisor.Modulus; + ulong nPrime = divisor.NPrime; ulong result = divisor.MontgomeryOne; ulong baseVal = divisor.MontgomeryTwo; - ulong nPrime = divisor.NPrime; ulong remainingExponent = exponent; - while (remainingExponent > 0UL) + while (remainingExponent != 0UL) { if ((remainingExponent & 1UL) != 0UL) { - // TODO: Switch to the ProcessEightBitWindows helper once the scalar implementation ships so CPU callers - // inherit the same 2× speedup that the benchmarked windowed ladder delivered for large exponents. result = result.MontgomeryMultiply(baseVal, modulus, nPrime); } @@ -324,68 +387,68 @@ public static class ULongExtensions baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime); } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModMontgomery(this ulong exponent, in MontgomeryDivisorData divisor) + private static int GetWindowSize(int bitLength) { - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) + if (bitLength <= 6) { - return 0UL; + return Math.Max(bitLength, 1); } - ulong result = divisor.MontgomeryOne; - ulong baseVal = divisor.MontgomeryTwo; - ulong remainingExponent = exponent; - - while (remainingExponent > 0UL) + if (bitLength <= 23) { - if ((remainingExponent & 1UL) != 0UL) - { - // TODO: Route this Montgomery-domain variant through the same windowed pow2 helper once available so delta - // stepping matches the ProcessEightBitWindows gains highlighted by GpuPow2ModBenchmarks. - result = result.MontgomeryMultiply(baseVal, modulus, divisor.NPrime); - } - - remainingExponent >>= 1; - if (remainingExponent == 0UL) - { - break; - } - - baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, divisor.NPrime); + return 4; } - return result; + if (bitLength <= 79) + { + return 5; + } + + if (bitLength <= 239) + { + return 6; + } + + if (bitLength <= 671) + { + return 7; + } + + return Pow2WindowSize; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void InitializeMontgomeryOddPowers(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span oddPowers) + { + oddPowers[0] = divisor.MontgomeryTwo; + if (oddPowers.Length == 1) + { + return; + } + + ulong square = divisor.MontgomeryTwoSquared; + for (int i = 1; i < oddPowers.Length; i++) + { + oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(square, modulus, nPrime); + } } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) { - if (cycleLength == 0UL) - { - return exponent.Pow2MontgomeryMod(divisor); - } - - // TODO: Swap this modulo with the divisor-cycle remainder helper from Pow2MontgomeryModCycleComputationBenchmarks so we - // advance the exponent using cached cycle deltas instead of paying for a fresh `%` on every query. ulong rotationCount = exponent % cycleLength; - return rotationCount.Pow2MontgomeryModFromCycleRemainder(divisor); + return Pow2MontgomeryModWindowed(rotationCount, divisor, keepMontgomery: false); } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor) { - if (reducedExponent == 0UL) - { - return 1UL; - } - - return reducedExponent.Pow2MontgomeryMod(divisor); + return Pow2MontgomeryModWindowed(reducedExponent, divisor, keepMontgomery: false); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 PowMod(this ulong exponent, UInt128 modulus) {