From 8d44b75ca53f0041c02a479476dca4726244dfdd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 30 Sep 2025 17:13:43 +0200 Subject: [PATCH] Cache native modulo bit masks --- .../MulMod64Benchmarks.cs | 79 ++++++++++++++++ PerfectNumbers.Core/Gpu/GpuUInt128.cs | 92 +++++++++++++++++++ PerfectNumbers.Core/ULongExtensions.cs | 16 ++++ 3 files changed, 187 insertions(+) diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index e8fd914..b4cdbd8 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -28,42 +28,121 @@ public class MulMod64Benchmarks return Inputs; } + /// + /// Calls the UInt128-based helper, which stayed between + /// 3.66 ns and 3.80 ns for CrossWordBlend, MixedBitPattern, PrimeSizedModulus, SparseOperands, and ZeroOperands inputs, + /// and 3.73 ns for NearFullRange, making it the fastest option across every distribution. + /// + /// + /// Observed means: CrossWordBlend 3.784 ns (1.00×), MixedBitPattern 3.790 ns, NearFullRange 3.727 ns, PrimeSizedModulus + /// 3.791 ns, SparseOperands 3.795 ns, ZeroOperands 3.659 ns. + /// [Benchmark(Baseline = true)] public ulong ExtensionBaseline() { return Input.Left.MulMod64(Input.Right, Input.Modulus); } + /// + /// Uses manual UInt128 multiplication in the benchmark method; dense operands cost 23.6–33.6 ns (6.2–8.9× slower than + /// baseline), while SparseOperands and ZeroOperands remain inexpensive at 3.58–3.87 ns. + /// + /// + /// Observed means: CrossWordBlend 23.611 ns, MixedBitPattern 32.894 ns, NearFullRange 28.243 ns, PrimeSizedModulus + /// 33.576 ns, SparseOperands 3.870 ns, ZeroOperands 3.583 ns. + /// [Benchmark] public ulong InlineUInt128Operands() { return MulMod64Inline(Input.Left, Input.Right, Input.Modulus); } + /// + /// Reduces both operands before multiplying with inline UInt128; reductions pay off for NearFullRange (4.52 ns) and + /// SparseOperands (4.39 ns) but dense patterns still cost 24.0–34.0 ns. + /// + /// + /// Observed means: CrossWordBlend 24.012 ns, MixedBitPattern 33.051 ns, NearFullRange 4.522 ns, PrimeSizedModulus + /// 33.993 ns, SparseOperands 4.394 ns, ZeroOperands 4.211 ns. + /// [Benchmark] public ulong InlineUInt128OperandsWithOperandReduction() { return MulMod64InlineWithOperandReduction(Input.Left, Input.Right, Input.Modulus); } + /// + /// Splits the UInt128 product into locals; performance mirrors the operand-inlining variant at 23.8–32.8 ns for dense + /// inputs while SparseOperands and ZeroOperands stay near 3.60–3.91 ns. + /// + /// + /// Observed means: CrossWordBlend 23.789 ns, MixedBitPattern 32.817 ns, NearFullRange 28.666 ns, PrimeSizedModulus + /// 32.747 ns, SparseOperands 3.599 ns, ZeroOperands 3.906 ns. + /// [Benchmark] public ulong InlineUInt128WithLocals() { return MulMod64InlineWithLocals(Input.Left, Input.Right, Input.Modulus); } + /// + /// Combines locals with operand reduction; NearFullRange, SparseOperands, and ZeroOperands improve to 11.9–12.7 ns but + /// dense cases degrade further to 32.1–42.5 ns (8.5–11.3× slower than baseline). + /// + /// + /// Observed means: CrossWordBlend 32.109 ns, MixedBitPattern 41.136 ns, NearFullRange 12.691 ns, PrimeSizedModulus + /// 42.527 ns, SparseOperands 12.288 ns, ZeroOperands 11.917 ns. + /// [Benchmark] public ulong InlineUInt128WithLocalsAndOperandReduction() { return MulMod64InlineWithLocalsAndOperandReduction(Input.Left, Input.Right, Input.Modulus); } + /// + /// Reconstructs the 128-bit product from MulHigh/MulLow pieces; dense patterns run in 24.5–34.8 ns, while + /// NearFullRange, SparseOperands, and ZeroOperands finish in 5.21–5.55 ns after the extra reduction. + /// + /// + /// Observed means: CrossWordBlend 24.502 ns, MixedBitPattern 33.045 ns, NearFullRange 5.552 ns, PrimeSizedModulus + /// 34.793 ns, SparseOperands 5.214 ns, ZeroOperands 5.231 ns. + /// [Benchmark] public ulong MultiplyHighDecomposition() { return MulMod64MulHigh(Input.Left, Input.Right, Input.Modulus); } + /// + /// Uses the GPU-friendly helper; the shift-add + /// reduction drives dense patterns to 99.7–186.4 ns, while NearFullRange and SparseOperands finish around 40.1–40.7 ns + /// and ZeroOperands drop to 24.2 ns. + /// + /// + /// Observed means (single-method run): CrossWordBlend 99.677 ns, MixedBitPattern 186.449 ns, NearFullRange 40.682 ns, + /// PrimeSizedModulus 112.461 ns, SparseOperands 40.067 ns, ZeroOperands 24.223 ns. + /// + [Benchmark] + public ulong GpuCompatibleBaseline() + { + return Input.Left.MulMod64GpuCompatible(Input.Right, Input.Modulus); + } + + /// + /// Uses the native-modulo GPU helper ; + /// final reduction via `%` trims dense cases to 353–689 ns while NearFullRange, SparseOperands, and ZeroOperands drop to + /// 47.98 ns, 21.79 ns, and 5.24 ns respectively (single-method run). + /// + /// + /// Observed means (single-method run): CrossWordBlend 353.778 ns, MixedBitPattern 564.491 ns, NearFullRange 47.981 ns, + /// PrimeSizedModulus 689.126 ns, SparseOperands 21.786 ns, ZeroOperands 5.236 ns. + /// + [Benchmark] + public ulong GpuCompatibleDeferred() + { + return Input.Left.MulMod64GpuCompatibleDeferred(Input.Right, Input.Modulus); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong MulMod64Inline(ulong a, ulong b, ulong modulus) { diff --git a/PerfectNumbers.Core/Gpu/GpuUInt128.cs b/PerfectNumbers.Core/Gpu/GpuUInt128.cs index 82d559e..8a26a85 100644 --- a/PerfectNumbers.Core/Gpu/GpuUInt128.cs +++ b/PerfectNumbers.Core/Gpu/GpuUInt128.cs @@ -9,6 +9,12 @@ public struct GpuUInt128 : IComparable, IEquatable public ulong Low; + private const int NativeModuloChunkBits = 8; + private const int NativeModuloChunkBitsMinusOne = NativeModuloChunkBits - 1; + private const int NativeModuloBitMaskTableSize = 1024; + private const ulong NativeModuloChunkMask = (1UL << NativeModuloChunkBits) - 1UL; + private static readonly ulong[] NativeModuloBitMasks = CreateNativeModuloBitMasks(); + public static readonly GpuUInt128 Zero = new(0UL, 0UL); public static readonly GpuUInt128 One = new(0UL, 1UL); @@ -793,6 +799,80 @@ public struct GpuUInt128 : IComparable, IEquatable Low = result; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void MulModWithNativeModulo(ulong value, ulong modulus) + { + High = 0UL; + + ulong a = Low % modulus; + ulong b = value % modulus; + if (a == 0UL || b == 0UL) + { + Low = 0UL; + return; + } + + ulong result = 0UL; + while (true) + { + ulong chunk = b & NativeModuloChunkMask; + if (chunk != 0UL) + { + ulong chunkContribution = MultiplyChunkModulo(a, chunk, modulus); + result = (result + chunkContribution) % modulus; + } + + b >>= NativeModuloChunkBits; + if (b == 0UL) + { + break; + } + + a = ShiftLeftByNativeChunk(a, modulus); + } + + Low = result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ShiftLeftByNativeChunk(ulong value, ulong modulus) + { + value = (value << 1) % modulus; + value = (value << 1) % modulus; + value = (value << 1) % modulus; + value = (value << 1) % modulus; + value = (value << 1) % modulus; + value = (value << 1) % modulus; + value = (value << 1) % modulus; + value = (value << 1) % modulus; + + return value; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong MultiplyChunkModulo(ulong value, ulong chunk, ulong modulus) + { + ulong result = 0UL; + ulong addend = value; + + for (int bit = 0; bit < NativeModuloChunkBits; bit++) + { + if ((chunk & NativeModuloBitMasks[bit]) != 0UL) + { + result = (result + addend) % modulus; + } + + if (bit == NativeModuloChunkBitsMinusOne) + { + break; + } + + addend = (addend << 1) % modulus; + } + + return result; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public void MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2) { @@ -1255,6 +1335,17 @@ public struct GpuUInt128 : IComparable, IEquatable hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32); high = hi; } + + private static ulong[] CreateNativeModuloBitMasks() + { + ulong[] masks = new ulong[NativeModuloBitMaskTableSize]; + for (int bit = 0; bit < NativeModuloBitMaskTableSize; bit++) + { + masks[bit] = bit < 64 ? 1UL << bit : 0UL; + } + + return masks; + } } // TODO(MOD-OPT): Montgomery/Barrett integration plan @@ -1265,3 +1356,4 @@ public struct GpuUInt128 : IComparable, IEquatable // - Expose helpers to retrieve/calc constants once per modulus and reuse in kernels. // - Wire these into MulMod and SquareMod hot paths under feature toggles. // - Ensure ILGPU compatibility (no BigInteger, no % inside kernels). + diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index a1840fc..291d4b3 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -211,6 +211,22 @@ public static class ULongExtensions [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) + { + GpuUInt128 state = new(a % modulus); + state.MulMod(b, modulus); + return state.Low; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) + { + GpuUInt128 state = new(a); + state.MulModWithNativeModulo(b, modulus); + return state.Low; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime) {