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;
     }
 
+    /// <summary>
+    /// Calls the UInt128-based <see cref="ULongExtensions.MulMod64(ulong, ulong, ulong)"/> 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.
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
     [Benchmark(Baseline = true)]
     public ulong ExtensionBaseline()
     {
         return Input.Left.MulMod64(Input.Right, Input.Modulus);
     }
 
+    /// <summary>
+    /// 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.
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
     [Benchmark]
     public ulong InlineUInt128Operands()
     {
         return MulMod64Inline(Input.Left, Input.Right, Input.Modulus);
     }
 
+    /// <summary>
+    /// 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.
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
     [Benchmark]
     public ulong InlineUInt128OperandsWithOperandReduction()
     {
         return MulMod64InlineWithOperandReduction(Input.Left, Input.Right, Input.Modulus);
     }
 
+    /// <summary>
+    /// 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.
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
     [Benchmark]
     public ulong InlineUInt128WithLocals()
     {
         return MulMod64InlineWithLocals(Input.Left, Input.Right, Input.Modulus);
     }
 
+    /// <summary>
+    /// 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).
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
     [Benchmark]
     public ulong InlineUInt128WithLocalsAndOperandReduction()
     {
         return MulMod64InlineWithLocalsAndOperandReduction(Input.Left, Input.Right, Input.Modulus);
     }
 
+    /// <summary>
+    /// 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.
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
     [Benchmark]
     public ulong MultiplyHighDecomposition()
     {
         return MulMod64MulHigh(Input.Left, Input.Right, Input.Modulus);
     }
 
+    /// <summary>
+    /// Uses the GPU-friendly <see cref="ULongExtensions.MulMod64GpuCompatible(ulong, ulong, ulong)"/> 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.
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
+    [Benchmark]
+    public ulong GpuCompatibleBaseline()
+    {
+        return Input.Left.MulMod64GpuCompatible(Input.Right, Input.Modulus);
+    }
+
+    /// <summary>
+    /// Uses the native-modulo GPU helper <see cref="ULongExtensions.MulMod64GpuCompatibleDeferred(ulong, ulong, ulong)"/>;
+    /// 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).
+    /// </summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
+    [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<GpuUInt128>, IEquatable<GpuUInt128>
 
     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<GpuUInt128>, IEquatable<GpuUInt128>
         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<GpuUInt128>, IEquatable<GpuUInt128>
         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<GpuUInt128>, IEquatable<GpuUInt128>
     // - 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)
     {