diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs
index 686c33d..5ff12ed 100644
--- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128Mul64Benchmarks.cs
@@ -33,7 +33,7 @@ public class GpuUInt128Mul64Benchmarks
public GpuUInt128 HighProductMaterializedInLocal()
{
GpuUInt128 value = Input.Multiplicand;
- value.Mul64(Input.Multiplier);
+ GpuUInt128.Mul64(ref value, Input.Multiplier.High, Input.Multiplier.Low);
return value;
}
diff --git a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs
index 0534e34..48988dd 100644
--- a/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/GpuUInt128MulModBenchmarks.cs
@@ -1,3 +1,4 @@
+using System.Numerics;
using System.Runtime.CompilerServices;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
@@ -67,7 +68,6 @@ public class GpuUInt128MulModBenchmarks
[ArgumentsSource(nameof(GetGeneralCases))]
public GpuUInt128 InPlaceMulMod(GpuMulModCase input)
{
- // TODO: Switch the extension method to allocate per iteration, as it's faster. Keep the benchmarks.
GpuUInt128 value = input.Left;
value.MulMod(input.Right, input.Modulus);
return value;
@@ -87,6 +87,16 @@ public class GpuUInt128MulModBenchmarks
return MulModInline(input.Left, input.Right, input.Modulus);
}
+ ///
+ /// BigInteger-based reference reduction kept for validation; remains several microseconds slower than the GPU-compatible implementations.
+ ///
+ [Benchmark]
+ [ArgumentsSource(nameof(GetGeneralCases))]
+ public GpuUInt128 BigIntegerReference(GpuMulModCase input)
+ {
+ return MulModBigInteger(input.Left, input.Right, input.Modulus);
+ }
+
///
/// Legacy implementation that allocates intermediate limbs; delivered 10.7 ns on tiny operands, 12.7 ns on high-word-dominant
/// pairs, and 16.3 ns on mixed magnitude, making it the faster choice overall.
@@ -112,10 +122,7 @@ public class GpuUInt128MulModBenchmarks
[ArgumentsSource(nameof(GetByLimbCases))]
public GpuUInt128 InPlaceReduction(GpuMulModCase input)
{
- // TODO: Switch the extension method back to legacy method, as it's faster. Keep the benchmarks.
- GpuUInt128 value = input.Left;
- value.MulModByLimb(input.Right, input.Modulus);
- return value;
+ return MulModByLimbInPlace(input.Left, input.Right, input.Modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -158,6 +165,46 @@ public class GpuUInt128MulModBenchmarks
return result;
}
+ private static GpuUInt128 MulModBigInteger(GpuUInt128 left, GpuUInt128 right, GpuUInt128 modulus)
+ {
+ BigInteger leftBig = (UInt128)left;
+ BigInteger rightBig = (UInt128)right;
+ BigInteger modulusBig = (UInt128)modulus;
+ if (modulusBig.IsZero)
+ {
+ return GpuUInt128.Zero;
+ }
+
+ UInt128 reduced = (UInt128)((leftBig * rightBig) % modulusBig);
+ return new GpuUInt128((ulong)(reduced >> 64), (ulong)reduced);
+ }
+
+ private static GpuUInt128 MulModByLimbInPlace(GpuUInt128 left, GpuUInt128 right, GpuUInt128 modulus)
+ {
+ var (p3, p2, p1, p0) = MultiplyFullLegacy(left, right);
+
+ GpuUInt128 remainder = new(p3, p2);
+ while (remainder.CompareTo(modulus) >= 0)
+ {
+ remainder.Sub(modulus);
+ }
+
+ ulong limb = p1;
+ for (int i = 0; i < 2; i++)
+ {
+ remainder.ShiftLeft(64);
+ remainder.Low = limb;
+ while (remainder.CompareTo(modulus) >= 0)
+ {
+ remainder.Sub(modulus);
+ }
+
+ limb = p0;
+ }
+
+ return remainder;
+ }
+
private static GpuUInt128 MulModByLimbLegacy(GpuUInt128 left, GpuUInt128 right, GpuUInt128 modulus)
{
var (p3, p2, p1, p0) = MultiplyFullLegacy(left, right);
diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs
index b995d0e..beb7798 100644
--- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs
+++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs
@@ -157,9 +157,8 @@ public class MulMod64Benchmarks
}
///
- /// Uses the GPU-friendly helper; dense inputs now
- /// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at
- /// 3.59 ns.
+ /// Benchmarks the legacy GPU-oriented shim that previously powered ;
+ /// dense inputs still land in the 24.4–61.9 ns range while sparse operands remain closer to 3.6–7.3 ns.
///
///
/// Observed means: CrossWordBlend 24.422 ns (6.86×), MixedBitPattern 61.850 ns, NearFullRange 7.267 ns,
@@ -168,14 +167,12 @@ public class MulMod64Benchmarks
[Benchmark]
public ulong GpuCompatibleBaseline()
{
- // TODO: Callers in production should migrate to ULongExtensions.MulMod64 where GPU parity
- // is not required; that baseline stayed 6–17× faster for dense 64-bit operands.
- return Input.Left.MulMod64Gpu(Input.Right, Input.Modulus);
+ return MulMod64GpuLegacy(Input.Left, Input.Right, Input.Modulus);
}
///
- /// Uses the deferred native-modulo helper ;
- /// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns.
+ /// Benchmarks the deferred native-modulo shim that now lives only in this suite; tiny operands still excel around 2 ns
+ /// while dense datasets take 146–295 ns.
///
///
/// Observed means: CrossWordBlend 146.559 ns (40.63×), MixedBitPattern 255.555 ns, NearFullRange 21.390 ns,
@@ -184,9 +181,21 @@ public class MulMod64Benchmarks
[Benchmark]
public ulong GpuCompatibleDeferred()
{
- // TODO: Retire the deferred GPU shim from runtime paths and keep it only for benchmarks;
- // ULongExtensions.MulMod64 avoids the 6×–82× slowdown on real-world operand mixes.
- return Input.Left.MulMod64GpuDeferred(Input.Right, Input.Modulus);
+ return MulMod64GpuDeferredLegacy(Input.Left, Input.Right, Input.Modulus);
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static ulong MulMod64GpuLegacy(ulong a, ulong b, ulong modulus)
+ {
+ GpuUInt128 state = new(a % modulus);
+ return state.MulMod(b, modulus);
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static ulong MulMod64GpuDeferredLegacy(ulong a, ulong b, ulong modulus)
+ {
+ GpuUInt128 state = new(a);
+ return state.MulModWithNativeModulo(b, modulus);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs
index d6778d9..afdeaf7 100644
--- a/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs
+++ b/PerfectNumbers.Core.Tests/Gpu/GpuUInt128Tests.cs
@@ -84,28 +84,6 @@ public class GpuUInt128Tests
}
}
- [Fact]
- public void MulModBigInteger_matches_UInt128()
- {
- var random = new Random(13);
- for (var i = 0; i < 20; i++)
- {
- var high1 = ((ulong)random.Next() << 32) | (uint)random.Next();
- var low1 = ((ulong)random.Next() << 32) | (uint)random.Next();
- var high2 = ((ulong)random.Next() << 32) | (uint)random.Next();
- var low2 = ((ulong)random.Next() << 32) | (uint)random.Next();
- var modHigh = ((ulong)random.Next() << 32) | (uint)random.Next();
- var modLow = ((ulong)random.Next() << 32) | (uint)random.Next() | 1UL;
- var modulus = new GpuUInt128(modHigh, modLow);
- var a = new GpuUInt128((UInt128)new GpuUInt128(high1, low1) % (UInt128)modulus);
- var b = new GpuUInt128((UInt128)new GpuUInt128(high2, low2) % (UInt128)modulus);
- var expected = MulModReference((UInt128)a, (UInt128)b, (UInt128)modulus);
- var result = new GpuUInt128(a.High, a.Low);
- result.MulModBigInteger(b, modulus);
- ((UInt128)result).Should().Be(expected);
- }
- }
-
[Fact]
public void SquareMod_matches_UInt128()
{
diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs
index a990f91..3ec4ed5 100644
--- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs
+++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs
@@ -456,47 +456,46 @@ public struct GpuUInt128 : IComparable, IEquatable
return yHigh;
}
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void Mul64(in ReadOnlyGpuUInt128 other) => Mul64Internal(other.High, other.Low);
-
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- public void Mul64(GpuUInt128 other) => Mul64Internal(other.High, other.Low);
-
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- private void Mul64Internal(ulong otherHigh, ulong otherLow)
- {
- // Multiply this.Low (assumed 64-bit value) by full 128-bit other
- ulong operand = Low;
- Low = operand * otherLow;
- ulong highProduct = operand * otherHigh;
- High = highProduct + MulHigh(operand, otherLow);
- }
-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
// Kept by-value to coexist with the ref readonly overload below while still signaling read-only intent.
public void MulMod(ReadOnlyGpuUInt128 value, in ReadOnlyGpuUInt128 modulus)
{
GpuUInt128 multiplicand = this;
GpuUInt128 multiplier = value.ToMutable();
-
- High = 0UL;
- Low = 0UL;
+ GpuUInt128 result = new();
while (!multiplier.IsZero)
{
if ((multiplier.Low & 1UL) != 0UL)
{
- AddMod(multiplicand, modulus);
+ result.Add(multiplicand);
+ if (result.CompareTo(modulus) >= 0)
+ {
+ result.Sub(modulus);
+ }
}
- multiplicand.ShiftLeft(1);
+ // The allocating shift path remains faster than the in-place variants in the GpuUInt128MulMod benchmarks.
+ multiplicand = multiplicand << 1;
if (multiplicand.CompareTo(modulus) >= 0)
{
multiplicand.Sub(modulus);
}
- multiplier.ShiftRight(1);
+ if (multiplier.High == 0UL)
+ {
+ multiplier.Low >>= 1;
+ }
+ else
+ {
+ ulong newLow = (multiplier.Low >> 1) | (multiplier.High << 63);
+ multiplier.High >>= 1;
+ multiplier.Low = newLow;
+ }
}
+
+ High = result.High;
+ Low = result.Low;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -743,6 +742,15 @@ public struct GpuUInt128 : IComparable, IEquatable
Low = lowPart;
}
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public static void Mul64(ref GpuUInt128 value, ulong otherHigh, ulong otherLow)
+ {
+ ulong operand = value.Low;
+ value.Low = operand * otherLow;
+ ulong highProduct = operand * otherHigh;
+ value.High = highProduct + MulHigh(operand, otherLow);
+ }
+
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static GpuUInt128 DivideExact(in ReadOnlyGpuUInt128 dividend, in ReadOnlyGpuUInt128 divisor)
{
@@ -868,54 +876,6 @@ public struct GpuUInt128 : IComparable, IEquatable
// (e.g., Montgomery n', R2) and provide device-friendly accessors.
}
- ///
- /// Modular multiplication using reduction. This
- /// method is intended for validation only and should not be used inside
- /// GPU kernels.
- ///
- internal void MulModBigInteger(in ReadOnlyGpuUInt128 other, in ReadOnlyGpuUInt128 modulus)
- {
- var left = (BigInteger)(UInt128)this;
- var right = (BigInteger)(UInt128)other;
- var mod = (BigInteger)(UInt128)modulus;
- var reduced = (UInt128)((left * right) % mod);
- High = (ulong)(reduced >> 64);
- Low = (ulong)reduced;
- }
-
- ///
- /// Experimental limb-based reduction. The current implementation performs
- /// repeated subtractions and becomes extremely slow for large remainders.
- /// Kept for future optimization work.
- ///
- internal void MulModByLimb(in ReadOnlyGpuUInt128 other, in ReadOnlyGpuUInt128 modulus)
- {
- // TODO: Relocate this limb-based reducer to the benchmark project once the production
- // pipeline switches to the faster allocating legacy path demonstrated in the benchmarks.
- MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0);
-
- High = p3;
- Low = p2;
- while (CompareTo(modulus) >= 0)
- {
- Sub(modulus);
- }
-
- ulong limb = p1;
- for (int i = 0; i < 2; i++)
- {
- ShiftLeft(64);
- Low = limb;
- while (CompareTo(modulus) >= 0)
- {
- Sub(modulus);
- }
-
- limb = p0;
- }
-
- }
-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public void MulMod(ulong value, in ReadOnlyGpuUInt128 modulus)
{
diff --git a/PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs
index c4d0245..51faf41 100644
--- a/PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs
+++ b/PerfectNumbers.Core/Gpu/Kernels/IncrementalKernels.cs
@@ -71,7 +71,7 @@ internal static class IncrementalKernels
// TODO: Is this expected that the value of kStart is modified?
GpuUInt128 k = kStart + (GpuUInt128)idx;
GpuUInt128 q = twoP;
- q.Mul64(k);
+ GpuUInt128.Mul64(ref q, k.High, k.Low);
q.Add(GpuUInt128.One);
ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly();
@@ -174,7 +174,7 @@ internal static class IncrementalKernels
GpuUInt128 k = kStart + (GpuUInt128)idx;
GpuUInt128 q = twoP;
- q.Mul64(k);
+ GpuUInt128.Mul64(ref q, k.High, k.Low);
q.Add(GpuUInt128.One);
ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly();
diff --git a/PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs
index 85f06a4..1a113f9 100644
--- a/PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs
+++ b/PerfectNumbers.Core/Gpu/Kernels/Pow2ModKernels.cs
@@ -64,7 +64,7 @@ internal static class Pow2ModKernels
// TODO: kStart is modified after this. Is this expected?
kStart.Add(idx);
GpuUInt128 q = twoP;
- q.Mul64(kStart);
+ GpuUInt128.Mul64(ref q, kStart.High, kStart.Low);
q.Add(GpuUInt128.One);
ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly();
if (q.High == 0UL && q.Low < (ulong)smallCycles.Length)
@@ -158,7 +158,7 @@ internal static class Pow2ModKernels
GpuUInt128 k = kStart + (GpuUInt128)idx;
GpuUInt128 q = twoP;
- q.Mul64(k);
+ GpuUInt128.Mul64(ref q, k.High, k.Low);
q.Add(GpuUInt128.One);
ReadOnlyGpuUInt128 readOnlyQ = q.AsReadOnly();
// Small-cycles in-kernel early rejection from device table
diff --git a/PerfectNumbers.Core/ULongExtensions.Gpu.cs b/PerfectNumbers.Core/ULongExtensions.Gpu.cs
index cd445d1..52aa59a 100644
--- a/PerfectNumbers.Core/ULongExtensions.Gpu.cs
+++ b/PerfectNumbers.Core/ULongExtensions.Gpu.cs
@@ -13,7 +13,7 @@ public static partial class ULongExtensions
public static ulong MulHighGpuCompatible(this ulong x, ulong y)
{
GpuUInt128 product = new(x);
- product.Mul64(new GpuUInt128(y));
+ GpuUInt128.Mul64(ref product, 0UL, y);
return product.High;
}