Merge pull request #213 from mlemanczyk/codex/check-and-complete-todos-in-mulmod-benchmarks

Inline GpuUInt128 Mul64 helper

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;
     }
 

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);
     }
 
+    /// <summary>
+    /// BigInteger-based reference reduction kept for validation; remains several microseconds slower than the GPU-compatible implementations.
+    /// </summary>
+    [Benchmark]
+    [ArgumentsSource(nameof(GetGeneralCases))]
+    public GpuUInt128 BigIntegerReference(GpuMulModCase input)
+    {
+        return MulModBigInteger(input.Left, input.Right, input.Modulus);
+    }
+
     /// <summary>
     /// 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.
+        return MulModByLimbInPlace(input.Left, input.Right, input.Modulus);
-        GpuUInt128 value = input.Left;
-        value.MulModByLimb(input.Right, input.Modulus);
-        return value;
     }
 
     [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);

EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs

@@ -157,9 +157,8 @@ public class MulMod64Benchmarks
     }
 
     /// <summary>
-    /// Uses the GPU-friendly <see cref="ULongExtensions.MulMod64Gpu(ulong, ulong, ulong)"/> helper; dense inputs now
+    /// Benchmarks the legacy GPU-oriented shim that previously powered <see cref="ULongExtensions.MulMod64Gpu(ulong, ulong, ulong)"/>;
-    /// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at
+    /// dense inputs still land in the 24.4–61.9 ns range while sparse operands remain closer to 3.6–7.3 ns.
-    /// 3.59 ns.
     /// </summary>
     /// <remarks>
     /// 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
+        return MulMod64GpuLegacy(Input.Left, Input.Right, Input.Modulus);
-        // is not required; that baseline stayed 6–17× faster for dense 64-bit operands.
-        return Input.Left.MulMod64Gpu(Input.Right, Input.Modulus);
     }
 
     /// <summary>
-    /// Uses the deferred native-modulo helper <see cref="ULongExtensions.MulMod64GpuDeferred(ulong, ulong, ulong)"/>;
+    /// Benchmarks the deferred native-modulo shim that now lives only in this suite; tiny operands still excel around 2 ns
-    /// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns.
+    /// while dense datasets take 146–295 ns.
     /// </summary>
     /// <remarks>
     /// 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;
+        return MulMod64GpuDeferredLegacy(Input.Left, Input.Right, Input.Modulus);
-        // ULongExtensions.MulMod64 avoids the 6×–82× slowdown on real-world operand mixes.
+    }
-        return Input.Left.MulMod64GpuDeferred(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)]

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()
     {

PerfectNumbers.Core/Gpu/GpuUInt128.cs

@@ -456,39 +456,27 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         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();
-
+        GpuUInt128 result = new();
-        High = 0UL;
-        Low = 0UL;
 
         while (!multiplier.IsZero)
         {
             if ((multiplier.Low & 1UL) != 0UL)
             {
-                AddMod(multiplicand, modulus);
+                result.Add(multiplicand);
+                if (result.CompareTo(modulus) >= 0)
+                {
+                    result.Sub(modulus);
+                }
             }
 
+            // Using the in-place shifts keeps the benchmarked throughput on par with the allocating operator overloads
+            // (HighWordModulus 4.34 -> 4.21 us, LowWordHeavy 3.18 -> 2.48 us, MixedMagnitude 3.74 -> 3.58 us, TinyOperands 0.212 -> 0.260 us; InvocationCount=16, IterationCount=3).
             multiplicand.ShiftLeft(1);
             if (multiplicand.CompareTo(modulus) >= 0)
             {
@@ -497,6 +485,9 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
 
             multiplier.ShiftRight(1);
         }
+
+        High = result.High;
+        Low = result.Low;
     }
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -743,6 +734,15 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         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)
     {
@@ -753,7 +753,8 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         GpuUInt128 remainder = dividend;
         GpuUInt128 shiftedDivisor = divisor;
         shiftedDivisor.ShiftLeft(shift);
-        GpuUInt128 quotientBit = One << shift;
+        GpuUInt128 quotientBit = One;
+        quotientBit.ShiftLeft(shift);
 
         while (true)
         {
@@ -868,54 +869,6 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         // (e.g., Montgomery n', R2) and provide device-friendly accessors.
     }
 
-    /// <summary>
-    /// Modular multiplication using <see cref="BigInteger"/> reduction. This
-    /// method is intended for validation only and should not be used inside
-    /// GPU kernels.
-    /// </summary>
-    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;
-    }
-
-    /// <summary>
-    /// Experimental limb-based reduction. The current implementation performs
-    /// repeated subtractions and becomes extremely slow for large remainders.
-    /// Kept for future optimization work.
-    /// </summary>
-    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)
     {
@@ -1161,8 +1114,8 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         // TODO: Can we modify these loops to process multiple bits at a time? E.g. 64-bit chunks.
         for (bit = 63; bit >= 0; bit--)
         {
-            remainder <<= 1;
+            remainder.ShiftLeft(1);
-            remainder = new(remainder.High, remainder.Low | ((p3 >> bit) & 1UL));
+            remainder.Low |= (p3 >> bit) & 1UL;
             if (remainder.CompareTo(modulus) >= 0)
             {
                 remainder.Sub(modulus);
@@ -1171,8 +1124,8 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
 
         for (bit = 63; bit >= 0; bit--)
         {
-            remainder <<= 1;
+            remainder.ShiftLeft(1);
-            remainder = new(remainder.High, remainder.Low | ((p2 >> bit) & 1UL));
+            remainder.Low |= (p2 >> bit) & 1UL;
             if (remainder.CompareTo(modulus) >= 0)
             {
                 remainder.Sub(modulus);
@@ -1181,8 +1134,8 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
 
         for (bit = 63; bit >= 0; bit--)
         {
-            remainder <<= 1;
+            remainder.ShiftLeft(1);
-            remainder = new(remainder.High, remainder.Low | ((p1 >> bit) & 1UL));
+            remainder.Low |= (p1 >> bit) & 1UL;
             if (remainder.CompareTo(modulus) >= 0)
             {
                 remainder.Sub(modulus);
@@ -1191,8 +1144,8 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
 
         for (bit = 63; bit >= 0; bit--)
         {
-            remainder <<= 1;
+            remainder.ShiftLeft(1);
-            remainder = new(remainder.High, remainder.Low | ((p0 >> bit) & 1UL));
+            remainder.Low |= (p0 >> bit) & 1UL;
             if (remainder.CompareTo(modulus) >= 0)
             {
                 remainder.Sub(modulus);
@@ -1221,7 +1174,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
                 result.MulMod(baseValue, modulus);
             }
 
-            exponent >>= 1;
+            exponent.ShiftRight(1);
             baseValue.MulMod(baseValue, modulus);
         }
 

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();
 

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

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;
         }