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

Merge pull request #247 from mlemanczyk/codex/fix-unhandled-exception-in-evenperfectbitscanner

Browse Source 
Fix GPU Montgomery multiply to avoid unsupported intrinsics
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-10-28 21:10:09 +01:00
2 parents 4db06cfe80 e5f5cd563e
commit b83bddf49a
9 files changed
+108 -70
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)] [Benchmark(Baseline = true)]
public ulong ExtensionMontgomeryMultiply() public ulong ExtensionMontgomeryMultiply()
{ {
return Input.Left.MontgomeryMultiply(Input.Right, Input.Modulus, Input.NPrime); return Input.Left.MontgomeryMultiplyCpu(Input.Right, Input.Modulus, Input.NPrime);
} }
/// <summary> /// <summary>
EvenPerfectBitScanner.Benchmarks/Pow2ModBenchmarks.cs
+21 -21
View File
@@ -640,7 +640,7 @@ public class Pow2ModBenchmarks
ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus); ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, 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( return new MontgomeryDivisorData(
modulus, modulus,
@@ -736,22 +736,22 @@ public class Pow2ModBenchmarks
if (exponent == 0UL) if (exponent == 0UL)
{ {
return result.MontgomeryMultiply(1UL, modulus, nPrime); return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent); int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent);
for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--) for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--)
{ {
result = result.MontgomeryMultiply(result, modulus, nPrime); result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
if (((exponent >> bitIndex) & 1UL) != 0UL) 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) private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor)
@@ -772,35 +772,35 @@ public class Pow2ModBenchmarks
ulong nibble = remainingExponent & 0xFUL; ulong nibble = remainingExponent & 0xFUL;
ulong base0 = baseVal; ulong base0 = baseVal;
ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime); ulong base1 = base0.MontgomeryMultiplyCpu(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime); ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime); ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime);
if ((nibble & 1UL) != 0UL) if ((nibble & 1UL) != 0UL)
{ {
result = result.MontgomeryMultiply(base0, modulus, nPrime); result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime);
} }
if ((nibble & 2UL) != 0UL) if ((nibble & 2UL) != 0UL)
{ {
result = result.MontgomeryMultiply(base1, modulus, nPrime); result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime);
} }
if ((nibble & 4UL) != 0UL) if ((nibble & 4UL) != 0UL)
{ {
result = result.MontgomeryMultiply(base2, modulus, nPrime); result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime);
} }
if ((nibble & 8UL) != 0UL) 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; remainingExponent >>= 4;
} }
return result.MontgomeryMultiply(1UL, modulus, nPrime); return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MontgomeryDivisorData divisor) private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MontgomeryDivisorData divisor)
@@ -817,16 +817,16 @@ public class Pow2ModBenchmarks
if (exponent == 0UL) if (exponent == 0UL)
{ {
return result.MontgomeryMultiply(1UL, modulus, nPrime); return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
const int WindowSize = 5; const int WindowSize = 5;
Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)]; Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)];
oddPowers[0] = baseVal; 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++) 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); int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent);
@@ -835,7 +835,7 @@ public class Pow2ModBenchmarks
{ {
if (((exponent >> bitIndex) & 1UL) == 0UL) if (((exponent >> bitIndex) & 1UL) == 0UL)
{ {
result = result.MontgomeryMultiply(result, modulus, nPrime); result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
bitIndex--; bitIndex--;
continue; continue;
} }
@@ -852,14 +852,14 @@ public class Pow2ModBenchmarks
int squares = bitIndex - windowStart + 1; int squares = bitIndex - windowStart + 1;
for (int i = 0; i < squares; i++) 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; bitIndex = windowStart - 1;
} }
return result.MontgomeryMultiply(1UL, modulus, nPrime); return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
private static ulong Pow2ModBinary(ulong exponent, ulong modulus) 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 nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus); ulong montgomeryOne = ComputeMontgomeryResidue(1UL, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(2UL, 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( return new MontgomeryDivisorData(
modulus, modulus,
@@ -300,22 +300,22 @@ public class Pow2MontgomeryModBenchmarks
if (exponent == 0UL) if (exponent == 0UL)
{ {
return result.MontgomeryMultiply(1UL, modulus, nPrime); return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent); int msbIndex = 63 - BitOperations.LeadingZeroCount(exponent);
for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--) for (int bitIndex = msbIndex; bitIndex >= 0; bitIndex--)
{ {
result = result.MontgomeryMultiply(result, modulus, nPrime); result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
if (((exponent >> bitIndex) & 1UL) != 0UL) 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) private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor)
@@ -336,35 +336,35 @@ public class Pow2MontgomeryModBenchmarks
ulong nibble = remainingExponent & 0xFUL; ulong nibble = remainingExponent & 0xFUL;
ulong base0 = baseVal; ulong base0 = baseVal;
ulong base1 = base0.MontgomeryMultiply(base0, modulus, nPrime); ulong base1 = base0.MontgomeryMultiplyCpu(base0, modulus, nPrime);
ulong base2 = base1.MontgomeryMultiply(base1, modulus, nPrime); ulong base2 = base1.MontgomeryMultiplyCpu(base1, modulus, nPrime);
ulong base3 = base2.MontgomeryMultiply(base2, modulus, nPrime); ulong base3 = base2.MontgomeryMultiplyCpu(base2, modulus, nPrime);
if ((nibble & 1UL) != 0UL) if ((nibble & 1UL) != 0UL)
{ {
result = result.MontgomeryMultiply(base0, modulus, nPrime); result = result.MontgomeryMultiplyCpu(base0, modulus, nPrime);
} }
if ((nibble & 2UL) != 0UL) if ((nibble & 2UL) != 0UL)
{ {
result = result.MontgomeryMultiply(base1, modulus, nPrime); result = result.MontgomeryMultiplyCpu(base1, modulus, nPrime);
} }
if ((nibble & 4UL) != 0UL) if ((nibble & 4UL) != 0UL)
{ {
result = result.MontgomeryMultiply(base2, modulus, nPrime); result = result.MontgomeryMultiplyCpu(base2, modulus, nPrime);
} }
if ((nibble & 8UL) != 0UL) 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; remainingExponent >>= 4;
} }
return result.MontgomeryMultiply(1UL, modulus, nPrime); return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
// TODO: Update it with the latest implementations // TODO: Update it with the latest implementations
@@ -382,16 +382,16 @@ public class Pow2MontgomeryModBenchmarks
if (exponent == 0UL) if (exponent == 0UL)
{ {
return result.MontgomeryMultiply(1UL, modulus, nPrime); return result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
const int WindowSize = 5; const int WindowSize = 5;
Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)]; Span<ulong> oddPowers = stackalloc ulong[1 << (WindowSize - 1)];
oddPowers[0] = baseVal; 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++) 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); int bitIndex = 63 - BitOperations.LeadingZeroCount(exponent);
@@ -400,7 +400,7 @@ public class Pow2MontgomeryModBenchmarks
{ {
if (((exponent >> bitIndex) & 1UL) == 0UL) if (((exponent >> bitIndex) & 1UL) == 0UL)
{ {
result = result.MontgomeryMultiply(result, modulus, nPrime); result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
bitIndex--; bitIndex--;
continue; continue;
} }
@@ -417,14 +417,14 @@ public class Pow2MontgomeryModBenchmarks
int squares = bitIndex - windowStart + 1; int squares = bitIndex - windowStart + 1;
for (int i = 0; i < squares; i++) 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; 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 // 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 // single-block divisor-cycle snapshot so we can skip the powmod entirely and reuse the
// cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks. // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
PreviousExponent = exponent; PreviousExponent = exponent;
return ReduceCurrent(); return ReduceCurrent();
} }
@@ -92,7 +92,7 @@ internal struct ExponentRemainderStepperCpu
// TODO: Reuse the divisor-cycle derived Montgomery delta once the cache exposes single-cycle // 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. // lookups so this branch also avoids recomputing powmods when the snapshot lacks the divisor.
ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(multiplier, _modulus, _nPrime);
PreviousExponent = exponent; PreviousExponent = exponent;
return _currentMontgomery == _montgomeryOne; 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 // 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 // of relying on the caller-provided delta so incremental scans remain in sync with the snapshot
// without computing additional cycles or mutating cache state. // without computing additional cycles or mutating cache state.
_currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime); _currentMontgomery = _currentMontgomery.MontgomeryMultiplyCpu(montgomeryDelta, _modulus, _nPrime);
PreviousExponent = exponent; PreviousExponent = exponent;
isUnity = _currentMontgomery == _montgomeryOne; isUnity = _currentMontgomery == _montgomeryOne;
return true; return true;
@@ -138,5 +138,5 @@ internal struct ExponentRemainderStepperCpu
PreviousExponent = exponent; 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 nPrime = ComputeMontgomeryNPrime(modulus);
ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus); ulong montgomeryOne = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64, modulus);
ulong montgomeryTwo = ComputeMontgomeryResidue(UInt128Numbers.OneShiftedLeft64x2, 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( return new MontgomeryDivisorData(
modulus, modulus,
PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs
+4 -6
View File
@@ -53,8 +53,7 @@ internal static partial class PrimeOrderCalculator
ulong phi = prime - 1UL; ulong phi = prime - 1UL;
// if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder))
if (PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder))
{ {
return gpuOrder; return gpuOrder;
} }
@@ -1070,7 +1069,7 @@ internal static partial class PrimeOrderCalculator
ulong remaining = value; ulong remaining = value;
bool gpuFactored = false; bool gpuFactored = false;
// if (IsGpuHeuristicDevice) if (IsGpuHeuristicDevice)
{ {
gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _);
} }
@@ -1079,7 +1078,6 @@ internal static partial class PrimeOrderCalculator
Stack<ulong>? compositeStack = null; Stack<ulong>? compositeStack = null;
PartialFactorResult result; PartialFactorResult result;
PrimeTester primeTester = PrimeTester.Exclusive; PrimeTester primeTester = PrimeTester.Exclusive;
CancellationToken ct = CancellationToken.None;
if (!gpuFactored) if (!gpuFactored)
{ {
@@ -1152,7 +1150,7 @@ internal static partial class PrimeOrderCalculator
// Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}"); // Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}");
// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); // 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); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
if (isPrime) if (isPrime)
@@ -1220,7 +1218,7 @@ internal static partial class PrimeOrderCalculator
// Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}"); // Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}");
// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); // 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); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
entry = entry.WithPrimality(isPrime); entry = entry.WithPrimality(isPrime);
PerfectNumbers.Core/ULongExtensions.Cpu.cs
+5 -5
View File
@@ -40,7 +40,7 @@ public static partial class ULongExtensions
{ {
if (((exponent >> index) & 1UL) == 0UL) if (((exponent >> index) & 1UL) == 0UL)
{ {
result = result.MontgomeryMultiply(result, modulus, nPrime); result = result.MontgomeryMultiplyCpu(result, modulus, nPrime);
index--; index--;
continue; continue;
} }
@@ -59,14 +59,14 @@ public static partial class ULongExtensions
int windowLength = index - windowStart + 1; int windowLength = index - windowStart + 1;
for (int square = 0; square < windowLength; square++) 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 mask = (1UL << windowLength) - 1UL;
ulong windowValue = (exponent >> windowStart) & mask; ulong windowValue = (exponent >> windowStart) & mask;
int tableIndex = (int)((windowValue - 1UL) >> 1); int tableIndex = (int)((windowValue - 1UL) >> 1);
ulong multiplier = oddPowers[tableIndex]; ulong multiplier = oddPowers[tableIndex];
result = result.MontgomeryMultiply(multiplier, modulus, nPrime); result = result.MontgomeryMultiplyCpu(multiplier, modulus, nPrime);
index = windowStart - 1; index = windowStart - 1;
} }
@@ -75,7 +75,7 @@ public static partial class ULongExtensions
if (!keepMontgomery) if (!keepMontgomery)
{ {
result = result.MontgomeryMultiply(1UL, modulus, nPrime); result = result.MontgomeryMultiplyCpu(1UL, modulus, nPrime);
} }
return result; return result;
@@ -96,7 +96,7 @@ public static partial class ULongExtensions
for (int i = 1; i < oddPowers.Length; i++) for (int i = 1; i < oddPowers.Length; i++)
{ {
previous = oddPowers[i - 1]; previous = oddPowers[i - 1];
oddPowers[i] = previous.MontgomeryMultiply(square, modulus, nPrime); oddPowers[i] = previous.MontgomeryMultiplyCpu(square, modulus, nPrime);
} }
} }
PerfectNumbers.Core/ULongExtensions.Gpu.cs
+47 -7
View File
@@ -46,6 +46,45 @@ public static partial class ULongExtensions
return state.MulMod(right, modulus); 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)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus) public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus)
{ {
@@ -204,14 +243,15 @@ public static partial class ULongExtensions
ulong modulus = divisor.Modulus; ulong modulus = divisor.Modulus;
ulong nPrime = divisor.NPrime; ulong nPrime = divisor.NPrime;
ulong montgomeryResult = Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent); ulong montgomeryResult = Pow2MontgomeryModWindowedGpuMontgomeryResult(divisor, exponent);
return MontgomeryMultiply(montgomeryResult, 1UL, modulus, nPrime); return MontgomeryMultiplyGpu(montgomeryResult, 1UL, modulus, nPrime);
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorData divisor, ulong exponent) private static ulong Pow2MontgomeryModWindowedGpuMontgomeryResult(in MontgomeryDivisorData divisor, ulong exponent)
{ {
ulong modulus = divisor.Modulus; 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. // We shouldn't hit it in production code. If we do, we're doing something wrong.
// if (exponent == 0UL) // if (exponent == 0UL)
// { // {
@@ -233,7 +273,7 @@ public static partial class ULongExtensions
while (index >= 0) while (index >= 0)
{ {
ulong currentBit = (exponent >> index) & 1UL; ulong currentBit = (exponent >> index) & 1UL;
ulong squared = MontgomeryMultiply(result, result, modulus, nPrime); ulong squared = MontgomeryMultiplyGpu(result, result, modulus, nPrime);
bool processWindow = currentBit != 0UL; bool processWindow = currentBit != 0UL;
result = processWindow ? result : squared; result = processWindow ? result : squared;
@@ -253,11 +293,11 @@ public static partial class ULongExtensions
int windowLength = processWindow ? index - windowStart + 1 : 0; int windowLength = processWindow ? index - windowStart + 1 : 0;
for (int square = 0; square < windowLength; square++) for (int square = 0; square < windowLength; square++)
{ {
result = MontgomeryMultiply(result, result, modulus, nPrime); result = MontgomeryMultiplyGpu(result, result, modulus, nPrime);
} }
result = processWindow result = processWindow
? MontgomeryMultiply( ? MontgomeryMultiplyGpu(
result, result,
ComputeMontgomeryOddPowerGpu( ComputeMontgomeryOddPowerGpu(
(exponent >> windowStart) & ((1UL << windowLength) - 1UL), (exponent >> windowStart) & ((1UL << windowLength) - 1UL),
@@ -330,7 +370,7 @@ public static partial class ULongExtensions
while (remaining != 0UL) while (remaining != 0UL)
{ {
ulong bit = remaining & 1UL; ulong bit = remaining & 1UL;
ulong multiplied = MontgomeryMultiply(power, baseValue, modulus, nPrime); ulong multiplied = MontgomeryMultiplyGpu(power, baseValue, modulus, nPrime);
ulong mask = (ulong)-(long)bit; ulong mask = (ulong)-(long)bit;
power = (power & ~mask) | (multiplied & mask); power = (power & ~mask) | (multiplied & mask);
@@ -340,7 +380,7 @@ public static partial class ULongExtensions
break; break;
} }
baseValue = MontgomeryMultiply(baseValue, baseValue, modulus, nPrime); baseValue = MontgomeryMultiplyGpu(baseValue, baseValue, modulus, nPrime);
} }
return power; return power;
PerfectNumbers.Core/ULongExtensions.cs
+4 -4
View File
@@ -553,7 +553,7 @@ public static partial class ULongExtensions
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [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 tLow = unchecked(a * b);
ulong m = unchecked(tLow * nPrime); ulong m = unchecked(tLow * nPrime);
@@ -581,7 +581,7 @@ public static partial class ULongExtensions
{ {
if ((remainingExponent & 1UL) != 0UL) if ((remainingExponent & 1UL) != 0UL)
{ {
result = result.MontgomeryMultiply(baseVal, modulus, nPrime); result = result.MontgomeryMultiplyCpu(baseVal, modulus, nPrime);
} }
remainingExponent >>= 1; remainingExponent >>= 1;
@@ -590,10 +590,10 @@ public static partial class ULongExtensions
break; 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);
} }