From 33043350724b8fd92d0bf8fc0f9fc503f77babc2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 28 Oct 2025 21:06:40 +0100 Subject: [PATCH] Rename CPU Montgomery multiply and align GPU helper --- .../GpuUInt128Montgomery64Benchmarks.cs | 2 +- .../Pow2ModBenchmarks.cs | 42 +++++++++---------- .../Pow2MontgomeryModBenchmarks.cs | 42 +++++++++---------- .../ExponentRemainderStepperCpu.cs | 8 ++-- PerfectNumbers.Core/MontgomeryDivisorData.cs | 2 +- PerfectNumbers.Core/ULongExtensions.Cpu.cs | 10 ++--- PerfectNumbers.Core/ULongExtensions.Gpu.cs | 37 ++++++++++++---- PerfectNumbers.Core/ULongExtensions.cs | 8 ++-- 8 files changed, 87 insertions(+), 64 deletions(-) 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/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 bc062c1..d4fc788 100644 --- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs +++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs @@ -49,17 +49,40 @@ public static partial class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static ulong MontgomeryMultiplyGpu(this ulong a, ulong b, ulong modulus, ulong nPrime) { - ulong tLow = unchecked(a * b); - ulong m = unchecked(tLow * nPrime); - ulong mTimesModulusLow = unchecked(m * modulus); + // 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 result = unchecked(a.MulHighGpu(b) + m.MulHighGpu(modulus) + (unchecked(tLow + mTimesModulusLow) < tLow ? 1UL : 0UL)); - if (result >= modulus) + ulong sumLow = unchecked(productLow + mTimesModulusLow); + ulong carry = sumLow < productLow ? 1UL : 0UL; + ulong sumHigh = unchecked(productHigh + mTimesModulusHigh + carry); + + if (sumHigh >= modulus) { - result -= modulus; + sumHigh -= modulus; } - return result; + 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)] 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); }