diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index 27fd86d..ce66960 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -257,10 +257,11 @@ public static class MersenneNumberDivisorByDivisorTester } else { + var scheduler = UnboundedTaskScheduler.Instance; ParallelOptions options = new() { MaxDegreeOfParallelism = workerCount, - TaskScheduler = UnboundedTaskScheduler.Instance + TaskScheduler = scheduler }; Parallel.ForEach(filteredPrimes, options, ProcessPrime); diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 11403d7..9da99ec 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -1,10 +1,8 @@ using System.Buffers; -using System.Collections.Generic; using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; -using System.Threading; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; @@ -12,17 +10,11 @@ namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - internal readonly struct PendingEntry - { - public readonly ulong Value; - public readonly bool KnownComposite; - - public PendingEntry(ulong value, bool knownComposite) - { - Value = value; - KnownComposite = knownComposite; - } - } + internal readonly struct PendingEntry(ulong value, bool knownComposite) + { + public readonly ulong Value = value; + public readonly bool KnownComposite = knownComposite; + } private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { @@ -70,8 +62,7 @@ internal static partial class PrimeOrderCalculator ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); - PartialFactorResult? orderFactors; - if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, out orderFactors)) + if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) { return candidateOrder; } @@ -104,6 +95,10 @@ internal static partial class PrimeOrderCalculator // continue; // } + // TODO: Can we optimize it thanks to knowing the factors, differences between factors, factors and / or phi + // being monotonically increasing etc. For example, if phi increases and factor increases, what does it tell us + // about the value of reduced? Or, maybe wee can use residue stepping to calculate the residue instead of calling + // Pow2equalsOneCpu? Maybe we don't need to calculate it from scratch? ulong reduced = phi / factor; if (Pow2EqualsOneCpu(reduced, prime, divisorData)) { @@ -237,8 +232,10 @@ internal static partial class PrimeOrderCalculator { ulong primeFactor = span[i].Value; ulong reduced = order; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) + int exponent = span[i].Exponent; + for (int iteration = 0; iteration < exponent; iteration++) { + // TODO: Benchmark and / or implement DivRem like solution or residue stepper? if ((reduced % primeFactor) != 0UL) { break; @@ -299,7 +296,7 @@ internal static partial class PrimeOrderCalculator orderFactors = extended; } - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks << 2; List candidates = ThreadStaticPools.RentUlongList(capacity); FactorEntry[] factorArray = orderFactors.Factors!; // DebugLog("Building candidates list"); @@ -669,7 +666,7 @@ internal static partial class PrimeOrderCalculator // GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); // if (status == GpuPow2ModStatus.Success) // { - return remainder == 1UL; + return remainder == 1UL; // } } @@ -755,14 +752,14 @@ internal static partial class PrimeOrderCalculator // exponentSlots.Clear(); counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(FactorSlotCount, 8)); - counts.Clear(); + // counts.Clear(); // bool gpuPopulated = TryPopulateSmallPrimeFactorsGpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining); // if (!gpuPopulated) // { // throw new InvalidOperationException($"GPU didn't populate factors for {value}"); - // counts.Clear(); - remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); + // counts.Clear(); + remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); // } int dictionaryCount = counts.Count; @@ -844,66 +841,68 @@ internal static partial class PrimeOrderCalculator } ulong cofactor = 1UL; - bool cofactorContainsComposite = false; - int pendingCount = pending.Count; - int index = 0; - for (; index < pendingCount; index++) - { - PendingEntry entry = pending[index]; - ulong composite = entry.Value; + bool cofactorContainsComposite = false; + int pendingCount = pending.Count; + int index = 0; + for (; index < pendingCount; index++) + { + PendingEntry entry = pending[index]; + ulong composite = entry.Value; - if (entry.KnownComposite) - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - continue; - } + if (entry.KnownComposite) + { + cofactor = checked(cofactor * composite); + cofactorContainsComposite = true; + continue; + } - bool isPrime = GetOrComputePrimality(composite); + bool isPrime = GetOrComputePrimality(composite); - if (isPrime) - { - AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); - } - // composite is never smaller on the execution path - // else if (composite > 1UL) - else - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - } - } + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + } + // composite is never smaller on the execution path + // else if (composite > 1UL) + else + { + cofactor = checked(cofactor * composite); + cofactorContainsComposite = true; + } + } + + bool cofactorIsPrime; + if (cofactor <= 1UL) + { + cofactorIsPrime = false; + } + else if (cofactorContainsComposite) + { + cofactorIsPrime = false; + } + else + { + cofactorIsPrime = GetOrComputePrimality(cofactor); + } - bool cofactorIsPrime; - if (cofactor <= 1UL) - { - cofactorIsPrime = false; - } - else if (cofactorContainsComposite) - { - cofactorIsPrime = false; - } - else - { - cofactorIsPrime = GetOrComputePrimality(cofactor); - } ArrayPool pool = ThreadStaticPools.FactorEntryPool; bool temp; if (useDictionary) { - if ((counts is null || counts.Count == 0) && cofactor == value) + temp = counts is null || counts.Count == 0; + if (temp && cofactor == value) { result = PartialFactorResult.Rent(null, cofactor, false, 0, cofactorIsPrime); goto ReturnResult; } - if (counts is null || counts.Count == 0) + if (temp) { result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0, cofactorIsPrime); goto ReturnResult; } - FactorEntry[] factors = pool.Rent(counts.Count); + FactorEntry[] factors = pool.Rent(counts!.Count); index = 0; foreach (KeyValuePair entry in counts) { @@ -1137,8 +1136,7 @@ internal static partial class PrimeOrderCalculator break; } - ulong primeValue = primeCandidate; - if ((remaining % primeValue) != 0UL) + if ((remaining % primeCandidate) != 0UL) { continue; } @@ -1146,12 +1144,13 @@ internal static partial class PrimeOrderCalculator int exponent = 0; do { - remaining /= primeValue; + // TODO: Implement DivRem like solution or residue stepper to re-use previous iteration results in new calculation. Add / use method, likely to ULongExtensions + remaining /= primeCandidate; exponent++; } - while ((remaining % primeValue) == 0UL); + while ((remaining % primeCandidate) == 0UL); - counts[primeValue] = exponent; + counts[primeCandidate] = exponent; } return remaining; @@ -1398,70 +1397,70 @@ internal static partial class PrimeOrderCalculator } private static void FactorCompletely(ulong value, Dictionary counts) - { - FactorCompletely(value, counts, knownComposite: false); - } + { + FactorCompletely(value, counts, knownComposite: false); + } - private static void FactorCompletely(ulong value, Dictionary counts, bool knownComposite) - { - if (value <= 1UL) - { - return; - } + private static void FactorCompletely(ulong value, Dictionary counts, bool knownComposite) + { + if (value <= 1UL) + { + return; + } - if (!knownComposite && GetOrComputePrimality(value)) - { - AddFactor(counts, value, 1); - return; - } + if (!knownComposite && GetOrComputePrimality(value)) + { + AddFactor(counts, value, 1); + return; + } - ulong factor = PollardRhoStrict(value); - ulong quotient = value / factor; + ulong factor = PollardRhoStrict(value); + ulong quotient = value / factor; - bool factorIsPrime = GetOrComputePrimality(factor); - if (factorIsPrime) - { - int exponent = 1; - ulong remaining = quotient; - while ((remaining % factor) == 0UL) - { - remaining /= factor; - exponent++; - } + bool factorIsPrime = GetOrComputePrimality(factor); + if (factorIsPrime) + { + int exponent = 1; + ulong remaining = quotient; + while ((remaining % factor) == 0UL) + { + remaining /= factor; + exponent++; + } - AddFactor(counts, factor, exponent); + AddFactor(counts, factor, exponent); - if (remaining > 1UL) - { - FactorCompletely(remaining, counts, knownComposite: false); - } + if (remaining > 1UL) + { + FactorCompletely(remaining, counts, knownComposite: false); + } - return; - } + return; + } - FactorCompletely(factor, counts, knownComposite: true); + FactorCompletely(factor, counts, knownComposite: true); - bool quotientIsPrime; - if (quotient == factor) - { - quotientIsPrime = false; - } - else - { - quotientIsPrime = GetOrComputePrimality(quotient); - } + bool quotientIsPrime; + if (quotient == factor) + { + quotientIsPrime = false; + } + else + { + quotientIsPrime = GetOrComputePrimality(quotient); + } - if (quotientIsPrime) - { - AddFactor(counts, quotient, 1); - } - else - { - FactorCompletely(quotient, counts, knownComposite: true); - } - } + if (quotientIsPrime) + { + AddFactor(counts, quotient, 1); + } + else + { + FactorCompletely(quotient, counts, knownComposite: true); + } + } - private static ulong PollardRhoStrict(ulong n) + private static ulong PollardRhoStrict(ulong n) { if ((n & 1UL) == 0UL) { diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index bd8a473..fd5bba1 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -2,79 +2,81 @@ namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - internal enum PrimeOrderMode - { - Heuristic, - Strict, - } + internal enum PrimeOrderMode + { + Heuristic, + Strict, + } - internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) - { - public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24, maxPowChecks: 24, PrimeOrderMode.Heuristic); - public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); + internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) + { + public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 40960, maxPowChecks: 24, PrimeOrderMode.Heuristic); + public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); - public readonly uint SmallFactorLimit = smallFactorLimit; - public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; - public readonly int MaxPowChecks = maxPowChecks; - public readonly PrimeOrderMode Mode = mode; - } + public readonly uint SmallFactorLimit = smallFactorLimit; + public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; + public readonly int MaxPowChecks = maxPowChecks; + public readonly PrimeOrderMode Mode = mode; + } - internal enum PrimeOrderHeuristicDevice - { - Cpu, - Gpu, - } + internal enum PrimeOrderHeuristicDevice + { + Cpu, + Gpu, + } - public static ulong Calculate( - ulong prime, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderSearchConfig config, - PrimeOrderHeuristicDevice device) - { - var scope = UsePow2Mode(device); - ulong order = CalculateInternal(prime, previousOrder, divisorData, config); - scope.Dispose(); - return order; - } + // TODO: Remove branching on the CPU / GPU paths + // TODO: Split big non-static functions into smaller / extract static code to limit JIT time + public static ulong Calculate( + ulong prime, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderSearchConfig config, + PrimeOrderHeuristicDevice device) + { + var scope = UsePow2Mode(device); + ulong order = CalculateInternal(prime, previousOrder, divisorData, config); + scope.Dispose(); + return order; + } - public static UInt128 Calculate( - in UInt128 prime, - in UInt128? previousOrder, - in PrimeOrderSearchConfig config, - in PrimeOrderHeuristicDevice device) - { - var scope = UsePow2Mode(device); - MontgomeryDivisorData divisorData; - UInt128 result; - if (prime <= ulong.MaxValue) - { - ulong? previous = null; - if (previousOrder.HasValue) - { - UInt128 previousValue = previousOrder.Value; - if (previousValue <= ulong.MaxValue) - { - previous = (ulong)previousValue; - } - else - { - previous = ulong.MaxValue; - } - } + public static UInt128 Calculate( + in UInt128 prime, + in UInt128? previousOrder, + in PrimeOrderSearchConfig config, + in PrimeOrderHeuristicDevice device) + { + var scope = UsePow2Mode(device); + MontgomeryDivisorData divisorData; + UInt128 result; + if (prime <= ulong.MaxValue) + { + ulong? previous = null; + if (previousOrder.HasValue) + { + UInt128 previousValue = previousOrder.Value; + if (previousValue <= ulong.MaxValue) + { + previous = (ulong)previousValue; + } + else + { + previous = ulong.MaxValue; + } + } - ulong prime64 = (ulong)prime; - divisorData = MontgomeryDivisorData.FromModulus(prime64); - ulong order64 = Calculate(prime64, previous, divisorData, config, device); - result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; - } - else - { - divisorData = default; - result = CalculateWideInternal(prime, previousOrder, divisorData, config); - } + ulong prime64 = (ulong)prime; + divisorData = MontgomeryDivisorData.FromModulus(prime64); + ulong order64 = Calculate(prime64, previous, divisorData, config, device); + result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; + } + else + { + divisorData = default; + result = CalculateWideInternal(prime, previousOrder, divisorData, config); + } - scope.Dispose(); - return result; - } + scope.Dispose(); + return result; + } } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 78ab92a..9351b65 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -10,1020 +10,1022 @@ namespace PerfectNumbers.Core; public static class ULongExtensions { - private const int Pow2WindowSize = 8; - private const ulong Pow2WindowFallbackThreshold = 32UL; - private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); + private const int Pow2WindowSize = 8; + private const ulong Pow2WindowFallbackThreshold = 32UL; + private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSize - 1); - public static ulong CalculateOrder(this ulong q) - { - if (q <= 2UL) - { - return 0UL; - } + public static ulong CalculateOrder(this ulong q) + { + if (q <= 2UL) + { + return 0UL; + } - ulong order = q - 1UL, prime, temp; - uint[] smallPrimes = PrimesGenerator.SmallPrimes; - ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; + ulong order = q - 1UL, prime, temp; + uint[] smallPrimes = PrimesGenerator.SmallPrimes; + ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; - int i = 0, primesLength = smallPrimes.Length; - UInt128 q128 = q, - cycle = MersenneDivisorCycles.GetCycle(q128); - // TODO: When the shared cycle snapshot cannot serve this divisor, trigger an on-demand - // GPU computation (respecting the configured device) without promoting the result into - // the cache so the order calculator still benefits from cycle stepping while keeping the - // single-block memory plan intact. + int i = 0, primesLength = smallPrimes.Length; + UInt128 q128 = q, + cycle = MersenneDivisorCycles.GetCycle(q128); + // TODO: When the shared cycle snapshot cannot serve this divisor, trigger an on-demand + // GPU computation (respecting the configured device) without promoting the result into + // the cache so the order calculator still benefits from cycle stepping while keeping the + // single-block memory plan intact. - for (; i < primesLength; i++) - { - if (smallPrimesPow2[i] > order) - { - break; - } + for (; i < primesLength; i++) + { + if (smallPrimesPow2[i] > order) + { + break; + } - prime = smallPrimes[i]; - // TODO: Replace this `%` driven factor peeling with the divisor-cycle aware - // factoring helper so large orders reuse the cached remainders highlighted in - // the latest divisor-cycle benchmarks instead of recomputing slow modulo checks. - while (order % prime == 0UL) - { - temp = order / prime; - if (temp.PowModWithCycle(q128, cycle) == UInt128.One) - { - order = temp; - } - else - { - break; - } - } - } + prime = smallPrimes[i]; + // TODO: Replace this `%` driven factor peeling with the divisor-cycle aware + // factoring helper so large orders reuse the cached remainders highlighted in + // the latest divisor-cycle benchmarks instead of recomputing slow modulo checks. + while (order % prime == 0UL) + { + temp = order / prime; + if (temp.PowModWithCycle(q128, cycle) == UInt128.One) + { + order = temp; + } + else + { + break; + } + } + } - return order; - } + return order; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static int GetBitLength(this ulong value) - { - return 64 - BitOperations.LeadingZeroCount(value); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static int GetBitLength(this ulong value) + { + return 64 - BitOperations.LeadingZeroCount(value); + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MultiplyShiftRight(ulong value, ulong multiplier, int shift) - { - UInt128 product = (UInt128)value * multiplier; - return (ulong)(product >> shift); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MultiplyShiftRight(ulong value, ulong multiplier, int shift) + { + UInt128 product = (UInt128)value * multiplier; + return (ulong)(product >> shift); + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift) - { - ulong high = value >> shift; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MultiplyShiftRightShiftFirst(ulong value, ulong multiplier, int shift) + { + ulong high = value >> shift; ulong mask = (1UL << shift) - 1UL; - ulong low = value & mask; + ulong low = value & mask; - UInt128 highContribution = (UInt128)high * multiplier; - UInt128 lowContribution = (UInt128)low * multiplier; + UInt128 highContribution = (UInt128)high * multiplier; + UInt128 lowContribution = (UInt128)low * multiplier; - UInt128 combined = highContribution + (lowContribution >> shift); - return (ulong)combined; - } + UInt128 combined = highContribution + (lowContribution >> shift); + return (ulong)combined; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong FastDiv64(this ulong value, ulong divisor, ulong mul) - { - ulong quotient = (ulong)(((UInt128)value * mul) >> 64); - UInt128 remainder = (UInt128)value - ((UInt128)quotient * divisor); - if (remainder >= divisor) - { - quotient++; - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong FastDiv64(this ulong value, ulong divisor, ulong mul) + { + ulong quotient = (ulong)(((UInt128)value * mul) >> 64); + UInt128 remainder = (UInt128)value - ((UInt128)quotient * divisor); + if (remainder >= divisor) + { + quotient++; + } - return quotient; - } + return quotient; + } - public const ulong WordBitMask = 0xFFFFUL; + public const ulong WordBitMask = 0xFFFFUL; - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static bool IsPrimeCandidate(this ulong n) - { - int i = 0; - uint[] smallPrimes = PrimesGenerator.SmallPrimes; - ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; - int len = smallPrimes.Length; - ulong p; - for (; i < len; i++) - { - if (smallPrimesPow2[i] > n) - { - break; - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool IsPrimeCandidate(this ulong n) + { + int i = 0; + uint[] smallPrimes = PrimesGenerator.SmallPrimes; + ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; + int len = smallPrimes.Length; + ulong p; + for (; i < len; i++) + { + if (smallPrimesPow2[i] > n) + { + break; + } - p = smallPrimes[i]; - // TODO: Swap this modulo check for the shared small-prime cycle filter once the - // divisor-cycle cache is mandatory, matching the PrimeTester improvements noted in - // the CPU sieve benchmarks. - if ((n % p) == 0UL) - { - return n == p; - } - } + p = smallPrimes[i]; + // TODO: Swap this modulo check for the shared small-prime cycle filter once the + // divisor-cycle cache is mandatory, matching the PrimeTester improvements noted in + // the CPU sieve benchmarks. + if ((n % p) == 0UL) + { + return n == p; + } + } - return true; - } + return true; + } - // Benchmarks (Mod5ULongBenchmarks) show the direct `% 5` is still cheaper (~0.26 ns vs 0.43 ns), so keep the modulo until a faster lookup is proven. - // (Mod8/Mod10 stay masked because they win; Mod5 currently does not.) - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod10(this ulong value) => (value & 1UL) == 0UL - ? (value % 5UL) switch - { - 0UL => 0UL, - 1UL => 6UL, - 2UL => 2UL, - 3UL => 8UL, - _ => 4UL, - } - : (value % 5UL) switch - { - 0UL => 5UL, - 1UL => 1UL, - 2UL => 7UL, - 3UL => 3UL, - _ => 9UL, - }; + // Benchmarks (Mod5ULongBenchmarks) show the direct `% 5` is still cheaper (~0.26 ns vs 0.43 ns), so keep the modulo until a faster lookup is proven. + // (Mod8/Mod10 stay masked because they win; Mod5 currently does not.) + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Mod10(this ulong value) => (value & 1UL) == 0UL + ? (value % 5UL) switch + { + 0UL => 0UL, + 1UL => 6UL, + 2UL => 2UL, + 3UL => 8UL, + _ => 4UL, + } + : (value % 5UL) switch + { + 0UL => 5UL, + 1UL => 1UL, + 2UL => 7UL, + 3UL => 3UL, + _ => 9UL, + }; - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod128(this ulong value) => value & 127UL; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Mod128(this ulong value) => value & 127UL; - // Benchmarks confirm `%` beats our current Mod5/Mod3 helpers for 64-bit inputs, so leave these modulo operations in place until a superior lookup is available. - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static void Mod10_8_5_3(this ulong value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3) - { - mod8 = value & 7UL; - // Benchmarks show `%` remains faster for the Mod5/Mod3 pair on ulong, so we keep the modulo path here for now. - mod5 = value % 5UL; - mod3 = value % 3UL; + // Benchmarks confirm `%` beats our current Mod5/Mod3 helpers for 64-bit inputs, so leave these modulo operations in place until a superior lookup is available. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void Mod10_8_5_3(this ulong value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3) + { + mod8 = value & 7UL; + // Benchmarks show `%` remains faster for the Mod5/Mod3 pair on ulong, so we keep the modulo path here for now. + mod5 = value % 5UL; + mod3 = value % 3UL; - mod10 = (mod8 & 1UL) == 0UL - ? mod5 switch - { - 0UL => 0UL, - 1UL => 6UL, - 2UL => 2UL, - 3UL => 8UL, - _ => 4UL, - } - : mod5 switch - { - 0UL => 5UL, - 1UL => 1UL, - 2UL => 7UL, - 3UL => 3UL, - _ => 9UL, - }; - } + mod10 = (mod8 & 1UL) == 0UL + ? mod5 switch + { + 0UL => 0UL, + 1UL => 6UL, + 2UL => 2UL, + 3UL => 8UL, + _ => 4UL, + } + : mod5 switch + { + 0UL => 5UL, + 1UL => 1UL, + 2UL => 7UL, + 3UL => 3UL, + _ => 9UL, + }; + } - // Mod5/Mod3 lookup tables are currently slower on 64-bit operands; keep the direct modulo until benchmarks flip. - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static void Mod10_8_5_3Steps(this ulong value, out ulong step10, out ulong step8, out ulong step5, out ulong step3) - { - ulong mod8 = value & 7UL; - // Same rationale: `%` wins in Mod5/Mod3 benches today, so avoid swapping until a faster lookup exists. - ulong mod5 = value % 5UL; - ulong mod3 = value % 3UL; - ulong mod10 = (mod8 & 1UL) == 0UL - ? mod5 switch - { - 0UL => 0UL, - 1UL => 6UL, - 2UL => 2UL, - 3UL => 8UL, - _ => 4UL, - } - : mod5 switch - { - 0UL => 5UL, - 1UL => 1UL, - 2UL => 7UL, - 3UL => 3UL, - _ => 9UL, - }; + // Mod5/Mod3 lookup tables are currently slower on 64-bit operands; keep the direct modulo until benchmarks flip. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void Mod10_8_5_3Steps(this ulong value, out ulong step10, out ulong step8, out ulong step5, out ulong step3) + { + ulong mod8 = value & 7UL; + // Same rationale: `%` wins in Mod5/Mod3 benches today, so avoid swapping until a faster lookup exists. + ulong mod5 = value % 5UL; + ulong mod3 = value % 3UL; + ulong mod10 = (mod8 & 1UL) == 0UL + ? mod5 switch + { + 0UL => 0UL, + 1UL => 6UL, + 2UL => 2UL, + 3UL => 8UL, + _ => 4UL, + } + : mod5 switch + { + 0UL => 5UL, + 1UL => 1UL, + 2UL => 7UL, + 3UL => 3UL, + _ => 9UL, + }; - step10 = mod10 + mod10; - if (step10 >= 10UL) - { - step10 -= 10UL; - } + step10 = mod10 + mod10; + if (step10 >= 10UL) + { + step10 -= 10UL; + } - step8 = (mod8 + mod8) & 7UL; + step8 = (mod8 + mod8) & 7UL; - step5 = mod5 + mod5; - if (step5 >= 5UL) - { - step5 -= 5UL; - } + step5 = mod5 + mod5; + if (step5 >= 5UL) + { + step5 -= 5UL; + } - step3 = mod3 + mod3; - if (step3 >= 3UL) - { - step3 -= 3UL; - } - } + step3 = mod3 + mod3; + if (step3 >= 3UL) + { + step3 -= 3UL; + } + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b); - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulHigh(this ulong x, ulong y) - { - // TODO: Investigate replacing this manual decomposition with the UInt128-based implementation - // for CPU callers; the latest benchmarks show the intrinsic path is an order of magnitude - // faster, while GPU code can keep using GpuUInt128.MulHigh. - ulong xLow = (uint)x; - ulong xHigh = x >> 32; - ulong yLow = (uint)y; - ulong yHigh = y >> 32; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulHigh(this ulong x, ulong y) + { + // TODO: Investigate replacing this manual decomposition with the UInt128-based implementation + // for CPU callers; the latest benchmarks show the intrinsic path is an order of magnitude + // faster, while GPU code can keep using GpuUInt128.MulHigh. + ulong xLow = (uint)x; + ulong xHigh = x >> 32; + ulong yLow = (uint)y; + ulong yHigh = y >> 32; - ulong w1 = xLow * yHigh; - ulong w2 = xHigh * yLow; - ulong w3 = xLow * yLow; + ulong w1 = xLow * yHigh; + ulong w2 = xHigh * yLow; + ulong w3 = xLow * yLow; - // Matching the layout used in GpuUInt128.MulHigh: introducing the - // intermediate result looks like one extra store, but it lets RyuJIT keep - // the accumulated high word entirely in registers. Without this explicit - // local the JIT spills the partial sum, which is where the performance - // regression in the benchmarks came from. - ulong result = (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32); - result += ((w3 >> 32) + (uint)w1 + (uint)w2) >> 32; - return result; - } + // Matching the layout used in GpuUInt128.MulHigh: introducing the + // intermediate result looks like one extra store, but it lets RyuJIT keep + // the accumulated high word entirely in registers. Without this explicit + // local the JIT spills the partial sum, which is where the performance + // regression in the benchmarks came from. + ulong result = (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32); + result += ((w3 >> 32) + (uint)w1 + (uint)w2) >> 32; + return result; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulHighGpuCompatible(this ulong x, ulong y) - { - GpuUInt128 product = new(x); - product.Mul64(new GpuUInt128(y)); - return product.High; - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulHighGpuCompatible(this ulong x, ulong y) + { + GpuUInt128 product = new(x); + product.Mul64(new GpuUInt128(y)); + return product.High; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus) - { - 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; + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus) + { + 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; - while (exponent != 0UL) - { - if ((exponent & 1UL) != 0UL) - { - result = MulMod64(result, value, modulus); - } + while (exponent != 0UL) + { + if ((exponent & 1UL) != 0UL) + { + result = MulMod64(result, value, modulus); + } - value = MulMod64(value, value, modulus); - exponent >>= 1; - } + value = MulMod64(value, value, modulus); + exponent >>= 1; + } - return result; - } + 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)] + // 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)] - public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) - { - // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64, - // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs. - GpuUInt128 state = new(a % modulus); - return state.MulMod(b, modulus); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) + { + // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64, + // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs. + GpuUInt128 state = new(a % modulus); + return state.MulMod(b, modulus); + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) - { - // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the - // 5-40× slowdown seen across real-world operand distributions. - GpuUInt128 state = new(a); - return state.MulModWithNativeModulo(b, modulus); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) + { + // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the + // 5-40× slowdown seen across real-world operand distributions. + GpuUInt128 state = new(a); + return state.MulModWithNativeModulo(b, modulus); + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime) - { - ulong tLow = unchecked(a * b); - ulong m = unchecked(tLow * nPrime); - ulong mTimesModulusLow = unchecked(m * modulus); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime) + { + ulong tLow = unchecked(a * b); + ulong m = unchecked(tLow * nPrime); + ulong mTimesModulusLow = unchecked(m * modulus); - ulong result = unchecked(a.MulHigh(b) + m.MulHigh(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL)); - if (result >= modulus) - { - result -= modulus; - } + ulong result = unchecked(a.MulHigh(b) + m.MulHigh(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL)); + if (result >= modulus) + { + result -= modulus; + } - return result; - } + return result; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) - { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWindowedCpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { ulong modulus = divisor.Modulus; // This should never happen in production code. If it does, we need to fix it. - // if (exponent == 0UL) - // { - // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; - // } - - if (exponent <= Pow2WindowFallbackThreshold) - { - return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); - } - - int bitLength = GetPortableBitLength(exponent); - int windowSize = GetWindowSize(bitLength); - int oddPowerCount = 1 << (windowSize - 1); - - ulong result = divisor.MontgomeryOne; - ulong nPrime = divisor.NPrime; - - ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); - Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); - InitializeMontgomeryOddPowersCpu(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; - } - - ThreadStaticPools.UlongPool.Return(oddPowersArray); - - if (!keepMontgomery) - { - result = result.MontgomeryMultiply(1UL, modulus, nPrime); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) - { - ulong modulus = divisor.Modulus; - if (exponent == 0UL) - { - return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; - } - - if (exponent <= Pow2WindowFallbackThreshold) - { - return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); - } - - if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult)) - { - return gpuResult; - } - - return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong Pow2MontgomeryModWindowedKernel(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) - { - ulong modulus = divisor.Modulus; - if (exponent == 0UL) - { - return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; - } - - if (exponent <= Pow2WindowFallbackThreshold) - { - return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); - } - - int bitLength = GetPortableBitLength(exponent); - int windowSize = GetWindowSize(bitLength); - ulong result = divisor.MontgomeryOne; - ulong nPrime = divisor.NPrime; - - 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; - ulong multiplier = ComputeMontgomeryOddPower(windowValue, divisor, modulus, nPrime); - result = result.MontgomeryMultiply(multiplier, modulus, nPrime); - - index = windowStart - 1; - } - - if (!keepMontgomery) - { - result = result.MontgomeryMultiply(1UL, modulus, nPrime); - } - - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeMontgomeryOddPower(ulong exponent, in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime) - { - ulong baseValue = divisor.MontgomeryTwo; - ulong power = divisor.MontgomeryOne; - ulong remaining = exponent; - - while (remaining != 0UL) - { - if ((remaining & 1UL) != 0UL) - { - power = power.MontgomeryMultiply(baseValue, modulus, nPrime); - } - - remaining >>= 1; - if (remaining == 0UL) - { - break; - } - - baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime); - } - - return power; - } - - private static class Pow2MontgomeryGpuExecutor - { - private static readonly ConcurrentDictionary, MontgomeryDivisorData, byte, ArrayView1D>> KernelCache = new(); - - public static bool TryExecute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery, out ulong result) - { - result = 0UL; - - GpuKernelLease lease; - try - { - lease = GpuKernelPool.GetKernel(useGpuOrder: true); - } - catch (Exception) - { - return false; - } - - var execution = lease.EnterExecutionScope(); - try - { - Accelerator accelerator = lease.Accelerator; - if (accelerator.AcceleratorType == AcceleratorType.CPU) - { - return false; - } - - var kernel = KernelCache.GetOrAdd(accelerator, static accel => - { - var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, byte, ArrayView1D>(Pow2MontgomeryKernel); - var launcher = KernelUtil.GetKernel(loaded); - return launcher.CreateLauncherDelegate, MontgomeryDivisorData, byte, ArrayView1D>>(); - }); - - using var exponentBuffer = accelerator.Allocate1D(1); - using var resultBuffer = accelerator.Allocate1D(1); - - exponentBuffer.View.CopyFromCPU(ref exponent, 1); - resultBuffer.MemSetToZero(); - - byte keepFlag = keepMontgomery ? (byte)1 : (byte)0; - AcceleratorStream stream = lease.Stream; - kernel(stream, 1, exponentBuffer.View, divisor, keepFlag, resultBuffer.View); - stream.Synchronize(); - - resultBuffer.View.CopyToCPU(ref result, 1); - return true; - } - catch (AcceleratorException) - { - return false; - } - catch (NotSupportedException) - { - return false; - } - finally - { - execution.Dispose(); - lease.Dispose(); - } - } - - private static void Pow2MontgomeryKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, byte keepMontgomery, ArrayView1D results) - { - results[index] = Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0); - } - } - - [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 remainingExponent = exponent; - - while (remainingExponent != 0UL) - { - if ((remainingExponent & 1UL) != 0UL) - { - result = result.MontgomeryMultiply(baseVal, modulus, nPrime); - } - - remainingExponent >>= 1; - if (remainingExponent == 0UL) - { - break; - } - - baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime); - } - - return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static int GetPortableBitLength(ulong value) - { - if (value == 0UL) - { - return 0; - } - - return 64 - XMath.LeadingZeroCount(value); - } - - private static int GetWindowSize(int bitLength) - { - if (bitLength <= 6) - { - return Math.Max(bitLength, 1); - } - - if (bitLength <= 23) - { - return 4; - } - - if (bitLength <= 79) - { - return 5; - } - - if (bitLength <= 239) - { - return 6; - } - - if (bitLength <= 671) - { - return 7; - } - - return Pow2WindowSize; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span oddPowers) - { - oddPowers[0] = divisor.MontgomeryTwo; - if (oddPowers.Length == 1) - { - return; - } - - ulong previous; - ulong square = divisor.MontgomeryTwoSquared; - - for (int i = 1; i < oddPowers.Length; i++) - { - previous = oddPowers[i - 1]; - oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) - { - ulong rotationCount = exponent % cycleLength; - return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) - { - return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) - { - ulong rotationCount = exponent % cycleLength; - return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) - { - return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static UInt128 PowMod(this ulong exponent, UInt128 modulus) - { - UInt128 result = UInt128.One; - ulong exponentLoopIndex = 0UL; - - // 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; - - // Any number mod 1 is 0 - if (modulus == UInt128.One) - return UInt128.Zero; - - // For small exponents, do classic method - if (exponent < 64 || modulus < 4) - { - for (; exponentLoopIndex < exponent; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } - - return result; - } - - // 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 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; - - // 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; - - // Any number mod 1 is 0 - if (modulus == UInt128.One) - return UInt128.Zero; - - // For small exponents, do classic method - if (exponent < 64 || modulus < 4) - { - for (; exponentLoopIndex < exponent; exponentLoopIndex++) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - } - - return result; - } - - // 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; - - // 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; - - // Any number mod 1 is 0 - if (modulus == UInt128.One) - return UInt128.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; - } - - return result; - } - - // 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); - } - - // 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; - - // 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; - - // Any number mod 1 is 0 - 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; - } - - return result; - } - - // 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); - } - - // 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; - - // We're reusing "zero" as rotation index for just a little better performance - while (zero < rotationCount) - { - result <<= 1; - if (result >= modulus) - result -= modulus; - - zero += one; - } - - return result; - } - - /// - /// Computes 2^exponent mod modulus using iterative CRT composition from mod 10 up to modulus. - /// Only for modulus >= 10 and reasonable size. - /// - public static UInt128 PowModCrt(this ulong exponent, UInt128 modulus, MersenneDivisorCycles cycles) - { - if (modulus < 10) - return PowMod(exponent, modulus); // fallback to classic - - // Use cycle length 4 for mod 10 - UInt128 currentModulus = 10, - cycle, - modulusCandidate = 11, - remainderForCandidate, - result = PowModWithCycle(exponent, 10, 4), - zero = UInt128.Zero; - - for (; modulusCandidate <= modulus; modulusCandidate++) - { - cycle = MersenneDivisorCycles.GetCycle(modulusCandidate); - remainderForCandidate = cycle > zero - ? PowModWithCycle(exponent, modulusCandidate, cycle) - : PowMod(exponent, modulusCandidate); - - // Solve x ≡ result mod currentModulus - // x ≡ remM mod m - // Find x mod (currentModulus * m) - // Since currentModulus and m are coprime, use CRT: - // x = result + currentModulus * t, where t ≡ (remM - result) * inv(currentModulus, m) mod m - - // TODO: Replace this `% modulusCandidate` with the cached residue helper derived from Mod10_8_5_3Benchmarks so CRT - // composition avoids repeated modulo divisions when combining residues for large divisor sets. - result += currentModulus * ((remainderForCandidate + modulusCandidate - (result % modulusCandidate)) * ModInverse(currentModulus, modulusCandidate) % modulusCandidate); - currentModulus *= modulusCandidate; - - if (currentModulus >= modulus) - break; - } - - // TODO: Swap this final `% modulus` with the pooled remainder cache so the CRT result write-back avoids one more division, - // aligning with the optimizations captured in Mod10_8_5_3Benchmarks. - return result % modulus; - } - - // Helper: modular inverse (extended Euclidean algorithm) - private static UInt128 ModInverse(UInt128 a, UInt128 m) - { - UInt128 m0 = m, - originalA, - originalM, - temp, - x0 = 0, - x1 = 1; - - if (m == 1) - { - return 0; - } - - while (a > 1) - { - originalA = a; - originalM = m; - m = originalA % originalM; - a = originalM; - temp = x0; - x0 = x1 - (originalA / originalM) * x0; - x1 = temp; - } - - if (x1 < 0) - { - x1 += m0; - } - - return x1; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static bool SharesFactorWithExponentMinusOne(this ulong exponent) - { - ulong prime, value = exponent - 1UL; - value >>= BitOperations.TrailingZeroCount(value); - uint[] smallPrimes = PrimesGenerator.SmallPrimes; - ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; - int i = 1, smallPrimesLength = smallPrimes.Length; - - for (; i < smallPrimesLength && smallPrimesPow2[i] <= value; i++) - { - prime = smallPrimes[i]; - if (value % prime != 0UL) - { - continue; - } - - if (exponent % prime.CalculateOrder() == 0UL) - { - return true; - } - - do - { - value /= prime; - } - while (value % prime == 0UL); - } - - if (value > 1UL && exponent % value.CalculateOrder() == 0UL) - { - return true; - } - - return false; - } + // if (exponent == 0UL) + // { + // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; + // } + + // There is no point in having extra branch, because in real use-cases it'll barely ever hit. + // Out of all 3_358_944_545 Pow2MontgomeryModWindowedCpu calls, we hit Pow2MontgomeryModSingleBit only 772 times. + // if (exponent <= Pow2WindowFallbackThreshold) + // { + // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + // } + + int bitLength = GetPortableBitLength(exponent); + int windowSize = GetWindowSize(bitLength); + int oddPowerCount = 1 << (windowSize - 1); + + ulong result = divisor.MontgomeryOne; + ulong nPrime = divisor.NPrime; + + ulong[] oddPowersArray = ThreadStaticPools.UlongPool.Rent(oddPowerCount); + Span oddPowers = oddPowersArray.AsSpan(0, oddPowerCount); + InitializeMontgomeryOddPowersCpu(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; + } + + ThreadStaticPools.UlongPool.Return(oddPowersArray); + + if (!keepMontgomery) + { + result = result.MontgomeryMultiply(1UL, modulus, nPrime); + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWindowedGpu(this ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { + // We shouldn't hit it in production code. If we do, we're doing something wrong. + // ulong modulus = divisor.Modulus; + // if (exponent == 0UL) + // { + // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; + // } + + // We barely ever hit it in production code. It's e.g. 772 calls out of billions + // if (exponent <= Pow2WindowFallbackThreshold) + // { + // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + // } + + if (Pow2MontgomeryGpuExecutor.TryExecute(exponent, divisor, keepMontgomery, out ulong gpuResult)) + { + return gpuResult; + } + + return Pow2MontgomeryModWindowedCpu(exponent, divisor, keepMontgomery); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong Pow2MontgomeryModWindowedKernel(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery) + { + ulong modulus = divisor.Modulus; + // This should never happen in production code + // if (exponent == 0UL) + // { + // return keepMontgomery ? divisor.MontgomeryOne : 1UL % modulus; + // } + + // We barely every hit it from production code. It's only 772 calls out of billions. + // if (exponent <= Pow2WindowFallbackThreshold) + // { + // return Pow2MontgomeryModSingleBit(exponent, divisor, keepMontgomery); + // } + + int bitLength = GetPortableBitLength(exponent); + int windowSize = GetWindowSize(bitLength); + ulong result = divisor.MontgomeryOne; + ulong nPrime = divisor.NPrime; + + 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; + ulong multiplier = ComputeMontgomeryOddPower(windowValue, divisor, modulus, nPrime); + result = result.MontgomeryMultiply(multiplier, modulus, nPrime); + + index = windowStart - 1; + } + + if (!keepMontgomery) + { + result = result.MontgomeryMultiply(1UL, modulus, nPrime); + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeMontgomeryOddPower(ulong exponent, in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime) + { + ulong baseValue = divisor.MontgomeryTwo; + ulong power = divisor.MontgomeryOne; + ulong remaining = exponent; + + while (remaining != 0UL) + { + if ((remaining & 1UL) != 0UL) + { + power = power.MontgomeryMultiply(baseValue, modulus, nPrime); + } + + remaining >>= 1; + if (remaining == 0UL) + { + break; + } + + baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime); + } + + return power; + } + + private static class Pow2MontgomeryGpuExecutor + { + private static readonly ConcurrentDictionary, MontgomeryDivisorData, byte, ArrayView1D>> KernelCache = new(); + + public static bool TryExecute(ulong exponent, in MontgomeryDivisorData divisor, bool keepMontgomery, out ulong result) + { + result = 0UL; + + GpuKernelLease lease = GpuKernelPool.GetKernel(useGpuOrder: true); + var execution = lease.EnterExecutionScope(); + try + { + Accelerator accelerator = lease.Accelerator; + // Keep this commented out. It should never happen in production code. + // if (accelerator.AcceleratorType == AcceleratorType.CPU) + // { + // return false; + // } + + var kernel = KernelCache.GetOrAdd(accelerator, static accel => + { + var loaded = accel.LoadAutoGroupedStreamKernel, MontgomeryDivisorData, byte, ArrayView1D>(Pow2MontgomeryKernel); + var launcher = KernelUtil.GetKernel(loaded); + return launcher.CreateLauncherDelegate, MontgomeryDivisorData, byte, ArrayView1D>>(); + }); + + using var exponentBuffer = accelerator.Allocate1D(1); + using var resultBuffer = accelerator.Allocate1D(1); + + exponentBuffer.View.CopyFromCPU(ref exponent, 1); + // We don't need to worry about any left-overs here. + // resultBuffer.MemSetToZero(); + + // TODO: Have 2 kernels for both cases to have branch-less solution. + byte keepFlag = keepMontgomery ? (byte)1 : (byte)0; + AcceleratorStream stream = lease.Stream; + kernel(stream, 1, exponentBuffer.View, divisor, keepFlag, resultBuffer.View); + stream.Synchronize(); + + resultBuffer.View.CopyToCPU(ref result, 1); + return true; + } + // Keep this commented. We don't want to catch any exceptions. All should crash the scanner. + // catch (AcceleratorException) + // { + // return false; + // } + // catch (NotSupportedException) + // { + // return false; + // } + finally + { + execution.Dispose(); + lease.Dispose(); + } + } + + private static void Pow2MontgomeryKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, byte keepMontgomery, ArrayView1D results) + { + results[index] = Pow2MontgomeryModWindowedKernel(exponents[index], divisor, keepMontgomery != 0); + } + } + + [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 remainingExponent = exponent; + + while (remainingExponent != 0UL) + { + if ((remainingExponent & 1UL) != 0UL) + { + result = result.MontgomeryMultiply(baseVal, modulus, nPrime); + } + + remainingExponent >>= 1; + if (remainingExponent == 0UL) + { + break; + } + + baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime); + } + + return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int GetPortableBitLength(ulong value) + { + // Keep this commented out. It will never happen in production code. + // if (value == 0UL) + // { + // return 0; + // } + + return 64 - XMath.LeadingZeroCount(value); + } + + private static int GetWindowSize(int bitLength) + { + if (bitLength <= 6) + { + return Math.Max(bitLength, 1); + } + + if (bitLength <= 23) + { + return 4; + } + + if (bitLength <= 79) + { + return 5; + } + + if (bitLength <= 239) + { + return 6; + } + + if (bitLength <= 671) + { + return 7; + } + + return Pow2WindowSize; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void InitializeMontgomeryOddPowersCpu(in MontgomeryDivisorData divisor, ulong modulus, ulong nPrime, Span oddPowers) + { + oddPowers[0] = divisor.MontgomeryTwo; + if (oddPowers.Length == 1) + { + return; + } + + ulong previous; + ulong square = divisor.MontgomeryTwoSquared; + + for (int i = 1; i < oddPowers.Length; i++) + { + previous = oddPowers[i - 1]; + oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWithCycleCpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) + { + ulong rotationCount = exponent % cycleLength; + return Pow2MontgomeryModWindowedCpu(rotationCount, divisor, keepMontgomery: false); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModFromCycleRemainderCpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) + { + return Pow2MontgomeryModWindowedCpu(reducedExponent, divisor, keepMontgomery: false); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWithCycleGpu(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) + { + ulong rotationCount = exponent % cycleLength; + return Pow2MontgomeryModWindowedGpu(rotationCount, divisor, keepMontgomery: false); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModFromCycleRemainderGpu(this ulong reducedExponent, in MontgomeryDivisorData divisor) + { + return Pow2MontgomeryModWindowedGpu(reducedExponent, divisor, keepMontgomery: false); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static UInt128 PowMod(this ulong exponent, UInt128 modulus) + { + UInt128 result = UInt128.One; + ulong exponentLoopIndex = 0UL; + + // 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; + + // Any number mod 1 is 0 + if (modulus == UInt128.One) + return UInt128.Zero; + + // For small exponents, do classic method + if (exponent < 64 || modulus < 4) + { + for (; exponentLoopIndex < exponent; exponentLoopIndex++) + { + result <<= 1; + if (result >= modulus) + result -= modulus; + } + + return result; + } + + // 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 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; + + // 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; + + // Any number mod 1 is 0 + if (modulus == UInt128.One) + return UInt128.Zero; + + // For small exponents, do classic method + if (exponent < 64 || modulus < 4) + { + for (; exponentLoopIndex < exponent; exponentLoopIndex++) + { + result <<= 1; + if (result >= modulus) + result -= modulus; + } + + return result; + } + + // 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; + + // 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; + + // Any number mod 1 is 0 + if (modulus == UInt128.One) + return UInt128.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; + } + + return result; + } + + // 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); + } + + // 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; + + // 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; + + // Any number mod 1 is 0 + 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; + } + + return result; + } + + // 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); + } + + // 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; + + // We're reusing "zero" as rotation index for just a little better performance + while (zero < rotationCount) + { + result <<= 1; + if (result >= modulus) + result -= modulus; + + zero += one; + } + + return result; + } + + /// + /// Computes 2^exponent mod modulus using iterative CRT composition from mod 10 up to modulus. + /// Only for modulus >= 10 and reasonable size. + /// + public static UInt128 PowModCrt(this ulong exponent, UInt128 modulus, MersenneDivisorCycles cycles) + { + if (modulus < 10) + return PowMod(exponent, modulus); // fallback to classic + + // Use cycle length 4 for mod 10 + UInt128 currentModulus = 10, + cycle, + modulusCandidate = 11, + remainderForCandidate, + result = PowModWithCycle(exponent, 10, 4), + zero = UInt128.Zero; + + for (; modulusCandidate <= modulus; modulusCandidate++) + { + cycle = MersenneDivisorCycles.GetCycle(modulusCandidate); + remainderForCandidate = cycle > zero + ? PowModWithCycle(exponent, modulusCandidate, cycle) + : PowMod(exponent, modulusCandidate); + + // Solve x ≡ result mod currentModulus + // x ≡ remM mod m + // Find x mod (currentModulus * m) + // Since currentModulus and m are coprime, use CRT: + // x = result + currentModulus * t, where t ≡ (remM - result) * inv(currentModulus, m) mod m + + // TODO: Replace this `% modulusCandidate` with the cached residue helper derived from Mod10_8_5_3Benchmarks so CRT + // composition avoids repeated modulo divisions when combining residues for large divisor sets. + result += currentModulus * ((remainderForCandidate + modulusCandidate - (result % modulusCandidate)) * ModInverse(currentModulus, modulusCandidate) % modulusCandidate); + currentModulus *= modulusCandidate; + + if (currentModulus >= modulus) + break; + } + + // TODO: Swap this final `% modulus` with the pooled remainder cache so the CRT result write-back avoids one more division, + // aligning with the optimizations captured in Mod10_8_5_3Benchmarks. + return result % modulus; + } + + // Helper: modular inverse (extended Euclidean algorithm) + private static UInt128 ModInverse(UInt128 a, UInt128 m) + { + UInt128 m0 = m, + originalA, + originalM, + temp, + x0 = 0, + x1 = 1; + + if (m == 1) + { + return 0; + } + + while (a > 1) + { + originalA = a; + originalM = m; + m = originalA % originalM; + a = originalM; + temp = x0; + x0 = x1 - (originalA / originalM) * x0; + x1 = temp; + } + + if (x1 < 0) + { + x1 += m0; + } + + return x1; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool SharesFactorWithExponentMinusOne(this ulong exponent) + { + ulong prime, value = exponent - 1UL; + value >>= BitOperations.TrailingZeroCount(value); + uint[] smallPrimes = PrimesGenerator.SmallPrimes; + ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; + int i = 1, smallPrimesLength = smallPrimes.Length; + + for (; i < smallPrimesLength && smallPrimesPow2[i] <= value; i++) + { + prime = smallPrimes[i]; + if (value % prime != 0UL) + { + continue; + } + + if (exponent % prime.CalculateOrder() == 0UL) + { + return true; + } + + do + { + value /= prime; + } + while (value % prime == 0UL); + } + + if (value > 1UL && exponent % value.CalculateOrder() == 0UL) + { + return true; + } + + return false; + } } diff --git a/PerfectNumbers.Core/UnboundedTaskScheduler.cs b/PerfectNumbers.Core/UnboundedTaskScheduler.cs index 82a8b1b..56e6f49 100644 --- a/PerfectNumbers.Core/UnboundedTaskScheduler.cs +++ b/PerfectNumbers.Core/UnboundedTaskScheduler.cs @@ -1,6 +1,3 @@ -using System.Threading; -using System.Threading.Tasks; - namespace PerfectNumbers.Core; public sealed class UnboundedTaskScheduler : TaskScheduler @@ -10,27 +7,21 @@ public sealed class UnboundedTaskScheduler : TaskScheduler private UnboundedTaskScheduler() { ThreadPool.SetMaxThreads(100_000, 100_000); - ThreadPool.SetMinThreads(10_300, 100); + ThreadPool.SetMinThreads(10_340, 100); } public override int MaximumConcurrencyLevel => int.MaxValue; protected override void QueueTask(Task task) { - ThreadPool.UnsafeQueueUserWorkItem(static state => + ThreadPool.UnsafeQueueUserWorkItem( (state) => { - (UnboundedTaskScheduler Scheduler, Task Task) tuple = ((UnboundedTaskScheduler, Task))state!; - tuple.Scheduler.TryExecuteTask(tuple.Task); + var (scheduler, task) = ((UnboundedTaskScheduler, Task))state!; + scheduler.TryExecuteTask(task); }, (this, task)); } - protected override bool TryExecuteTaskInline(Task task, bool taskWasPreviouslyQueued) - { - return TryExecuteTask(task); - } + protected override bool TryExecuteTaskInline(Task task, bool taskWasPreviouslyQueued) => TryExecuteTask(task); - protected override IEnumerable? GetScheduledTasks() - { - return null; - } + protected override IEnumerable? GetScheduledTasks() => null; }