diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs index f480e12..ecf9022 100644 --- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs +++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs @@ -105,7 +105,7 @@ public static class MersennePrimeFactorTester } candidate = order / prime; - if (ModPow128(2, candidate, q, ct) == 1) + if (ModPowWithCancellation(UInt128Numbers.Two, candidate, q, ct) == UInt128.One) { order = candidate; } @@ -131,80 +131,42 @@ public static class MersennePrimeFactorTester } [MethodImpl(MethodImplOptions.AggressiveInlining)] - // TODO: Replace this duplicate powmod with ULongExtensions.ModPow64 so we inherit the optimized MulMod64 path that led the MulMod64Benchmarks on large operands. - private static ulong ModPow64(ulong value, ulong exponent, ulong modulus, CancellationToken ct) + private static UInt128 ModPowWithCancellation(UInt128 baseValue, UInt128 exponent, UInt128 modulus, CancellationToken ct) { - ulong result = 1UL; - ulong baseValue = value % modulus; - ulong exp = exponent; + // Keep this helper local to the factoring path so we can preserve the token-aware early outs + // without carrying a separate overload on the shared UInt128 extensions. + if (modulus <= UInt128.One) + { + return UInt128.Zero; + } - while (exp != 0UL) + UInt128 result = UInt128.One; + UInt128 baseResidue = baseValue % modulus; + + while (exponent != UInt128.Zero) { if (ct.IsCancellationRequested) { - return 0UL; + return UInt128.Zero; } - if ((exp & 1UL) != 0UL) + if ((exponent & UInt128.One) != UInt128.Zero) { - result = MulMod64(result, baseValue, modulus); + result = result.MulMod(baseResidue, modulus); } if (ct.IsCancellationRequested) { - return 0UL; + return UInt128.Zero; } - baseValue = MulMod64(baseValue, baseValue, modulus); - exp >>= 1; + baseResidue = baseResidue.MulMod(baseResidue, modulus); + exponent >>= 1; } return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong MulMod64(ulong a, ulong b, ulong mod) => (ulong)((UInt128)a * b % mod); // TODO: Route this through ULongExtensions.MulMod64 so the powmod above adopts the inline UInt128 implementation that dominated the MulMod64Benchmarks. - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - // TODO: Switch to UInt128Extensions.ModPow once its MulMod backend adopts the faster UInt128BuiltIn path highlighted in MulHighBenchmarks to avoid the BigInteger fallback. - private static UInt128 ModPow128(UInt128 value, UInt128 exponent, UInt128 modulus, CancellationToken ct) - { - UInt128 result = 1; - UInt128 baseValue = value % modulus; - UInt128 exp = exponent; - - while (exp != 0) - { - if (ct.IsCancellationRequested) - { - return 0; - } - - if ((exp & 1) != 0) - { - result = MulMod128(result, baseValue, modulus); - } - - if (ct.IsCancellationRequested) - { - return 0; - } - - baseValue = MulMod128(baseValue, baseValue, modulus); - exp >>= 1; - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 mod) - { - BigInteger product = (BigInteger)a * b; // TODO: Replace this BigInteger reduction with the forthcoming UInt128 intrinsic helper measured faster in Mul64Benchmarks and MulHighBenchmarks for huge operands. - product %= (BigInteger)mod; - return (UInt128)product; - } - private static readonly ConcurrentBag> ListPool = new(); private static readonly ConcurrentBag> DictPool = new(); private static readonly Dictionary FactorCache = new(); @@ -391,24 +353,19 @@ public static class MersennePrimeFactorTester return 0UL; } - x = (MulMod64(x, x, n) + c) % n; // TODO: Swap the `% n` with the Montgomery folding - // helper from the PollardRho benchmarks once the - // specialized reducer lands. + x = x.MulMod64(x, n).AddMod64(c, n); if (ct.IsCancellationRequested) { return 0UL; } - y = (MulMod64(y, y, n) + c) % n; // TODO: Use the same Montgomery folding helper here - // to keep the tortoise sequence on the optimized - // path. + y = y.MulMod64(y, n).AddMod64(c, n); if (ct.IsCancellationRequested) { return 0UL; } - y = (MulMod64(y, y, n) + c) % n; // TODO: Replace this modulo with the optimized helper - // so both steps share the fast reduction. + y = y.MulMod64(y, n).AddMod64(c, n); diff = x > y ? x - y : y - x; if (ct.IsCancellationRequested) { @@ -539,22 +496,19 @@ public static class MersennePrimeFactorTester return 0; } - x = (MulMod128(x, x, n) + c) % n; // TODO: Swap this modulo for the UInt128 Montgomery - // reducer once the optimized helper lands. + x = x.MulMod(x, n).AddMod(c, n); if (ct.IsCancellationRequested) { return 0; } - y = (MulMod128(y, y, n) + c) % n; // TODO: Use the same optimized reducer for the - // tortoise step. + y = y.MulMod(y, n).AddMod(c, n); if (ct.IsCancellationRequested) { return 0; } - y = (MulMod128(y, y, n) + c) % n; // TODO: Replace with the optimized reducer so both - // steps avoid BigInteger-based modulo. + y = y.MulMod(y, n).AddMod(c, n); diff = x > y ? x - y : y - x; if (ct.IsCancellationRequested) { diff --git a/PerfectNumbers.Core/ModResidueTracker.cs b/PerfectNumbers.Core/ModResidueTracker.cs index 9b0393c..0043699 100644 --- a/PerfectNumbers.Core/ModResidueTracker.cs +++ b/PerfectNumbers.Core/ModResidueTracker.cs @@ -44,7 +44,7 @@ public sealed class ModResidueTracker if (divisor == UInt128.Zero) { throw new ArgumentOutOfRangeException(nameof(divisor)); - } + } EnsureAt(number); @@ -58,15 +58,15 @@ public sealed class ModResidueTracker if (_hasLastAppended && idx == divisors.Count - 1 && divisor >= _lastAppended) { _lastAppended = divisor; - } + } else if (!_hasLastAppended) { _hasLastAppended = true; _lastAppended = divisor; - } + } return residue == UInt128.Zero; - } + } return _residues[idx] == UInt128.Zero; } @@ -90,7 +90,7 @@ public sealed class ModResidueTracker EnsureAt(number); // After EnsureAt we can keep _cursor (still valid monotonic pointer) if (_cursor > divisors.Count) _cursor = divisors.Count; - } + } UInt128 residue; // Fast path: append to tail when candidates are strictly increasing. @@ -105,21 +105,21 @@ public sealed class ModResidueTracker _cursor = divisors.Count; // position after last divisible = residue == UInt128.Zero; return true; - } + } // Merge pointer: advance until current >= divisor var cursor = _cursor; while (cursor < divisors.Count && divisors[cursor] < divisor) { cursor++; - } + } if (cursor < divisors.Count && divisors[cursor] == divisor) { divisible = residues[cursor] == UInt128.Zero; _cursor = cursor + 1; // position after matched return true; - } + } // Insert at cursor to keep order @@ -134,7 +134,7 @@ public sealed class ModResidueTracker { _hasLastAppended = true; _lastAppended = divisor; - } + } return true; } @@ -150,17 +150,17 @@ public sealed class ModResidueTracker // Initialize residues lazily upon first check for each divisor. _initialized = true; return; - } + } if (number == _currentNumber) { return; - } + } if (number < _currentNumber) { throw new InvalidOperationException("Number must be non-decreasing."); - } + } UInt128 delta = number - _currentNumber; List divisors = _divisors; @@ -169,7 +169,7 @@ public sealed class ModResidueTracker { _currentNumber = number; return; - } + } // Advance all tracked residues from _currentNumber to 'number'. // For Identity: r' = (r + delta) mod d @@ -206,17 +206,15 @@ public sealed class ModResidueTracker // TODO: Replace this PowMod128 call with the upcoming eight-bit window helper once the // ProcessEightBitWindows scalar implementation lands so residue updates benefit from the // ~2× win recorded in GpuPow2ModBenchmarks for large moduli. - value = PowMod128(UInt128Numbers.Two, delta, d); + UInt128 powDelta = UInt128Numbers.Two.ModPow(delta, d); // r' = pow2Delta*r + (pow2Delta - 1) mod d - value = MulMod128(value, residues[i], d) + value - UInt128.One; - if (value >= d) - { - value -= d; - } + value = powDelta.MulMod(residues[i], d); + value = value.AddMod(powDelta, d); + value = value.SubtractOneMod(d); - residues[i] = value; - } - } + residues[i] = value; + } + } _currentNumber = number; } @@ -231,12 +229,12 @@ public sealed class ModResidueTracker if (list[mid] < value) { lo = mid + 1; - } + } else { hi = mid; - } - } + } + } return lo; } @@ -264,70 +262,18 @@ public sealed class ModResidueTracker // instead of recomputing powmods for every divisor; the divisor-cycle benchmarks showed large wins once the // cached orders were used across scans, and a miss should trigger the configured device to compute the cycle // immediately without storing it back or requesting additional cache blocks. - UInt128 pow = PowMod128(UInt128Numbers.Two, number, divisor); + UInt128 pow = UInt128Numbers.Two.ModPow(number, divisor); // q divides M_p ⇔ 2^p ≡ 1 (mod q) AND ord_q(2) | p if (pow != UInt128.One) return (pow - UInt128.One) % divisor; ulong ord = divisor.CalculateOrder(); - if (ord == 0UL || divisor % ord != UInt128.Zero) + if (ord == 0UL || number % (UInt128)ord != UInt128.Zero) return UInt128.One; // does NOT divide return UInt128.Zero; // q really divides M_p } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static UInt128 PowMod128(UInt128 baseValue, UInt128 exponent, UInt128 modulus) - { - if (modulus == UInt128.One) - { - return UInt128.Zero; - } - - UInt128 result = UInt128.One; - UInt128 value = baseValue % modulus; - - while (exponent != UInt128.Zero) - { - if ((exponent & UInt128.One) != UInt128.Zero) - { - // TODO: Route this through the planned UInt128 windowed powmod so the inner MulMod adopts the - // UInt128BuiltIn path that dominated the MulHighBenchmarks on wide operands instead of the - // current square-and-multiply loop. - result = MulMod128(result, value, modulus); - } - - exponent >>= 1; - if (exponent != UInt128.Zero) - { - value = MulMod128(value, value, modulus); - } - } - - return result; - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 modulus) - { - // Same approach as MulMod in MersenneNumberTester: double-and-add to avoid overflow. - // TODO: Replace this double-and-add fallback with the UInt128 intrinsic-backed multiplier once the - // MulHighBenchmarks-guided implementation lands; the intrinsic path measured dozens of times faster for - // dense 128-bit workloads. - UInt128 result = UInt128.Zero; - a %= modulus; - - while (b != UInt128.Zero) - { - if ((b & UInt128.One) != UInt128.Zero) - { - result += result < modulus ? a : a - modulus; - } - - a = a < modulus ? a << 1 : (a << 1) - modulus; - b >>= 1; - } - - return result; - } + } diff --git a/PerfectNumbers.Core/UInt128Extensions.cs b/PerfectNumbers.Core/UInt128Extensions.cs index 42d5c2a..91e596a 100644 --- a/PerfectNumbers.Core/UInt128Extensions.cs +++ b/PerfectNumbers.Core/UInt128Extensions.cs @@ -1,5 +1,6 @@ using System.Numerics; using System.Runtime.CompilerServices; +using System.Threading; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -630,6 +631,57 @@ public static class UInt128Extensions return left; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 AddMod(this UInt128 value, UInt128 addend, UInt128 modulus) + { + // All current CPU callers (Pollard Rho and the Mersenne residue tracker) only supply moduli + // of at least three, so leave this guard documented but disabled. + // if (modulus <= UInt128.One) + // { + // return UInt128.Zero; + // } + + // Both operands arrive pre-reduced on these paths, letting us skip the redundant folds the old + // helper performed. + // value %= modulus; + // addend %= modulus; + // Pollard Rho keeps its polynomial constant inside [1, modulus - 1], and the Mersenne pow-delta + // term never hits zero for odd moduli, leaving this branch inactive outside targeted tests. + // if (addend == UInt128.Zero) + // { + // return value; + // } + + UInt128 threshold = modulus - addend; + if (value >= threshold) + { + return value - threshold; + } + + UInt128 sum = value + addend; + // The sum stays below the modulus under the threshold guard above, so no wraparound occurs. + return sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 SubtractOneMod(this UInt128 value, UInt128 modulus) + { + // Only the Mersenne residue tracker calls this helper, and it never observes moduli below three, + // so leave this guard documented but inactive. + // if (modulus <= UInt128.One) + // { + // return UInt128.Zero; + // } + + if (value == UInt128.Zero) + { + // Wrapping occurs when the prior AddMod call produced one, so step back to modulus - 1. + return modulus - UInt128.One; + } + + return value - UInt128.One; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 Pow(this UInt128 value, ulong exponent) { diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 33c7585..fdd9f07 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -455,13 +455,73 @@ public static partial class ULongExtensions return result; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong AddMod64(this ulong value, ulong addend, ulong modulus) + { + // Every CPU caller (Pollard Rho and the ModPow64 ladder) supplies moduli of at least three, so leave this guard + // documented but disabled. + // if (modulus <= 1UL) + // { + // return 0UL; + // } + + // Both flows keep their operands below the modulus ahead of time, letting us skip the re-fold the old helper performed. + UInt128 sum = (UInt128)value + addend; + if (sum < modulus) + { + return (ulong)sum; + } + + sum -= modulus; + // After one subtraction the sum always falls below the modulus on this path, making the old guard redundant. + // if (sum >= modulus) + // { + // sum -= modulus; + // } + + return (ulong)sum; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong FoldMontgomery(ulong value, ulong modulus) + { + // Pollard Rho and ModPow64 both feed moduli of at least three, so retain the documentation and skip the runtime guard. + // if (modulus <= 1UL) + // { + // return 0UL; + // } + + if (value < modulus) + { + return value; + } + + ulong folded = 0UL; + int leadingZeros = BitOperations.LeadingZeroCount(value); + int bitIndex = 63 - leadingZeros; + for (; bitIndex >= 0; bitIndex--) + { + folded = folded.AddMod64(folded, modulus); + if (((value >> bitIndex) & 1UL) != 0UL) + { + folded = folded.AddMod64(1UL, modulus); + } + } + + return folded; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus) { + // Pollard Rho and the Mersenne testers never pass moduli below two, so keep the guard commented for clarity. + // if (modulus <= 1UL) + // { + // return 0UL; + // } + ulong result = 1UL; - // TODO: Replace this `%` with the Montgomery folding helper highlighted in MulMod64Benchmarks so the - // modular exponentiation avoids the slow integer division before the ladder even starts. - value %= modulus; + value = FoldMontgomery(value, modulus); while (exponent != 0UL) { @@ -477,10 +537,13 @@ public static partial class ULongExtensions return result; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - // TODO: Replace this fallback with the UInt128 Montgomery helper measured fastest in - // MulMod64Benchmarks so CPU callers stop paying for triple modulo operations. - public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + // Reuse the direct UInt128 reduction path that dominated MulMod64Benchmarks to avoid + // the redundant operand modulo operations in the older fallback. + public static ulong MulMod64(this ulong a, ulong b, ulong modulus) + { + return (ulong)(((UInt128)a * b) % modulus); + } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime)