diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 91e596a..2f52e08 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -65,19 +65,16 @@ public static class UInt128Extensions [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static UInt128 Pow2MontgomeryModWindowed(this UInt128 exponent, UInt128 modulus) + public static UInt128 Pow2MontgomeryModWindowed(this ulong exponent, UInt128 modulus) { - // TODO: Prototype a UInt128-native Pow2MontgomeryModWindowed path that matches the GPU implementation once a faster modular multiplication helper is available, so we can revisit removing the GPU struct dependency without regressing benchmarks. - if (modulus == UInt128.One) - { - return UInt128.Zero; - } - - if (exponent == UInt128.Zero) - { - return UInt128.One % modulus; - } - + // EvenPerfectBitScanner keeps CPU exponents above Pow2WindowFallbackThreshold once p >= 138,000,000, + // so production never routes the tiny-exponent single-bit fallback through this helper. Leave the + // guard commented out for targeted microbenchmarks. + // This helper executes only on CPU fallbacks (PrimeOrderCalculator.Shared and PrimeOrderGpuHeuristics + // call it after the GPU heuristics decline), but we still use GpuUInt128 because its mutable Montgomery + // operations beat the pure-UInt128 prototype by roughly 12% in Pow2MontgomeryModBenchmarks. + // TODO: Swap in a CPU-native Montgomery ladder once it surpasses the shared struct so GPU and CPU tuning + // can diverge without compromising either path. GpuUInt128 modulusGpu = new(modulus); GpuUInt128 baseValue = new(2UL); if (baseValue.CompareTo(modulusGpu) >= 0) @@ -85,11 +82,11 @@ public static class UInt128Extensions baseValue.Sub(modulusGpu); } - if (ShouldUseSingleBit(exponent)) - { - GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponent, modulusGpu, baseValue); - return (UInt128)singleBitResult; - } + // if (exponent <= Pow2WindowFallbackThreshold) + // { + // GpuUInt128 tinyResult = Pow2MontgomeryModSingleBit((UInt128)exponent, modulusGpu, baseValue); + // return (UInt128)tinyResult; + // } GpuUInt128 exponentGpu = new(exponent); int bitLength = exponentGpu.GetBitLength(); @@ -139,6 +136,176 @@ public static class UInt128Extensions return (UInt128)result; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 Pow2MontgomeryModWindowed(this UInt128 exponent, UInt128 modulus) + { + // TODO: Prototype a UInt128-native Pow2MontgomeryModWindowed path that matches the GPU implementation once a faster + // modular multiplication helper is available, so we can replace the shared GpuUInt128 usage on CPU fallbacks and let + // the GPU kernels continue using their dedicated ladder without cross-contaminating optimizations. + // All GPU kernels already route through Pow2MontgomeryModWindowedGpu, and Pow2MontgomeryModBenchmarks confirmed this + // windowed ladder stays the fastest CPU option with the current Montgomery helpers. + if (modulus == UInt128.One) + { + return UInt128.Zero; + } + + if (exponent == UInt128.Zero) + { + return UInt128.One % modulus; + } + + // This helper still executes on CPU-only call sites, so we lean on the shared GpuUInt128 ladder until + // the TODO above lands and gives the CPU path its own Montgomery struct. + GpuUInt128 modulusGpu = new(modulus); + GpuUInt128 baseValue = new(2UL); + if (baseValue.CompareTo(modulusGpu) >= 0) + { + baseValue.Sub(modulusGpu); + } + + // EvenPerfectBitScanner CPU exponents exceed Pow2WindowFallbackThreshold once p >= 138,000,000, + // so the single-bit fallback never triggers during production scans. Leave it commented out for targeted benchmarks. + // if (ShouldUseSingleBit(exponent)) + // { + // GpuUInt128 singleBitResult = Pow2MontgomeryModSingleBit(exponent, modulusGpu, baseValue); + // return (UInt128)singleBitResult; + // } + + GpuUInt128 exponentGpu = new(exponent); + int bitLength = exponentGpu.GetBitLength(); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + Span oddPowers = stackalloc GpuUInt128[oddPowerCount]; + InitializeOddPowers(baseValue, modulusGpu, oddPowers); + + GpuUInt128 result = GpuUInt128.One; + int index = bitLength - 1; + + while (index >= 0) + { + if (!IsBitSet(exponentGpu, index)) + { + result.MulMod(result, modulusGpu); + index--; + continue; + } + + int windowStart = index - windowSize + 1; + if (windowStart < 0) + { + windowStart = 0; + } + + while (!IsBitSet(exponentGpu, windowStart)) + { + windowStart++; + } + + int windowBitCount = index - windowStart + 1; + for (int square = 0; square < windowBitCount; square++) + { + result.MulMod(result, modulusGpu); + } + + ulong windowValue = ExtractWindowValue(exponentGpu, windowStart, windowBitCount); + int tableIndex = (int)((windowValue - 1UL) >> 1); + GpuUInt128 factor = oddPowers[tableIndex]; + result.MulMod(factor, modulusGpu); + + index = windowStart - 1; + } + + return (UInt128)result; + } + + /// + /// Uses the known cycle length to reduce the exponent and reuse cached remainders before + /// running the windowed pow2 ladder. + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 Pow2MontgomeryModWindowedWithCycle(this UInt128 exponent, UInt128 modulus, ulong cycleLength, UInt128? precomputedRemainder = null) + { + // The EvenPerfectBitScanner CPU path receives cycle lengths from MersenneDivisorCycles, which never + // produces a zero-length order for admissible divisors. Keep the fallback documented for tests. + // if (cycleLength == 0UL) + // { + // return exponent.Pow2MontgomeryModWindowed(modulus); + // } + + UInt128 remainder; + if (precomputedRemainder.HasValue) + { + // PowModWithCycle caches the reduced rotation count before reaching this helper, so production + // scans always reuse the precomputed remainder instead of recomputing the modulus here. + remainder = precomputedRemainder.Value; + } + else if (cycleLength != 0UL && exponent >= cycleLength) + { + if ((exponent >> 64) == UInt128.Zero) + { + remainder = (UInt128)((ulong)exponent % cycleLength); + } + else + { + remainder = exponent % (UInt128)cycleLength; + } + } + else + { + remainder = exponent; + } + + // When the cycle divides the exponent, 2^exponent ≡ 1 (mod modulus) and the caller expects an + // immediate success signal. + if (remainder == UInt128.Zero) + { + return UInt128.One; + } + + return remainder.Pow2MontgomeryModWindowed(modulus); + } + + /// + /// Uses the known cycle length to reduce the exponent and reuse cached remainders before + /// running the windowed pow2 ladder. + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 Pow2MontgomeryModWindowedWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength, UInt128? precomputedRemainder = null) + { + // Divisor-cycle lookups never hand the EvenPerfectBitScanner CPU flow a zero cycle length. Keep the + // fallback documented for targeted benchmarks and tests. + // if (cycleLength == UInt128.Zero) + // { + // return exponent.Pow2MontgomeryModWindowed(modulus); + // } + + UInt128 remainder; + if (precomputedRemainder.HasValue) + { + // PowModWithCycle forwards the cached rotation count, so production scans always reuse the + // supplied remainder instead of recomputing exponent % cycleLength. + remainder = precomputedRemainder.Value; + } + else if (cycleLength != UInt128.Zero && exponent >= cycleLength) + { + remainder = exponent % cycleLength; + } + else + { + remainder = exponent; + } + + // A zero remainder means the divisor order divides the exponent, so 2^exponent ≡ 1 (mod modulus) + // and the sieves rely on the early exit. + if (remainder == UInt128.Zero) + { + return UInt128.One; + } + + return remainder.Pow2MontgomeryModWindowed(modulus); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool ShouldUseSingleBit(UInt128 exponent) => (exponent >> 64) == UInt128.Zero && (ulong)exponent <= Pow2WindowFallbackThreshold; diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index fdd9f07..08af224 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -10,6 +10,8 @@ public static partial class ULongExtensions private const ulong Pow2WindowFallbackThreshold = 32UL; private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); + private static readonly UInt128 Pow2WindowedModulusThreshold = (UInt128)PerfectNumberConstants.MaxQForDivisorCycles; + private static readonly byte[] BitStatsZeroCount = new byte[256]; private static readonly byte[] BitStatsPrefixZero = new byte[256]; private static readonly byte[] BitStatsSuffixZero = new byte[256]; @@ -591,224 +593,257 @@ public static partial class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static UInt128 PowMod(this ulong exponent, UInt128 modulus) - { - UInt128 result = UInt128.One; - ulong exponentLoopIndex = 0UL; + public static UInt128 PowMod(this ulong exponent, UInt128 modulus) + { + // EvenPerfectBitScanner starts CPU scans at p >= 138,000,000 (Program.cs), + // so production never routes exponent == 0 here. Keep the guard commented out + // for targeted micro-benchmarks. + // if (exponent == 0UL) + // { + // return UInt128.One; + // } - // TODO: Port this scalar PowMod fallback to the ProcessEightBitWindows helper so CPU callers get the - // eight-bit window wins measured against the classic square-and-subtract implementation. - // Return 1 because 2^0 = 1 - if (exponent == 0UL) - return result; + // MersenneResidueAutomaton seeds q = 2 * p * k + 1 (ModuloAutomata.cs), + // making modulus >= 2 * p + 1 > 276,000,000 and odd on the CPU path. The + // modulus == 1 guard is therefore unreachable during EvenPerfectBitScanner + // runs. + // if (modulus == UInt128.One) + // { + // return UInt128.Zero; + // } - // Any number mod 1 is 0 - if (modulus == UInt128.One) - return UInt128.Zero; + // Production exponents stay above 64 bits and modulus never drops below 4 + // for admissible divisors, so the classic doubling fallback does not + // trigger in scans. + // if (exponent < 64 || modulus < 4) + // { + // return Pow2ByDoubling(exponent, modulus); + // } - // For small exponents, do classic method - if (exponent < 64 || modulus < 4) - { - for (; exponentLoopIndex < exponent; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } + // Every admissible divisor q is odd (2 * p * k + 1), preventing the + // power-of-two mask path from activating in production. + // if ((modulus & (modulus - 1)) == 0) + // { + // UInt128 masked = UInt128.One << (int)(exponent & 127); + // return masked & (modulus - 1); + // } - return result; - } + if (!ShouldUseWindowedPow2(modulus)) + { + return Pow2ByDoubling(exponent, modulus); + } - // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness - if ((modulus & (modulus - 1)) == 0) - { - result = UInt128.One << (int)(exponent & 127); - return result & (modulus - 1); - } + return exponent.Pow2MontgomeryModWindowed(modulus); + } - // Reusing exponentLoopIndex to iterate again over the exponent range for the general-case accumulation. - exponentLoopIndex = 0UL; - // Reusing result after resetting it for the general modulus accumulation phase. - result = UInt128.One; - for (; exponentLoopIndex < exponent; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } - - return result; - } /// /// Computes 2^exponent mod modulus using a known cycle length. /// - public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, ulong cycleLength) - { - UInt128 result = UInt128.One; - ulong exponentLoopIndex = 0UL; + + public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, ulong cycleLength) + { + // EvenPerfectBitScanner keeps exponent >= 138,000,000 on CPU (Program.cs), + // so the zero-exponent guard stays disabled outside tests. + // if (exponent == 0UL) + // { + // return UInt128.One; + // } - // TODO: Wire this cycle-aware overload into the ProcessEightBitWindows helper so the reduced exponent path - // inherits the faster windowed pow2 routine highlighted in the Pow2Montgomery benchmarks. - // Return 1 because 2^0 = 1 - if (exponent == 0UL) - return result; + // The divisor generator produces odd q = 2 * p * k + 1 (ModuloAutomata.cs), + // so modulus never reaches 1 in production. + // if (modulus == UInt128.One) + // { + // return UInt128.Zero; + // } - // Any number mod 1 is 0 - if (modulus == UInt128.One) - return UInt128.Zero; + // Every admissible divisor is odd, so the power-of-two shortcut cannot + // trigger while scanning real candidates. + // if ((modulus & (modulus - 1)) == 0) + // { + // UInt128 masked = UInt128.One << (int)(exponent & 127); + // return masked & (modulus - 1); + // } - // For small exponents, do classic method - if (exponent < 64 || modulus < 4) - { - for (; exponentLoopIndex < exponent; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } + bool hasCycle = cycleLength != 0UL; + // TODO: Replace this modulo with the cached cycle remainder produced by the divisor-cycle cache so PowModWithCycle avoids + // repeated `%` work, matching the ProcessEightBitWindows wins captured in Pow2MontgomeryModCycleComputationBenchmarks. + ulong rotationCount = hasCycle && exponent >= cycleLength ? exponent % cycleLength : exponent; - return result; - } + // Keep this fast path: valid divisors satisfy cycle | exponent, producing + // a zero remainder that means 2^exponent ≡ 1 (mod q). + if (hasCycle && rotationCount == 0UL) + { + return UInt128.One; + } + + // FastDiv64 can return tiny divisors (see MersenneNumberIncrementalCpuTester), + // so the small-rotation fallback still fires on the CPU path. + if (modulus < 4 || rotationCount < 64) + { + return Pow2ByDoubling(rotationCount, modulus); + } + + if (ShouldUseWindowedPow2(modulus)) + { + UInt128 exponent128 = exponent; + UInt128 rotationCount128 = rotationCount; + return exponent128.Pow2MontgomeryModWindowedWithCycle(modulus, cycleLength, rotationCount128); + } + + return Pow2ByDoubling(rotationCount, modulus); + } - // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness - if ((modulus & (modulus - 1)) == 0) - { - result = UInt128.One << (int)(exponent & 127); - return result & (modulus - 1); - } - // Reusing exponentLoopIndex to iterate over the rotation count for the cycle-aware path. - exponentLoopIndex = 0UL; - // Reusing result after resetting it for the rotation-accumulation pass. - result = UInt128.One; - // TODO: Replace this modulo with the cached cycle remainder produced by the divisor-cycle cache so PowModWithCycle avoids - // repeated `%` work, matching the ProcessEightBitWindows wins captured in Pow2MontgomeryModCycleComputationBenchmarks. - ulong rotationCount = exponent % cycleLength; - for (; exponentLoopIndex < rotationCount; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } - return result; - } /// /// Computes 2^exponent mod modulus using a known cycle length. /// - public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, UInt128 cycleLength) - { - UInt128 result = UInt128.One; - ulong exponentLoopIndex = 0UL; + + public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, UInt128 cycleLength) + { + // CPU scans only feed exponents >= 138,000,000 (Program.cs), so keep the + // zero-exponent shortcut commented out for tests. + // if (exponent == 0UL) + // { + // return UInt128.One; + // } - // TODO: Replace this UInt128-cycle overload with the ProcessEightBitWindows helper so large-exponent CPU scans - // reuse the faster windowed pow2 ladder instead of the manual rotation loop measured to lag behind in benchmarks. - // Return 1 because 2^0 = 1 - if (exponent == UInt128.Zero) - return result; + // Modulus originates from q = 2 * p * k + 1 (ModuloAutomata.cs), therefore + // modulus == 1 never occurs outside dedicated benchmarks. + // if (modulus == UInt128.One) + // { + // return UInt128.Zero; + // } - // Any number mod 1 is 0 - if (modulus == UInt128.One) - return UInt128.Zero; + // Admissible q are odd, rendering the power-of-two shortcut unreachable. + // if ((modulus & (modulus - UInt128.One)) == UInt128.Zero) + // { + // UInt128 masked = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven); + // return masked & (modulus - UInt128.One); + // } - // For small exponents, do classic method - if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four) - { - for (; exponentLoopIndex < exponent; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } + bool hasCycle = cycleLength != UInt128.Zero; + UInt128 exponent128 = exponent; + // TODO: Swap this modulo with the upcoming UInt128 cycle remainder helper so large-exponent scans reuse cached + // reductions instead of recomputing `%` for every lookup, as demonstrated in Pow2MontgomeryModCycleComputationBenchmarks. + UInt128 rotationCount = hasCycle && exponent128 >= cycleLength ? exponent128 % cycleLength : exponent128; - return result; - } + if (hasCycle && rotationCount == UInt128.Zero) + { + return UInt128.One; + } - // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness - if ((modulus & (modulus - UInt128.One)) == UInt128.Zero) - { - result = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven); - return result & (modulus - 1); - } + if (modulus < UInt128Numbers.Four || rotationCount < UInt128Numbers.SixtyFour) + { + return Pow2ByDoubling(rotationCount, modulus); + } + + if (ShouldUseWindowedPow2(modulus)) + { + return exponent128.Pow2MontgomeryModWindowedWithCycle(modulus, cycleLength, rotationCount); + } + + return Pow2ByDoubling(rotationCount, modulus); + } - // Reusing result after resetting it for the rotation-driven accumulation phase. - result = UInt128.One; - // TODO: Swap this modulo with the upcoming UInt128 cycle remainder helper so large-exponent scans reuse cached - // reductions instead of recomputing `%` for every lookup, as demonstrated in Pow2MontgomeryModCycleComputationBenchmarks. - UInt128 rotationCount = exponent % cycleLength; - UInt128 rotationIndex = UInt128.Zero; - while (rotationIndex < rotationCount) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - rotationIndex += UInt128.One; - } - return result; - } /// /// Computes 2^exponent mod modulus using a known cycle length. /// - public static UInt128 PowModWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength) - { - UInt128 one = UInt128.One, - result = one, - zero = UInt128.Zero; - ulong exponentLoopIndex = 0UL; + + public static UInt128 PowModWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength) + { + UInt128 one = UInt128.One; + UInt128 zero = UInt128.Zero; - // TODO: Migrate this UInt128 exponent overload to ProcessEightBitWindows so the large-cycle reductions drop the - // slow manual loop that underperforms the windowed pow2 helper in the Pow2 benchmark suite. - // Return 1 because 2^0 = 1 - if (exponent == zero) - return result; + // EvenPerfectBitScanner drives this path with exponents derived from p, phi, + // or divisor quotients, none of which are zero on the production scan. + // if (exponent == zero) + // { + // return one; + // } - // Any number mod 1 is 0 - if (modulus == one) - return zero; + // Modulus again comes from q = 2 * p * k + 1, so modulus == 1 is unreachable. + // if (modulus == one) + // { + // return zero; + // } - // For small exponents, do classic method - if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four) - { - for (; exponentLoopIndex < exponent; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } + // Candidate divisors stay odd, therefore the power-of-two shortcut stays off. + // if ((modulus & (modulus - one)) == zero) + // { + // UInt128 masked = one << (int)(exponent & UInt128Numbers.OneHundredTwentySeven); + // return masked & (modulus - one); + // } - return result; - } + bool hasCycle = cycleLength != zero; + // TODO: Swap this modulo with the shared UInt128 cycle remainder helper once available so CRT powmods reuse cached + // reductions in the windowed ladder, avoiding the `%` cost highlighted in Pow2MontgomeryModCycleComputationBenchmarks. + UInt128 rotationCount = hasCycle && exponent >= cycleLength ? exponent % cycleLength : exponent; - // Special case: if modulus is a power of two, use bitmasking for efficiency and correctness - if ((modulus & (modulus - one)) == zero) - { - result = one << (int)(exponent & UInt128Numbers.OneHundredTwentySeven); - return result & (modulus - 1); - } + if (hasCycle && rotationCount == zero) + { + return one; + } - // Reusing result after resetting it for the rotation-driven accumulation phase. - result = one; - // TODO: Swap this modulo with the shared UInt128 cycle remainder helper once available so CRT powmods reuse cached - // reductions in the windowed ladder, avoiding the `%` cost highlighted in Pow2MontgomeryModCycleComputationBenchmarks. - UInt128 rotationCount = exponent % cycleLength; + if (modulus < UInt128Numbers.Four || rotationCount < UInt128Numbers.SixtyFour) + { + return Pow2ByDoubling(rotationCount, modulus); + } - // We're reusing "zero" as rotation index for just a little better performance - while (zero < rotationCount) - { - result <<= 1; - if (result >= modulus) - result -= modulus; + if (ShouldUseWindowedPow2(modulus)) + { + return exponent.Pow2MontgomeryModWindowedWithCycle(modulus, cycleLength, rotationCount); + } - zero += one; - } + return Pow2ByDoubling(rotationCount, modulus); + } - return result; - } + + + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool ShouldUseWindowedPow2(UInt128 modulus) => modulus > Pow2WindowedModulusThreshold; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static UInt128 Pow2ByDoubling(ulong exponent, UInt128 modulus) + { + UInt128 result = UInt128.One; + for (ulong i = 0; i < exponent; i++) + { + result <<= 1; + while (result >= modulus) + { + result -= modulus; + } + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static UInt128 Pow2ByDoubling(UInt128 exponent, UInt128 modulus) + { + UInt128 result = UInt128.One; + UInt128 index = UInt128.Zero; + + while (index < exponent) + { + result <<= 1; + while (result >= modulus) + { + result -= modulus; + } + + index += UInt128.One; + } + + return result; + } /// /// Computes 2^exponent mod modulus using iterative CRT composition from mod 10 up to modulus.