mlemanczyk/even-perfect-numbers-scanner
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Rename CPU Montgomery multiply and align GPU helper

Browse Source
This commit is contained in:
Marek Lemańczyk committed 2025-10-28 21:06:40 +01:00
1 parent eaee599f25
commit 3304335072
8 files changed
+87 -64
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

EvenPerfectBitScanner.Benchmarks/GpuUInt128Montgomery64Benchmarks.cs
+1 -1
View File
@@ -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);
}
/// <summary>
EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs
+21 -21
View File
@@ -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<ulong> 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)
EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs
+21 -21
View File
@@ -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<ulong> 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
PerfectNumbers.Core/ExponentRemainderStepperCpu.cs
+4 -4
View File
@@ -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);
}
PerfectNumbers.Core/MontgomeryDivisorData.cs
+1 -1
View File
@@ -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,
PerfectNumbers.Core/ULongExtensions.Cpu.cs
+5 -5
View File
@@ -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);
}
}
PerfectNumbers.Core/ULongExtensions.Gpu.cs
+30 -7
View File
@@ -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)]
PerfectNumbers.Core/ULongExtensions.cs
+4 -4
View File
@@ -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);
}