diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs index 26b8e08..52647f1 100644 --- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs @@ -32,7 +32,7 @@ public class GpuUInt128Montgomery64Benchmarks [Benchmark(Baseline = true)] public ulong ExtensionMontgomeryMultiply() { - return Input.Left.MontgomeryMultiply(Input.Right, Input.Modulus, Input.NPrime); + return Input.Left.MontgomeryMultiplyCpu(Input.Right, Input.Modulus, Input.NPrime); } /// diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs index 5637499..5bd59d7 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs @@ -640,7 +640,7 @@ public class Pow2ModBenchmarks ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus); ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus); - ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); + ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime); return new MontgomeryDivisorData( modulus, @@ -736,22 +736,22 @@ public class Pow2ModBenchmarks if (exponent == 0UL) { - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent); for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(result, modulus, nPrime); if (((exponent >> bitIndex) & 1UL) != 0UL) { - result = result.MontgomeryMultiply(baseVal, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime); } } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor) @@ -772,35 +772,35 @@ public class Pow2ModBenchmarks ulong nibble = remainingExponent & 0xFUL; ulong base0 = baseVal; - ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime); - ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime); - ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime); + ulong base1 = base0.MontgomeryMultiplyCpu(base0, modulus, nPrime); + ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime); + ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime); if ((nibble & 1UL) != 0UL) { - result = result.MontgomeryMultiply(base0, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime); } if ((nibble & 2UL) != 0UL) { - result = result.MontgomeryMultiply(base1, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime); } if ((nibble & 4UL) != 0UL) { - result = result.MontgomeryMultiply(base2, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime); } if ((nibble & 8UL) != 0UL) { - result = result.MontgomeryMultiply(base3, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base3, modulus, nPrime); } - baseVal = base3.MontgomeryMultiply(base3, modulus, nPrime); + baseVal = base3.MontgomeryMultiplyCpu(base3, modulus, nPrime); remainingExponent >>= 4; } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MontgomeryDivisorData divisor) @@ -817,16 +817,16 @@ public class Pow2ModBenchmarks if (exponent == 0UL) { - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } const int WindowSize = 5; Span oddPowers = stackalloc ulong[1 << (WindowSize - 1)]; oddPowers[0] = baseVal; - ulong baseSquared = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime); + ulong baseSquared = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime); for (int i = 1; i < oddPowers.Length; i++) { - oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(baseSquared, modulus, nPrime); + oddPowers[i] = oddPowers[i - 1].MontgomeryMultiplyCpu(baseSquared, modulus, nPrime); } int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent); @@ -835,7 +835,7 @@ public class Pow2ModBenchmarks { if (((exponent >> bitIndex) & 1UL) == 0UL) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(result, modulus, nPrime); bitIndex--; continue; } @@ -852,14 +852,14 @@ public class Pow2ModBenchmarks int squares = bitIndex - windowStart + 1; for (int i = 0; i < squares; i++) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(result, modulus, nPrime); } - result = result.MontgomeryMultiply(oddPowers[(int)(windowValue >> 1)], modulus, nPrime); + result = result.MontgomeryMultiplyCpu(oddPowers[(int)(windowValue >> 1)], modulus, nPrime); bitIndex = windowStart - 1; } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } private static ulong Pow2ModBinary(ulong exponent, ulong modulus) diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs index 3c4f26e..37946e8 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs @@ -262,7 +262,7 @@ public class Pow2MontgomeryModBenchmarks ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus); ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, modulus); - ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); + ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime); return new MontgomeryDivisorData( modulus, @@ -300,22 +300,22 @@ public class Pow2MontgomeryModBenchmarks if (exponent == 0UL) { - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent); for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(result, modulus, nPrime); if (((exponent >> bitIndex) & 1UL) != 0UL) { - result = result.MontgomeryMultiply(baseVal, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime); } } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor) @@ -336,35 +336,35 @@ public class Pow2MontgomeryModBenchmarks ulong nibble = remainingExponent & 0xFUL; ulong base0 = baseVal; - ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime); - ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime); - ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime); + ulong base1 = base0.MontgomeryMultiplyCpu(base0, modulus, nPrime); + ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime); + ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime); if ((nibble & 1UL) != 0UL) { - result = result.MontgomeryMultiply(base0, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime); } if ((nibble & 2UL) != 0UL) { - result = result.MontgomeryMultiply(base1, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime); } if ((nibble & 4UL) != 0UL) { - result = result.MontgomeryMultiply(base2, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime); } if ((nibble & 8UL) != 0UL) { - result = result.MontgomeryMultiply(base3, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(base3, modulus, nPrime); } - baseVal = base3.MontgomeryMultiply(base3, modulus, nPrime); + baseVal = base3.MontgomeryMultiplyCpu(base3, modulus, nPrime); remainingExponent >>= 4; } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } // TODO: Update it with the latest implementations @@ -382,16 +382,16 @@ public class Pow2MontgomeryModBenchmarks if (exponent == 0UL) { - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } const int WindowSize = 5; Span oddPowers = stackalloc ulong[1 << (WindowSize - 1)]; oddPowers[0] = baseVal; - ulong baseSquared = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime); + ulong baseSquared = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime); for (int i = 1; i < oddPowers.Length; i++) { - oddPowers[i] = oddPowers[i - 1].MontgomeryMultiply(baseSquared, modulus, nPrime); + oddPowers[i] = oddPowers[i - 1].MontgomeryMultiplyCpu(baseSquared, modulus, nPrime); } int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent); @@ -400,7 +400,7 @@ public class Pow2MontgomeryModBenchmarks { if (((exponent >> bitIndex) & 1UL) == 0UL) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(result, modulus, nPrime); bitIndex--; continue; } @@ -417,14 +417,14 @@ public class Pow2MontgomeryModBenchmarks int squares = bitIndex - windowStart + 1; for (int i = 0; i < squares; i++) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(result, modulus, nPrime); } - result = result.MontgomeryMultiply(oddPowers[(int)(windowValue >> 1)], modulus, nPrime); + result = result.MontgomeryMultiplyCpu(oddPowers[(int)(windowValue >> 1)], modulus, nPrime); bitIndex = windowStart - 1; } - return result.MontgomeryMultiply(1UL, modulus, nPrime); + return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } // TODO: Update it with the latest implementations, if it exists in production code. Otherwise just remove TODO diff --git a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs index dc22238..afa681b 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs @@ -74,7 +74,7 @@ internal struct ExponentRemainderStepperCpu // single-block divisor-cycle snapshot so we can skip the powmod entirely and reuse the // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks. ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); - _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); + _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime); PreviousExponent = exponent; return ReduceCurrent(); } @@ -92,7 +92,7 @@ internal struct ExponentRemainderStepperCpu // TODO: Reuse the divisor-cycle derived Montgomery delta once the cache exposes single-cycle // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor. ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); - _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); + _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime); PreviousExponent = exponent; return _currentMontgomery == _montgomeryOne; } @@ -126,7 +126,7 @@ internal struct ExponentRemainderStepperCpu // TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead // of relying on the caller-provided delta so incremental scans remain in sync with the snapshot // without computing additional cycles or mutating cache state. - _currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime); + _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(montgomeryDelta, _modulus, _nPrime); PreviousExponent = exponent; isUnity = _currentMontgomery == _montgomeryOne; return true; @@ -138,5 +138,5 @@ internal struct ExponentRemainderStepperCpu PreviousExponent = exponent; } - private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime); + private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiplyCpu(1UL, _modulus, _nPrime); } diff --git a/PerfectNumbers.Core/MontgomeryDivisorData.cs b/PerfectNumbers.Core/MontgomeryDivisorData.cs index 9167c30..97ba2d8 100644 --- a/PerfectNumbers.Core/MontgomeryDivisorData.cs +++ b/PerfectNumbers.Core/MontgomeryDivisorData.cs @@ -34,7 +34,7 @@ public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus); ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, modulus); - ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiply(montgomeryTwo, montgomeryTwo, modulus, nPrime); + ulong montgomeryTwoSquared = ULongExtensions.MontgomeryMultiplyCpu(montgomeryTwo, montgomeryTwo, modulus, nPrime); return new MontgomeryDivisorData( modulus, diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index e070135..31ba9da 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -53,8 +53,7 @@ internal static partial class PrimeOrderCalculator ulong phi = prime - 1UL; - // if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) - if (PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) + if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) { return gpuOrder; } @@ -1070,7 +1069,7 @@ internal static partial class PrimeOrderCalculator ulong remaining = value; bool gpuFactored = false; - // if (IsGpuHeuristicDevice) + if (IsGpuHeuristicDevice) { gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); } @@ -1079,7 +1078,6 @@ internal static partial class PrimeOrderCalculator Stack? compositeStack = null; PartialFactorResult result; PrimeTester primeTester = PrimeTester.Exclusive; - CancellationToken ct = CancellationToken.None; if (!gpuFactored) { @@ -1152,7 +1150,7 @@ internal static partial class PrimeOrderCalculator // Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}"); // bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeCpu(composite, ct); + bool isPrime = PrimeTester.IsPrimeGpu(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) @@ -1220,7 +1218,7 @@ internal static partial class PrimeOrderCalculator // Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}"); // bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeCpu(composite, ct); + bool isPrime = PrimeTester.IsPrimeGpu(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); diff --git a/PerfectNumbers.Core/ULongExtensions.Cpu.cs b/PerfectNumbers.Core/ULongExtensions.Cpu.cs index 24e415c..5be09dc 100644 --- a/PerfectNumbers.Core/ULongExtensions.Cpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Cpu.cs @@ -40,7 +40,7 @@ public static partial class ULongExtensions { if (((exponent >> index) & 1UL) == 0UL) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(result, modulus, nPrime); index--; continue; } @@ -59,14 +59,14 @@ public static partial class ULongExtensions int windowLength = index - windowStart + 1; for (int square = 0; square < windowLength; square++) { - result = result.MontgomeryMultiply(result, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(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); + result = result.MontgomeryMultiplyCpu(multiplier, modulus, nPrime); index = windowStart - 1; } @@ -75,7 +75,7 @@ public static partial class ULongExtensions if (!keepMontgomery) { - result = result.MontgomeryMultiply(1UL, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); } return result; @@ -96,7 +96,7 @@ public static partial class ULongExtensions for (int i = 1; i < oddPowers.Length; i++) { previous = oddPowers[i - 1]; - oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); + oddPowers[i] = previous.MontgomeryMultiplyCpu(square, modulus, nPrime); } } diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs index 085d41c..d4fc788 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -46,6 +46,45 @@ public static partial class ULongExtensions return state.MulMod(right, modulus); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong MontgomeryMultiplyGpu(this ulong a, ulong b, ulong modulus, ulong nPrime) + { + // Mirror the struct-based MontMul64 helper measured as the second-fastest GPU-compatible + // option in GpuUInt128Montgomery64Benchmarks so accelerator kernels reuse the validated path. + MultiplyPartsGpu(a, b, out ulong productHigh, out ulong productLow); + MultiplyPartsGpu(productLow, nPrime, out _, out ulong mLow); + MultiplyPartsGpu(mLow, modulus, out ulong mTimesModulusHigh, out ulong mTimesModulusLow); + + ulong sumLow = unchecked(productLow + mTimesModulusLow); + ulong carry = sumLow < productLow ? 1UL : 0UL; + ulong sumHigh = unchecked(productHigh + mTimesModulusHigh + carry); + + if (sumHigh >= modulus) + { + sumHigh -= modulus; + } + + return sumHigh; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void MultiplyPartsGpu(ulong left, ulong right, out ulong high, out ulong low) + { + ulong leftLow = (uint)left; + ulong leftHigh = left >> 32; + ulong rightLow = (uint)right; + ulong rightHigh = right >> 32; + + ulong lowProduct = unchecked(leftLow * rightLow); + ulong cross1 = unchecked(leftHigh * rightLow); + ulong cross2 = unchecked(leftLow * rightHigh); + ulong highProduct = unchecked(leftHigh * rightHigh); + + ulong carry = unchecked((lowProduct >> 32) + (uint)cross1 + (uint)cross2); + low = unchecked((lowProduct & 0xFFFFFFFFUL) | (carry << 32)); + high = unchecked(highProduct + (cross1 >> 32) + (cross2 >> 32) + (carry >> 32)); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) { @@ -204,14 +243,15 @@ public static partial class ULongExtensions ulong modulus = divisor.Modulus; ulong nPrime = divisor.NPrime; ulong montgomeryResult = Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent); - return MontgomeryMultiply(montgomeryResult, 1UL, modulus, nPrime); + return MontgomeryMultiplyGpu(montgomeryResult, 1UL, modulus, nPrime); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorData divisor, ulong exponent) { ulong modulus = divisor.Modulus; - // Reuse the optimized MontgomeryMultiply extension highlighted in MontgomeryMultiplyBenchmarks so the GPU path matches the fastest CPU reduction. + // Reuse the GPU-compatible MontgomeryMultiply helper highlighted in MontgomeryMultiplyBenchmarks so the accelerator path + // mirrors the fastest CPU reduction without relying on UInt128 intrinsics. // We shouldn't hit it in production code. If we do, we're doing something wrong. // if (exponent == 0UL) // { @@ -233,7 +273,7 @@ public static partial class ULongExtensions while (index >= 0) { ulong currentBit = (exponent >> index) & 1UL; - ulong squared = MontgomeryMultiply(result, result, modulus, nPrime); + ulong squared = MontgomeryMultiplyGpu(result, result, modulus, nPrime); bool processWindow = currentBit != 0UL; result = processWindow ? result : squared; @@ -253,11 +293,11 @@ public static partial class ULongExtensions int windowLength = processWindow ? index - windowStart + 1 : 0; for (int square = 0; square < windowLength; square++) { - result = MontgomeryMultiply(result, result, modulus, nPrime); + result = MontgomeryMultiplyGpu(result, result, modulus, nPrime); } result = processWindow - ? MontgomeryMultiply( + ? MontgomeryMultiplyGpu( result, ComputeMontgomeryOddPowerGpu( (exponent >> windowStart) & ((1UL << windowLength) - 1UL), @@ -330,7 +370,7 @@ public static partial class ULongExtensions while (remaining != 0UL) { ulong bit = remaining & 1UL; - ulong multiplied = MontgomeryMultiply(power, baseValue, modulus, nPrime); + ulong multiplied = MontgomeryMultiplyGpu(power, baseValue, modulus, nPrime); ulong mask = (ulong)-(long)bit; power = (power & ~mask) | (multiplied & mask); @@ -340,7 +380,7 @@ public static partial class ULongExtensions break; } - baseValue = MontgomeryMultiply(baseValue, baseValue, modulus, nPrime); + baseValue = MontgomeryMultiplyGpu(baseValue, baseValue, modulus, nPrime); } return power; diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 4c057ef..a0726e8 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -553,7 +553,7 @@ public static partial class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime) + public static ulong MontgomeryMultiplyCpu(this ulong a, ulong b, ulong modulus, ulong nPrime) { ulong tLow = unchecked(a * b); ulong m = unchecked(tLow * nPrime); @@ -581,7 +581,7 @@ public static partial class ULongExtensions { if ((remainingExponent & 1UL) != 0UL) { - result = result.MontgomeryMultiply(baseVal, modulus, nPrime); + result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime); } remainingExponent >>= 1; @@ -590,10 +590,10 @@ public static partial class ULongExtensions break; } - baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, nPrime); + baseVal = baseVal.MontgomeryMultiplyCpu(baseVal, modulus, nPrime); } - return keepMontgomery ? result : result.MontgomeryMultiply(1UL, modulus, nPrime); + return keepMontgomery ? result : result.MontgomeryMultiplyCpu(1UL, modulus, nPrime); }