Merge pull request #229 from mlemanczyk/codex/analyze-and-implement-exponentremainderstepper-for-gpu-lbfhg2

Optimize GPU divisor heuristics pipeline

Benchmarks.txt

@@ -23,7 +23,7 @@ Job=.NET 8.0  Runtime=.NET 8.0
 |------------------------------------------- |------------------ |-----------:|----------:|----------:|------:|--------:|----------:|------------:|
 | ExtensionBaseline                          | CrossWordBlend    |   3.606 ns | 0.0030 ns | 0.0025 ns |  1.00 |    0.00 |         - |          NA |
 | GpuCompatibleMulModExtension               | CrossWordBlend    |  22.870 ns | 0.1463 ns | 0.1297 ns |  6.34 |    0.04 |         - |          NA |
-| GpuCompatibleMulModSimplifiedExtension     | CrossWordBlend    |  22.780 ns | 0.4463 ns | 0.8600 ns |  6.32 |    0.23 |         - |          NA |
+| MulModGpu     | CrossWordBlend    |  22.780 ns | 0.4463 ns | 0.8600 ns |  6.32 |    0.23 |         - |          NA |
 | InlineUInt128Operands                      | CrossWordBlend    |  22.218 ns | 0.0319 ns | 0.0266 ns |  6.16 |    0.01 |         - |          NA |
 | InlineUInt128OperandsWithReductionFirst    | CrossWordBlend    |  22.626 ns | 0.0553 ns | 0.0517 ns |  6.27 |    0.02 |         - |          NA |
 | InlineUInt128OperandsWithOperandReduction  | CrossWordBlend    |  22.548 ns | 0.0487 ns | 0.0432 ns |  6.25 |    0.01 |         - |          NA |
@@ -35,7 +35,7 @@ Job=.NET 8.0  Runtime=.NET 8.0
 |                                            |                   |            |           |           |       |         |           |             |
 | ExtensionBaseline                          | MixedBitPattern   |   3.564 ns | 0.0059 ns | 0.0052 ns |  1.00 |    0.00 |         - |          NA |
 | GpuCompatibleMulModExtension               | MixedBitPattern   |  64.800 ns | 0.3399 ns | 0.3179 ns | 18.17 |    0.09 |         - |          NA |
-| GpuCompatibleMulModSimplifiedExtension     | MixedBitPattern   |  59.165 ns | 0.4243 ns | 0.3969 ns | 16.60 |    0.11 |         - |          NA |
+| MulModGpu     | MixedBitPattern   |  59.165 ns | 0.4243 ns | 0.3969 ns | 16.60 |    0.11 |         - |          NA |
 | InlineUInt128Operands                      | MixedBitPattern   |  30.539 ns | 0.1396 ns | 0.1306 ns |  8.57 |    0.03 |         - |          NA |
 | InlineUInt128OperandsWithReductionFirst    | MixedBitPattern   |  31.237 ns | 0.1013 ns | 0.0898 ns |  8.76 |    0.03 |         - |          NA |
 | InlineUInt128OperandsWithOperandReduction  | MixedBitPattern   |  30.210 ns | 0.1323 ns | 0.1237 ns |  8.48 |    0.04 |         - |          NA |
@@ -47,7 +47,7 @@ Job=.NET 8.0  Runtime=.NET 8.0
 |                                            |                   |            |           |           |       |         |           |             |
 | ExtensionBaseline                          | NearFullRange     |   3.583 ns | 0.0105 ns | 0.0098 ns |  1.00 |    0.00 |         - |          NA |
 | GpuCompatibleMulModExtension               | NearFullRange     |   5.708 ns | 0.0520 ns | 0.0486 ns |  1.59 |    0.02 |         - |          NA |
-| GpuCompatibleMulModSimplifiedExtension     | NearFullRange     |   5.952 ns | 0.1253 ns | 0.1172 ns |  1.66 |    0.04 |         - |          NA |
+| MulModGpu     | NearFullRange     |   5.952 ns | 0.1253 ns | 0.1172 ns |  1.66 |    0.04 |         - |          NA |
 | InlineUInt128Operands                      | NearFullRange     |  27.159 ns | 0.0936 ns | 0.0781 ns |  7.58 |    0.03 |         - |          NA |
 | InlineUInt128OperandsWithReductionFirst    | NearFullRange     |   4.357 ns | 0.0205 ns | 0.0182 ns |  1.22 |    0.00 |         - |          NA |
 | InlineUInt128OperandsWithOperandReduction  | NearFullRange     |   4.319 ns | 0.0158 ns | 0.0140 ns |  1.21 |    0.00 |         - |          NA |
@@ -59,7 +59,7 @@ Job=.NET 8.0  Runtime=.NET 8.0
 |                                            |                   |            |           |           |       |         |           |             |
 | ExtensionBaseline                          | PrimeSizedModulus |   3.593 ns | 0.0039 ns | 0.0037 ns |  1.00 |    0.00 |         - |          NA |
 | GpuCompatibleMulModExtension               | PrimeSizedModulus |  52.766 ns | 0.3086 ns | 0.2887 ns | 14.69 |    0.08 |         - |          NA |
-| GpuCompatibleMulModSimplifiedExtension     | PrimeSizedModulus |  50.490 ns | 0.1977 ns | 0.1752 ns | 14.05 |    0.05 |         - |          NA |
+| MulModGpu     | PrimeSizedModulus |  50.490 ns | 0.1977 ns | 0.1752 ns | 14.05 |    0.05 |         - |          NA |
 | InlineUInt128Operands                      | PrimeSizedModulus |  31.409 ns | 0.0984 ns | 0.0920 ns |  8.74 |    0.03 |         - |          NA |
 | InlineUInt128OperandsWithReductionFirst    | PrimeSizedModulus |  32.176 ns | 0.3427 ns | 0.3206 ns |  8.96 |    0.09 |         - |          NA |
 | InlineUInt128OperandsWithOperandReduction  | PrimeSizedModulus |  31.786 ns | 0.1483 ns | 0.1387 ns |  8.85 |    0.04 |         - |          NA |
@@ -71,7 +71,7 @@ Job=.NET 8.0  Runtime=.NET 8.0
 |                                            |                   |            |           |           |       |         |           |             |
 | ExtensionBaseline                          | SparseOperands    |   3.555 ns | 0.0087 ns | 0.0082 ns |  1.00 |    0.00 |         - |          NA |
 | GpuCompatibleMulModExtension               | SparseOperands    |   5.671 ns | 0.0043 ns | 0.0036 ns |  1.60 |    0.00 |         - |          NA |
-| GpuCompatibleMulModSimplifiedExtension     | SparseOperands    |   5.239 ns | 0.0103 ns | 0.0096 ns |  1.47 |    0.00 |         - |          NA |
+| MulModGpu     | SparseOperands    |   5.239 ns | 0.0103 ns | 0.0096 ns |  1.47 |    0.00 |         - |          NA |
 | InlineUInt128Operands                      | SparseOperands    |   3.563 ns | 0.0217 ns | 0.0203 ns |  1.00 |    0.01 |         - |          NA |
 | InlineUInt128OperandsWithReductionFirst    | SparseOperands    |   4.617 ns | 0.0095 ns | 0.0089 ns |  1.30 |    0.00 |         - |          NA |
 | InlineUInt128OperandsWithOperandReduction  | SparseOperands    |   4.282 ns | 0.0193 ns | 0.0181 ns |  1.20 |    0.01 |         - |          NA |
@@ -83,7 +83,7 @@ Job=.NET 8.0  Runtime=.NET 8.0
 |                                            |                   |            |           |           |       |         |           |             |
 | ExtensionBaseline                          | ZeroOperands      |   3.599 ns | 0.0069 ns | 0.0065 ns |  1.00 |    0.00 |         - |          NA |
 | GpuCompatibleMulModExtension               | ZeroOperands      |   1.976 ns | 0.0053 ns | 0.0049 ns |  0.55 |    0.00 |         - |          NA |
-| GpuCompatibleMulModSimplifiedExtension     | ZeroOperands      |   1.986 ns | 0.0048 ns | 0.0045 ns |  0.55 |    0.00 |         - |          NA |
+| MulModGpu     | ZeroOperands      |   1.986 ns | 0.0048 ns | 0.0045 ns |  0.55 |    0.00 |         - |          NA |
 | InlineUInt128Operands                      | ZeroOperands      |   3.542 ns | 0.0125 ns | 0.0117 ns |  0.98 |    0.00 |         - |          NA |
 | InlineUInt128OperandsWithReductionFirst    | ZeroOperands      |   4.286 ns | 0.0178 ns | 0.0167 ns |  1.19 |    0.01 |         - |          NA |
 | InlineUInt128OperandsWithOperandReduction  | ZeroOperands      |   4.256 ns | 0.0116 ns | 0.0097 ns |  1.18 |    0.00 |         - |          NA |
@@ -97,7 +97,7 @@ Job=.NET 8.0  Runtime=.NET 8.0
 Outliers
   MulMod64Benchmarks.ExtensionBaseline: .NET 8.0                          -> 2 outliers were removed (4.92 ns, 5.01 ns)
   MulMod64Benchmarks.GpuCompatibleMulModExtension: .NET 8.0               -> 1 outlier  was  removed, 2 outliers were detected (23.75 ns, 25.37 ns)
-  MulMod64Benchmarks.GpuCompatibleMulModSimplifiedExtension: .NET 8.0     -> 2 outliers were removed, 3 outliers were detected (22.07 ns, 26.53 ns, 26.65 ns)
+  MulMod64Benchmarks.MulModGpu: .NET 8.0     -> 2 outliers were removed, 3 outliers were detected (22.07 ns, 26.53 ns, 26.65 ns)
   MulMod64Benchmarks.InlineUInt128Operands: .NET 8.0                      -> 2 outliers were removed (23.62 ns, 23.64 ns)
   MulMod64Benchmarks.InlineUInt128OperandsWithOperandReduction: .NET 8.0  -> 1 outlier  was  removed (24.42 ns)
   MulMod64Benchmarks.MultiplyHighDecomposition: .NET 8.0                  -> 2 outliers were removed (25.12 ns, 25.31 ns)
@@ -112,7 +112,7 @@ Outliers
   MulMod64Benchmarks.InlineUInt128Operands: .NET 8.0                      -> 2 outliers were removed (29.01 ns, 29.49 ns)
   MulMod64Benchmarks.InlineUInt128OperandsWithReductionFirst: .NET 8.0    -> 1 outlier  was  removed (5.72 ns)
   MulMod64Benchmarks.InlineUInt128OperandsWithOperandReduction: .NET 8.0  -> 1 outlier  was  removed (5.73 ns)
-  MulMod64Benchmarks.GpuCompatibleMulModSimplifiedExtension: .NET 8.0     -> 1 outlier  was  removed, 3 outliers were detected (51.41 ns, 51.47 ns, 52.10 ns)
+  MulMod64Benchmarks.MulModGpu: .NET 8.0     -> 1 outlier  was  removed, 3 outliers were detected (51.41 ns, 51.47 ns, 52.10 ns)
   MulMod64Benchmarks.InlineUInt128Operands: .NET 8.0                      -> 1 outlier  was  detected (32.44 ns)
   MulMod64Benchmarks.InlineUInt128WithLocals: .NET 8.0                    -> 2 outliers were removed (33.47 ns, 34.46 ns)
   MulMod64Benchmarks.MultiplyHighDecomposition: .NET 8.0                  -> 1 outlier  was  removed (33.93 ns)

EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.Gpu.cs

@@ -0,0 +1,461 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Jobs;
+using ILGPU;
+using ILGPU.Runtime;
+using PerfectNumbers.Core;
+using PerfectNumbers.Core.Gpu;
+
+namespace EvenPerfectBitScanner.Benchmarks;
+
+/// <summary>
+/// GPU benchmarks for 64-bit modular multiplication helpers.
+/// Windows 11 results ranged from 0.62 ms to 1.63 ms, and MulModGpu led the PrimeSizedModulus case at 0.641 ms.
+/// Recommended helper for p ≥ 138,000,000 is MulModGpu.
+/// </summary>
+[SimpleJob(RuntimeMoniker.Net80)]
+[MemoryDiagnoser]
+public class MulMod64BenchmarksGpu : IDisposable
+{
+    private static readonly MulMod64Input[] Inputs =
+    [
+        new(0UL, 0UL, 3UL, "ZeroOperands"),
+        new(ulong.MaxValue, ulong.MaxValue - 1UL, ulong.MaxValue - 58UL, "NearFullRange"),
+        new(0xFFFF_FFFF_0000_0000UL, 0x0000_0000_FFFF_FFFFUL, 0xFFFF_FFFF_0000_002BUL, "CrossWordBlend"),
+        new(0x1234_5678_9ABC_DEF0UL, 0x0FED_CBA9_8765_4321UL, 0x1FFF_FFFF_FFFF_FFFBUL, "MixedBitPattern"),
+        new(0x0000_0001_0000_0000UL, 0x0000_0000_0000_0003UL, 0x0000_0000_FFFF_FFC3UL, "SparseOperands"),
+        new(0x7FFF_FFFF_FFFF_FFA3UL, 0x6FFF_FFFF_FFFF_FF81UL, 0x7FFF_FFFF_FFFF_FFE7UL, "PrimeSizedModulus"),
+    ];
+
+    [ParamsSource(nameof(GetInputs))]
+    public MulMod64Input Input { get; set; }
+
+    private Context? _context;
+    private Accelerator? _accelerator;
+    private MemoryBuffer1D<MulMod64GpuKernelInput, Stride1D.Dense>? _inputBuffer;
+    private MemoryBuffer1D<ulong, Stride1D.Dense>? _resultBuffer;
+
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _extensionBaselineKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _gpuCompatibleExtensionKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _gpuCompatibleSimplifiedKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _inlineKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _inlineReductionFirstKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _inlineOperandReductionKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _inlineLocalsKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _inlineLocalsOperandReductionKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _mulHighKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _legacyKernel;
+    private Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? _legacyDeferredKernel;
+
+    public static IEnumerable<MulMod64Input> GetInputs()
+    {
+        return Inputs;
+    }
+
+    [GlobalSetup]
+    public void GlobalSetup()
+    {
+        _context = Context.Create(builder => builder.Default().EnableAlgorithms());
+        _accelerator = _context.GetPreferredDevice(preferCPU: false).CreateAccelerator(_context);
+
+        MulMod64GpuKernelInput[] gpuInputs = Inputs
+            .Select(i => new MulMod64GpuKernelInput(i.Left, i.Right, i.Modulus))
+            .ToArray();
+
+        _inputBuffer = _accelerator.Allocate1D(gpuInputs);
+        _resultBuffer = _accelerator.Allocate1D<ulong>(gpuInputs.Length);
+
+        _extensionBaselineKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(ExtensionBaselineKernel);
+        _gpuCompatibleExtensionKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(GpuCompatibleExtensionKernel);
+        _gpuCompatibleSimplifiedKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(GpuCompatibleSimplifiedKernel);
+        _inlineKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(InlineKernel);
+        _inlineReductionFirstKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(InlineReductionFirstKernel);
+        _inlineOperandReductionKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(InlineOperandReductionKernel);
+        _inlineLocalsKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(InlineLocalsKernel);
+        _inlineLocalsOperandReductionKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(InlineLocalsOperandReductionKernel);
+        _mulHighKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(MulHighKernel);
+        _legacyKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(LegacyKernel);
+        _legacyDeferredKernel = _accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>(LegacyDeferredKernel);
+    }
+
+    /// <summary>
+    /// Inline wide multiply/reduce baseline: ZeroOperands 0.618 ms, MixedBitPattern 0.640 ms, CrossWordBlend 0.928 ms, NearFullRange 0.875 ms, PrimeSizedModulus 0.759 ms, SparseOperands 0.804 ms.
+    /// Use primarily as a reference; optimized GPU-specific helpers beat it across all datasets.
+    /// </summary>
+    [Benchmark(Baseline = true)]
+    public void ExtensionBaseline()
+    {
+        LaunchKernel(_extensionBaselineKernel);
+    }
+
+    /// <summary>
+    /// Original GPU-compatible extension path: means from 0.614 ms (MixedBitPattern) to 1.023 ms (CrossWordBlend) with PrimeSizedModulus at 0.697 ms.
+    /// Retain only when mirroring the CPU extension logic outweighs the speed gap to the simplified helper.
+    /// </summary>
+    [Benchmark]
+    public void GpuCompatibleMulModExtension()
+    {
+        LaunchKernel(_gpuCompatibleExtensionKernel);
+    }
+
+    /// <summary>
+    /// Simplified GPU extension helper: 0.708 ms (ZeroOperands), 0.905 ms (CrossWordBlend), 0.798 ms (NearFullRange), and class-leading 0.641 ms (PrimeSizedModulus).
+    /// Recommended for production GPU scans, especially for large p ≥ 138,000,000.
+    /// </summary>
+    [Benchmark]
+    public void MulModGpu()
+    {
+        LaunchKernel(_gpuCompatibleSimplifiedKernel);
+    }
+
+    /// <summary>
+    /// Inline UInt128 multiplication without operand reduction: CrossWordBlend 0.838 ms, MixedBitPattern 0.865 ms, NearFullRange 0.786 ms, PrimeSizedModulus 0.709 ms, SparseOperands 0.751 ms, ZeroOperands 0.827 ms.
+    /// Suitable when explicit UInt128 math is acceptable, though it trails the simplified GPU helper on every dataset.
+    /// </summary>
+    [Benchmark]
+    public void InlineUInt128Operands()
+    {
+        LaunchKernel(_inlineKernel);
+    }
+
+    /// <summary>
+    /// Reduces operands before multiplying: ZeroOperands 0.926 ms, CrossWordBlend 0.847 ms, MixedBitPattern 0.843 ms, NearFullRange 0.689 ms, PrimeSizedModulus 0.809 ms, SparseOperands 0.746 ms.
+    /// Only competitive on NearFullRange; the extra reductions hurt other scenarios.
+    /// </summary>
+    [Benchmark]
+    public void InlineUInt128OperandsWithReductionFirst()
+    {
+        LaunchKernel(_inlineReductionFirstKernel);
+    }
+
+    /// <summary>
+    /// Performs operand reduction and then multiplies inline: CrossWordBlend 1.041 ms, MixedBitPattern 0.913 ms, NearFullRange 0.847 ms, PrimeSizedModulus 0.812 ms, SparseOperands 0.744 ms, ZeroOperands 0.802 ms.
+    /// Consistently slower than both the reduction-first and simplified GPU helpers.
+    /// </summary>
+    [Benchmark]
+    public void InlineUInt128OperandsWithOperandReduction()
+    {
+        LaunchKernel(_inlineOperandReductionKernel);
+    }
+
+    /// <summary>
+    /// Uses local temporaries for the wide product: ZeroOperands 0.691 ms, CrossWordBlend 0.914 ms, MixedBitPattern 0.874 ms, NearFullRange 0.702 ms, PrimeSizedModulus 0.817 ms, SparseOperands 0.728 ms.
+    /// Balanced choice when sticking with manual wide products, but still slower than the simplified GPU path.
+    /// </summary>
+    [Benchmark]
+    public void InlineUInt128WithLocals()
+    {
+        LaunchKernel(_inlineLocalsKernel);
+    }
+
+    /// <summary>
+    /// Combines local temporaries with operand reduction: ZeroOperands 0.698 ms, CrossWordBlend 0.858 ms, MixedBitPattern 0.877 ms, NearFullRange 0.851 ms, PrimeSizedModulus 0.761 ms, SparseOperands 0.779 ms.
+    /// Adds little benefit over the non-reduced local variant while increasing cost on most datasets.
+    /// </summary>
+    [Benchmark]
+    public void InlineUInt128WithLocalsAndOperandReduction()
+    {
+        LaunchKernel(_inlineLocalsOperandReductionKernel);
+    }
+
+    /// <summary>
+    /// Multiply-high decomposition: ZeroOperands 0.651 ms, CrossWordBlend 0.848 ms, MixedBitPattern 0.758 ms, NearFullRange 0.812 ms, PrimeSizedModulus 0.728 ms, SparseOperands 0.699 ms.
+    /// Useful when avoiding UInt128 temporaries on sparse inputs, but still behind the simplified GPU helper for prime-sized moduli.
+    /// </summary>
+    [Benchmark]
+    public void MultiplyHighDecomposition()
+    {
+        LaunchKernel(_mulHighKernel);
+    }
+
+    /// <summary>
+    /// Legacy GPU baseline: ZeroOperands 0.759 ms, CrossWordBlend 0.971 ms, MixedBitPattern 0.827 ms, NearFullRange 0.870 ms, PrimeSizedModulus 0.737 ms, SparseOperands 0.620 ms.
+    /// Retain only for regression comparisons; newer helpers are faster except on the sparse sample.
+    /// </summary>
+    [Benchmark]
+    public void GpuCompatibleBaseline()
+    {
+        LaunchKernel(_legacyKernel);
+    }
+
+    /// <summary>
+    /// Deferred legacy helper: ZeroOperands 1.561 ms, CrossWordBlend 1.631 ms, MixedBitPattern 1.540 ms, NearFullRange 1.438 ms, PrimeSizedModulus 1.598 ms, SparseOperands 1.096 ms.
+    /// Slowest option in every scenario; keep only for historical context.
+    /// </summary>
+    [Benchmark]
+    public void GpuCompatibleDeferred()
+    {
+        LaunchKernel(_legacyDeferredKernel);
+    }
+
+    [GlobalCleanup]
+    public void GlobalCleanup()
+    {
+        _resultBuffer?.Dispose();
+        _inputBuffer?.Dispose();
+        _accelerator?.Dispose();
+        _context?.Dispose();
+    }
+
+    public void Dispose()
+    {
+        GlobalCleanup();
+        GC.SuppressFinalize(this);
+    }
+
+    private void LaunchKernel(Action<Index1D, ArrayView<MulMod64GpuKernelInput>, ArrayView<ulong>>? kernel)
+    {
+        if (kernel == null || _accelerator == null || _inputBuffer == null || _resultBuffer == null)
+        {
+            throw new InvalidOperationException("Benchmark kernels are not initialized.");
+        }
+
+        kernel(new Index1D((int)_inputBuffer.Length), _inputBuffer.View, _resultBuffer.View);
+        _accelerator.Synchronize();
+    }
+
+    private static void ExtensionBaselineKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        results[index] = MulMod64Baseline(input.Left, input.Right, input.Modulus);
+    }
+
+    private static void GpuCompatibleExtensionKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        GpuUInt128 state = new(input.Left);
+        results[index] = state.MulMod(input.Right, input.Modulus);
+    }
+
+    private static void GpuCompatibleSimplifiedKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        GpuUInt128 state = new(input.Left);
+        results[index] = state.MulMod(input.Right, input.Modulus);
+    }
+
+    private static void InlineKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        results[index] = MulMod64Baseline(input.Left, input.Right, input.Modulus);
+    }
+
+    private static void InlineReductionFirstKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        ulong left = input.Left % input.Modulus;
+        ulong right = input.Right % input.Modulus;
+        results[index] = MulMod64Baseline(left, right, input.Modulus);
+    }
+
+    private static void InlineOperandReductionKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        ulong left = input.Left % input.Modulus;
+        ulong right = input.Right % input.Modulus;
+        results[index] = MulMod64Baseline(left, right, input.Modulus);
+    }
+
+    private static void InlineLocalsKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        MultiplyWide(input.Left, input.Right, out ulong productHigh, out ulong productLow);
+        results[index] = ReduceWide(productHigh, productLow, input.Modulus);
+    }
+
+    private static void InlineLocalsOperandReductionKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        ulong left = input.Left % input.Modulus;
+        ulong right = input.Right % input.Modulus;
+        MultiplyWide(left, right, out ulong productHigh, out ulong productLow);
+        results[index] = ReduceWide(productHigh, productLow, input.Modulus);
+    }
+
+    private static void MulHighKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        ulong left = input.Left % input.Modulus;
+        ulong right = input.Right % input.Modulus;
+        MultiplyWide(left, right, out ulong productHigh, out ulong productLow);
+        results[index] = ReduceWide(productHigh, productLow, input.Modulus);
+    }
+
+    private static void LegacyKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        results[index] = MulMod64GpuLegacy(input.Left, input.Right, input.Modulus);
+    }
+
+    private static void LegacyDeferredKernel(Index1D index, ArrayView<MulMod64GpuKernelInput> inputs, ArrayView<ulong> results)
+    {
+        if (index >= inputs.Length)
+        {
+            return;
+        }
+
+        MulMod64GpuKernelInput input = inputs[index];
+        results[index] = MulMod64GpuDeferredLegacy(input.Left, input.Right, input.Modulus);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong MulMod64Baseline(ulong left, ulong right, ulong modulus)
+    {
+        MultiplyWide(left, right, out ulong high, out ulong low);
+        return ReduceWide(high, low, modulus);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void MultiplyWide(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 productLL = leftLow * rightLow;
+        ulong productLH = leftLow * rightHigh;
+        ulong productHL = leftHigh * rightLow;
+        ulong productHH = leftHigh * rightHigh;
+
+        ulong carry = (productLL >> 32) + (productLH & 0xFFFF_FFFFUL) + (productHL & 0xFFFF_FFFFUL);
+        low = (carry << 32) | (uint)productLL;
+        high = productHH + (productLH >> 32) + (productHL >> 32) + (carry >> 32);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong ReduceWide(ulong high, ulong low, ulong modulus)
+    {
+        if (modulus <= 1UL)
+        {
+            return 0UL;
+        }
+
+        ulong remainder = 0UL;
+        if (high != 0UL)
+        {
+            for (int bit = 63; bit >= 0; --bit)
+            {
+                remainder = ShiftLeftAddBit(remainder, (high >> bit) & 1UL, modulus);
+            }
+        }
+
+        for (int bit = 63; bit >= 0; --bit)
+        {
+            remainder = ShiftLeftAddBit(remainder, (low >> bit) & 1UL, modulus);
+        }
+
+        return remainder;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong ShiftLeftAddBit(ulong value, ulong bit, ulong modulus)
+    {
+        value <<= 1;
+        if (value >= modulus)
+        {
+            value -= modulus;
+        }
+
+        if (bit != 0UL)
+        {
+            value += 1UL;
+            if (value >= modulus)
+            {
+                value -= modulus;
+            }
+        }
+
+        return value;
+    }
+
+    [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);
+    }
+
+    public readonly record struct MulMod64Input(ulong Left, ulong Right, ulong Modulus, string Name)
+    {
+        public override string ToString()
+        {
+            return Name;
+        }
+    }
+
+    private readonly struct MulMod64GpuKernelInput
+    {
+        public readonly ulong Left;
+        public readonly ulong Right;
+        public readonly ulong Modulus;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public MulMod64GpuKernelInput(ulong left, ulong right, ulong modulus)
+        {
+            Left = left;
+            Right = right;
+            Modulus = modulus;
+        }
+    }
+}

EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs

@@ -58,7 +58,7 @@ public class MulMod64Benchmarks
     }
 
     /// <summary>
-    /// Calls the UInt128-based <see cref="ULongExtensions.MulModSimplified(ulong, ulong, ulong)"/> helper; performance mirrors
+    /// Calls the UInt128-based <see cref="ULongExtensions.MulMod(ulong, ulong, ulong)"/> helper; performance mirrors
     /// the full helper with 5.24–5.95 ns on sparse inputs and 1.99 ns when zeros dominate, but 22.8–59.2 ns on dense blends.
     /// </summary>
     /// <remarks>
@@ -66,10 +66,10 @@ public class MulMod64Benchmarks
     /// PrimeSizedModulus 50.490 ns, SparseOperands 5.239 ns, ZeroOperands 1.986 ns.
     /// </remarks>
     [Benchmark]
-    public ulong GpuCompatibleMulModSimplifiedExtension()
+    public ulong MulModGpu()
     {
         GpuUInt128 gpuUInt128 = new(Input.Left);
-        return gpuUInt128.MulModSimplified(Input.Right, Input.Modulus);
+        return gpuUInt128.MulMod(Input.Right, Input.Modulus);
     }
 
     /// <summary>

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs

@@ -0,0 +1,94 @@
+using System;
+using PerfectNumbers.Core;
+
+namespace EvenPerfectBitScanner.Benchmarks;
+
+internal sealed class Pow2MontgomeryModWindowedBenchmarkInputs
+{
+    public Pow2MontgomeryModWindowedBenchmarkInputs(Pow2MontgomeryModWindowedBenchmarkCase[] cases)
+    {
+        Cases = cases;
+    }
+
+    public Pow2MontgomeryModWindowedBenchmarkCase[] Cases { get; }
+}
+
+internal readonly struct Pow2MontgomeryModWindowedBenchmarkCase
+{
+    public Pow2MontgomeryModWindowedBenchmarkCase(
+        ulong exponent,
+        ulong modulus,
+        in MontgomeryDivisorData divisor,
+        ulong cycleLength,
+        ulong reducedExponent)
+    {
+        Exponent = exponent;
+        Modulus = modulus;
+        Divisor = divisor;
+        CycleLength = cycleLength;
+        ReducedExponent = reducedExponent;
+    }
+
+    public ulong Exponent { get; }
+
+    public ulong Modulus { get; }
+
+    public MontgomeryDivisorData Divisor { get; }
+
+    public ulong CycleLength { get; }
+
+    public ulong ReducedExponent { get; }
+}
+
+internal static class Pow2MontgomeryModWindowedBenchmarkInputsProvider
+{
+    private static readonly Lazy<Pow2MontgomeryModWindowedBenchmarkInputs> Cache = new(Create);
+
+    public static Pow2MontgomeryModWindowedBenchmarkInputs Instance => Cache.Value;
+
+    private static Pow2MontgomeryModWindowedBenchmarkInputs Create()
+    {
+        Pow2MontgomeryModWindowedBenchmarkSeed[] seeds =
+        {
+            new Pow2MontgomeryModWindowedBenchmarkSeed(138_000_001UL, 1UL),
+            new Pow2MontgomeryModWindowedBenchmarkSeed(150_000_013UL, 2UL),
+            new Pow2MontgomeryModWindowedBenchmarkSeed(175_000_019UL, 3UL),
+            new Pow2MontgomeryModWindowedBenchmarkSeed(190_000_151UL, 4UL),
+            new Pow2MontgomeryModWindowedBenchmarkSeed(210_000_089UL, 5UL),
+            new Pow2MontgomeryModWindowedBenchmarkSeed(230_000_039UL, 6UL),
+            new Pow2MontgomeryModWindowedBenchmarkSeed(260_000_111UL, 7UL),
+            new Pow2MontgomeryModWindowedBenchmarkSeed(300_000_007UL, 8UL)
+        };
+
+        var cases = new Pow2MontgomeryModWindowedBenchmarkCase[seeds.Length];
+        for (int i = 0; i < seeds.Length; i++)
+        {
+            Pow2MontgomeryModWindowedBenchmarkSeed seed = seeds[i];
+            ulong modulus = 2UL * seed.Multiplier * seed.Exponent + 1UL;
+            MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus);
+            ulong cycleLength = MersenneDivisorCycles.CalculateCycleLength(modulus, divisor);
+            ulong reducedExponent = seed.Exponent % cycleLength;
+            cases[i] = new Pow2MontgomeryModWindowedBenchmarkCase(
+                seed.Exponent,
+                modulus,
+                divisor,
+                cycleLength,
+                reducedExponent);
+        }
+
+        return new Pow2MontgomeryModWindowedBenchmarkInputs(cases);
+    }
+
+    private readonly struct Pow2MontgomeryModWindowedBenchmarkSeed
+    {
+        public Pow2MontgomeryModWindowedBenchmarkSeed(ulong exponent, ulong multiplier)
+        {
+            Exponent = exponent;
+            Multiplier = multiplier;
+        }
+
+        public ulong Exponent { get; }
+
+        public ulong Multiplier { get; }
+    }
+}

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks.cs

@@ -0,0 +1,94 @@
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Jobs;
+using PerfectNumbers.Core;
+
+namespace EvenPerfectBitScanner.Benchmarks;
+
+/// <summary>
+/// Benchmarks CPU Montgomery exponentiation helpers for exponents at or above 138 million.
+/// Reuses the shared benchmark input cache so every method processes identical samples.
+/// </summary>
+[MemoryDiagnoser]
+[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
+public class Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks
+{
+    private Pow2MontgomeryModWindowedBenchmarkInputs _inputs = null!;
+
+    [GlobalSetup]
+    public void Setup()
+    {
+        _inputs = Pow2MontgomeryModWindowedBenchmarkInputsProvider.Instance;
+    }
+
+    /// <summary>
+    /// Baseline CPU windowed Montgomery exponentiation that converts the residue back to standard form.
+    /// </summary>
+    [Benchmark(Baseline = true)]
+    public ulong WindowedMontgomeryConvert()
+    {
+        ulong checksum = 0UL;
+        Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases;
+
+        for (int i = 0; i < cases.Length; i++)
+        {
+            ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i];
+            checksum ^= current.Exponent.Pow2MontgomeryModWindowedCpu(current.Divisor, keepMontgomery: false);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// CPU windowed Montgomery exponentiation that keeps the residue in the Montgomery domain.
+    /// </summary>
+    [Benchmark]
+    public ulong WindowedMontgomeryKeep()
+    {
+        ulong checksum = 0UL;
+        Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases;
+
+        for (int i = 0; i < cases.Length; i++)
+        {
+            ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i];
+            checksum ^= current.Exponent.Pow2MontgomeryModWindowedCpu(current.Divisor, keepMontgomery: true);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// Uses the known cycle length to fold the exponent before running the windowed Montgomery ladder.
+    /// </summary>
+    [Benchmark]
+    public ulong WindowedMontgomeryWithCycle()
+    {
+        ulong checksum = 0UL;
+        Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases;
+
+        for (int i = 0; i < cases.Length; i++)
+        {
+            ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i];
+            checksum ^= current.Exponent.Pow2MontgomeryModWithCycleCpu(current.CycleLength, current.Divisor);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// Runs the windowed Montgomery ladder on the reduced exponent that the divisor cycle exposes.
+    /// </summary>
+    [Benchmark]
+    public ulong WindowedMontgomeryFromCycleRemainder()
+    {
+        ulong checksum = 0UL;
+        Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases;
+
+        for (int i = 0; i < cases.Length; i++)
+        {
+            ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i];
+            checksum ^= current.ReducedExponent.Pow2MontgomeryModFromCycleRemainderCpu(current.Divisor);
+        }
+
+        return checksum;
+    }
+}

EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs

@@ -0,0 +1,78 @@
+using System;
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Jobs;
+using PerfectNumbers.Core;
+
+namespace EvenPerfectBitScanner.Benchmarks;
+
+/// <summary>
+/// Benchmarks GPU Montgomery exponentiation helpers for exponents at or above 138 million.
+/// The sample set is generated once and reused across all benchmark methods so every variant
+/// processes the same data without incurring repeated initialization costs.
+/// </summary>
+[MemoryDiagnoser]
+[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)]
+public class Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks
+{
+    private Pow2MontgomeryModWindowedBenchmarkInputs _inputs = null!;
+
+    [GlobalSetup]
+    public void Setup()
+    {
+        _inputs = Pow2MontgomeryModWindowedBenchmarkInputsProvider.Instance;
+    }
+
+    /// <summary>
+    /// Baseline GPU windowed Montgomery exponentiation that converts back to the standard residue.
+    /// </summary>
+    [Benchmark(Baseline = true)]
+    public ulong WindowedMontgomeryConvert()
+    {
+        ulong checksum = 0UL;
+        Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases;
+
+        for (int i = 0; i < cases.Length; i++)
+        {
+            ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i];
+            checksum ^= current.Exponent.Pow2MontgomeryModWindowedGpu(current.Divisor, keepMontgomery: false);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// GPU windowed Montgomery exponentiation that keeps the result in the Montgomery domain.
+    /// </summary>
+    [Benchmark]
+    public ulong WindowedMontgomeryKeep()
+    {
+        ulong checksum = 0UL;
+        Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases;
+
+        for (int i = 0; i < cases.Length; i++)
+        {
+            ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i];
+            checksum ^= current.Exponent.Pow2MontgomeryModWindowedGpu(current.Divisor, keepMontgomery: true);
+        }
+
+        return checksum;
+    }
+
+    /// <summary>
+    /// GPU windowed modular exponentiation without Montgomery reduction for comparison.
+    /// </summary>
+    [Benchmark]
+    public ulong WindowedStandardMod()
+    {
+        ulong checksum = 0UL;
+        Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases;
+
+        for (int i = 0; i < cases.Length; i++)
+        {
+            ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i];
+            checksum ^= current.Exponent.Pow2ModWindowedGpu(current.Modulus);
+        }
+
+        return checksum;
+    }
+}

PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs

@@ -74,18 +74,6 @@ public class GpuKernelCompilationTests
         typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetKernel", BindingFlags.NonPublic | BindingFlags.Instance)
         ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetKernel not found.");
 
-    private static readonly MethodInfo DivisorByDivisorGetExponentKernelMethod =
-        typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetKernelByPrimeExponent", BindingFlags.NonPublic | BindingFlags.Instance)
-        ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetKernelByPrimeExponent not found.");
-
-    private static readonly MethodInfo DivisorByDivisorGetRemainderDeltaKernelMethod =
-        typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetRemainderDeltaKernel", BindingFlags.NonPublic | BindingFlags.Instance)
-        ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetRemainderDeltaKernel not found.");
-
-    private static readonly MethodInfo DivisorByDivisorGetRemainderScanKernelMethod =
-        typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetRemainderScanKernel", BindingFlags.NonPublic | BindingFlags.Instance)
-        ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetRemainderScanKernel not found.");
-
     private static readonly MethodInfo DivisorTesterGetKernelMethod =
         typeof(MersenneNumberDivisorGpuTester).GetMethod("GetKernel", BindingFlags.NonPublic | BindingFlags.Instance)
         ?? throw new InvalidOperationException("MersenneNumberDivisorGpuTester.GetKernel not found.");
@@ -157,9 +145,6 @@ public class GpuKernelCompilationTests
         yield return Loader("NttGpuMath.ForwardKernel", CompileNttForwardKernel);
         yield return Loader("NttGpuMath.InverseKernel", CompileNttInverseKernel);
         yield return Loader("MersenneNumberDivisorByDivisorGpuTester.CheckKernel", CompileDivisorByDivisorCheckKernel);
-        yield return Loader("MersenneNumberDivisorByDivisorGpuTester.ComputeMontgomeryExponentKernel", CompileDivisorByDivisorExponentKernel);
-        yield return Loader("MersenneNumberDivisorByDivisorGpuTester.ComputeRemainderDeltasKernel", CompileDivisorByDivisorRemainderDeltaKernel);
-        yield return Loader("MersenneNumberDivisorByDivisorGpuTester.AccumulateRemaindersKernel", CompileDivisorByDivisorRemainderScanKernel);
         yield return Loader("MersenneNumberDivisorGpuTester.Kernel", CompileDivisorTesterKernel);
         yield return Loader("MersenneNumberLucasLehmerGpuTester.Kernel", CompileLucasLehmerKernel);
         yield return Loader("MersenneNumberLucasLehmerGpuTester.AddSmallKernel", CompileLucasLehmerAddSmallKernel);
@@ -281,24 +266,6 @@ public class GpuKernelCompilationTests
         _ = DivisorByDivisorGetKernelMethod.Invoke(tester, new object[] { accelerator });
     }
 
-    private static void CompileDivisorByDivisorExponentKernel(Accelerator accelerator)
-    {
-        var tester = new MersenneNumberDivisorByDivisorGpuTester();
-        _ = DivisorByDivisorGetExponentKernelMethod.Invoke(tester, new object[] { accelerator });
-    }
-
-    private static void CompileDivisorByDivisorRemainderDeltaKernel(Accelerator accelerator)
-    {
-        var tester = new MersenneNumberDivisorByDivisorGpuTester();
-        _ = DivisorByDivisorGetRemainderDeltaKernelMethod.Invoke(tester, new object[] { accelerator });
-    }
-
-    private static void CompileDivisorByDivisorRemainderScanKernel(Accelerator accelerator)
-    {
-        var tester = new MersenneNumberDivisorByDivisorGpuTester();
-        _ = DivisorByDivisorGetRemainderScanKernelMethod.Invoke(tester, new object[] { accelerator });
-    }
-
     private static void CompileDivisorTesterKernel(Accelerator accelerator)
     {
         var tester = new MersenneNumberDivisorGpuTester();

PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs

@@ -85,11 +85,17 @@ public class MersenneNumberDivisorGpuTesterTests
 
     [Fact]
     [Trait("Category", "Fast")]
-    public void ByDivisor_tester_requires_configuration()
+    public void ByDivisor_tester_defaults_to_zero_limit_before_configuration()
     {
         var tester = new MersenneNumberDivisorByDivisorGpuTester();
-        Action act = () => tester.IsPrime(31UL, out _);
+
-        act.Should().Throw<InvalidOperationException>();
+        tester.DivisorLimit.Should().Be(0UL);
+
+        ulong[] primes = { 31UL, 37UL };
+        ulong[] allowed = new ulong[primes.Length];
+
+        tester.PrepareCandidates(primes, allowed);
+        allowed.Should().OnlyContain(value => value == 0UL);
     }
 
     [Fact]
@@ -136,6 +142,22 @@ public class MersenneNumberDivisorGpuTesterTests
         divisorsExhausted.Should().BeTrue();
     }
 
+    [Fact]
+    [Trait("Category", "Fast")]
+    public void ByDivisor_session_computes_cycle_when_zero_is_provided()
+    {
+        var tester = new MersenneNumberDivisorByDivisorGpuTester();
+        tester.ConfigureFromMaxPrime(43UL);
+
+        using var session = tester.CreateDivisorSession();
+        ulong[] primes = { 31UL, 37UL, 41UL };
+        byte[] hits = new byte[primes.Length];
+
+        session.CheckDivisor(223UL, MontgomeryDivisorData.FromModulus(223UL), 0UL, primes, hits);
+
+        hits.Should().ContainInOrder(new byte[] { 0, 1, 0 });
+    }
+
     [Fact]
     [Trait("Category", "Fast")]
     public void ByDivisor_session_checks_divisors_across_primes()

PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs

@@ -42,7 +42,7 @@ internal sealed class MersenneCpuDivisorScanSession : IMersenneNumberDivisorByDi
 
         // Keep these remainder steppers in place so future updates continue reusing the previously computed residues.
         // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors.
-        var exponentStepper = new ExponentRemainderStepper(cachedData);
+        var exponentStepper = new ExponentRemainderStepperCpu(cachedData);
 
         var cycleStepper = new CycleRemainderStepper(divisorCycle);
 

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -175,17 +175,23 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         // Keep the divisibility filters aligned with the divisor-cycle generator so the
         // CPU path never requests cycles that were skipped during cache creation.
         bool lastIsSeven = (prime & 3UL) == 3UL;
+        const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9);
+        const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9);
 
         while (divisor.CompareTo(limit) <= 0)
         {
             ulong candidate = divisor.Low;
             processedCount++;
 
-            bool admissible = lastIsSeven
+            ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise;
-                ? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9)
+            bool admissible = ((decimalMask >> remainder10) & 1) != 0;
-                : (remainder10 == 1 || remainder10 == 3 || remainder10 == 9);
 
-            if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0)
+            if (admissible
+                && (remainder8 == 1 || remainder8 == 7)
+                && remainder3 != 0
+                && remainder5 != 0
+                && remainder7 != 0
+                && remainder11 != 0)
             {
                 MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate);
                 ulong divisorCycle;
@@ -196,7 +202,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 // }
                 // else
                 {
-                    if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(
+                    if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(
                             candidate,
                             prime,
                             divisorData,

PerfectNumbers.Core/DivisorCycleCache.cs

@@ -79,7 +79,7 @@ public sealed class DivisorCycleCache
                 throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
             }
 
-            ObserveCycleCacheDivisor(divisor);
+            // ObserveCycleCacheDivisor(divisor);
             return cached;
         }
 
@@ -114,7 +114,7 @@ public sealed class DivisorCycleCache
                     throw new InvalidDataException($"Divisor cycle is missing for {divisor}");
                 }
 
-                ObserveCycleCacheDivisor(divisor);
+                // ObserveCycleCacheDivisor(divisor);
                 cycles[i] = cached;
             }
             else

PerfectNumbers.Core/ExponentRemainderStepperCpu.cs

@@ -1,30 +1,28 @@
 namespace PerfectNumbers.Core;
 
 /// <summary>
-/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates.
+/// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path.
 /// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order.
 /// </summary>
-internal struct ExponentRemainderStepper
+internal struct ExponentRemainderStepperCpu
 {
     private readonly MontgomeryDivisorData _divisor;
     private readonly ulong _modulus;
     private readonly ulong _nPrime;
     private readonly ulong _montgomeryOne;
-    private ulong _previousExponent;
+    public ulong PreviousExponent;
     private ulong _currentMontgomery;
 
-    public ExponentRemainderStepper(in MontgomeryDivisorData divisor)
+    public ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor)
     {
         _divisor = divisor;
         _modulus = divisor.Modulus;
         _nPrime = divisor.NPrime;
         _montgomeryOne = divisor.MontgomeryOne;
-        _previousExponent = 0UL;
+        PreviousExponent = 0UL;
         _currentMontgomery = divisor.MontgomeryOne;
     }
 
-    public ulong PreviousExponent => _previousExponent;
-
     public bool MatchesDivisor(in MontgomeryDivisorData divisor)
     {
         return _modulus == divisor.Modulus && _nPrime == divisor.NPrime && _montgomeryOne == divisor.MontgomeryOne;
@@ -36,7 +34,7 @@ internal struct ExponentRemainderStepper
     /// </summary>
     public void Reset()
     {
-        _previousExponent = 0UL;
+        PreviousExponent = 0UL;
         _currentMontgomery = _montgomeryOne;
     }
 
@@ -71,38 +69,16 @@ internal struct ExponentRemainderStepper
     /// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
     public ulong ComputeNextCpu(ulong exponent)
     {
-        ulong delta = exponent - _previousExponent;
+        ulong delta = exponent - PreviousExponent;
         // TODO: Once divisor cycle lengths are mandatory, pull the delta multiplier from the
         // single-block divisor-cycle snapshot so we can skip the powmod entirely and reuse the
         // cached Montgomery residue ladder highlighted in MersenneDivisorCycleLengthGpuBenchmarks.
         ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
-        _previousExponent = exponent;
+        PreviousExponent = exponent;
         return ReduceCurrent();
     }
 
-    /// <summary>
-    /// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence.
-    /// </summary>
-    /// <param name="exponent">First exponent evaluated on the GPU path. The same divisor metadata must be reused for future calls.</param>
-    /// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
-    public ulong InitializeGpu(ulong exponent)
-    {
-        InitializeGpuState(exponent);
-        return ReduceCurrent();
-    }
-
-    /// <summary>
-    /// Initializes the GPU stepping state and reports whether the seed exponent produces unity.
-    /// </summary>
-    /// <param name="exponent">First exponent evaluated on the GPU path. Must match the divisor supplied during construction.</param>
-    /// <returns><see langword="true"/> when the seed exponent produces a Montgomery residue equal to one.</returns>
-    public bool InitializeGpuIsUnity(ulong exponent)
-    {
-        InitializeGpuState(exponent);
-        return _currentMontgomery == _montgomeryOne;
-    }
-
     /// <summary>
     /// Advances the CPU stepping state to the next exponent and reports whether it evaluates to unity.
     /// Callers must seed the stepper through <see cref="InitializeCpu(ulong)"/> or <see cref="InitializeCpuIsUnity(ulong)"/> and then
@@ -112,12 +88,12 @@ internal struct ExponentRemainderStepper
     /// <returns><see langword="true"/> when the canonical residue equals one.</returns>
     public bool ComputeNextIsUnity(ulong exponent)
     {
-        ulong delta = exponent - _previousExponent;
+        ulong delta = exponent - PreviousExponent;
         // 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.
         ulong multiplier = delta.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
         _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime);
-        _previousExponent = exponent;
+        PreviousExponent = exponent;
         return _currentMontgomery == _montgomeryOne;
     }
 
@@ -132,7 +108,7 @@ internal struct ExponentRemainderStepper
     public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity)
     {
         _currentMontgomery = montgomeryResult;
-        _previousExponent = exponent;
+        PreviousExponent = exponent;
         isUnity = _currentMontgomery == _montgomeryOne;
         return true;
     }
@@ -151,7 +127,7 @@ internal struct ExponentRemainderStepper
         // of relying on the caller-provided delta so incremental scans remain in sync with the snapshot
         // without computing additional cycles or mutating cache state.
         _currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime);
-        _previousExponent = exponent;
+        PreviousExponent = exponent;
         isUnity = _currentMontgomery == _montgomeryOne;
         return true;
     }
@@ -159,13 +135,7 @@ internal struct ExponentRemainderStepper
     private void InitializeCpuState(ulong exponent)
     {
         _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
-        _previousExponent = exponent;
+        PreviousExponent = exponent;
-    }
-
-    private void InitializeGpuState(ulong exponent)
-    {
-        _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true);
-        _previousExponent = exponent;
     }
 
     private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime);

PerfectNumbers.Core/Gpu/ExponentRemainderStepperGpu.cs

@@ -0,0 +1,112 @@
+using System.Runtime.CompilerServices;
+using PerfectNumbers.Core;
+
+namespace PerfectNumbers.Core.Gpu;
+
+/// <summary>
+/// Tracks successive residues for a single divisor while kernels advance exponent candidates on the GPU path.
+/// Callers must seed the stepper once per divisor and then supply exponents in strictly ascending order.
+/// </summary>
+internal struct ExponentRemainderStepperGpu
+{
+    private readonly ulong _modulus;
+    private readonly ReadOnlyGpuUInt128 _modulusWide;
+    public ulong PreviousExponent;
+    private GpuUInt128 _currentResidue;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public ExponentRemainderStepperGpu(in GpuDivisorPartialData divisor)
+    {
+        ulong modulus = divisor.Modulus;
+        _modulus = modulus;
+        _modulusWide = divisor.ModulusWide;
+        PreviousExponent = 0UL;
+        _currentResidue = GpuUInt128.One;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public bool MatchesDivisor(in GpuDivisorPartialData divisor)
+    {
+        return _modulus == divisor.Modulus;
+    }
+
+    /// <summary>
+    /// Clears the cached residue so the stepper can be reused for another divisor.
+    /// Callers must reinitialize the stepper before consuming any residues after a reset.
+    /// </summary>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public void Reset()
+    {
+        PreviousExponent = 0UL;
+        _currentResidue = GpuUInt128.One;
+    }
+
+    /// <summary>
+    /// Initializes the stepper for GPU calculations using the first exponent in a strictly ascending sequence.
+    /// The returned value is the canonical residue corresponding to the provided exponent.
+    /// </summary>
+    /// <param name="exponent">First exponent evaluated for the current divisor. The caller must reuse the same divisor metadata.</param>
+    /// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public ulong InitializeGpu(ulong exponent)
+    {
+        InitializeStateGpu(exponent);
+        return _currentResidue.Low;
+    }
+
+    /// <summary>
+    /// Initializes the GPU stepping state and reports whether the seed exponent already yields unity.
+    /// </summary>
+    /// <param name="exponent">First exponent evaluated for the current divisor. Must match the divisor passed to the constructor.</param>
+    /// <returns><see langword="true"/> when the seed exponent produces a residue equal to one.</returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public bool InitializeIsUnityGpu(ulong exponent)
+    {
+        InitializeStateGpu(exponent);
+        return _currentResidue.Low == 1UL;
+    }
+
+    /// <summary>
+    /// Advances the GPU stepping state to the next exponent in a strictly ascending sequence that shares the same divisor.
+    /// Callers must invoke <see cref="InitializeGpu(ulong)"/> beforehand.
+    /// </summary>
+    /// <param name="exponent">Next exponent in the ascending sequence. Callers must never move backwards or reuse an older exponent.</param>
+    /// <returns>The canonical residue produced by <c>2^exponent mod divisor</c>.</returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public ulong ComputeNextGpu(ulong exponent)
+    {
+        AdvanceStateGpu(exponent);
+        return _currentResidue.Low;
+    }
+
+    /// <summary>
+    /// Advances the GPU stepping state to the next exponent and reports whether it evaluates to unity.
+    /// Callers must seed the stepper through <see cref="InitializeGpu(ulong)"/> or <see cref="InitializeIsUnityGpu(ulong)"/> and then
+    /// pass exponents in strictly ascending order.
+    /// </summary>
+    /// <param name="exponent">Next exponent in the ascending sequence.</param>
+    /// <returns><see langword="true"/> when the canonical residue equals one.</returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public bool ComputeNextIsUnityGpu(ulong exponent)
+    {
+        AdvanceStateGpu(exponent);
+        return _currentResidue.Low == 1UL;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private void InitializeStateGpu(ulong exponent)
+    {
+        ulong residue = ULongExtensions.Pow2ModWindowedGpuKernel(exponent, _modulus);
+        _currentResidue = new GpuUInt128(residue);
+        PreviousExponent = exponent;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private void AdvanceStateGpu(ulong exponent)
+    {
+        ulong delta = exponent - PreviousExponent;
+        ulong multiplier = ULongExtensions.Pow2ModWindowedGpuKernel(delta, _modulus);
+        _currentResidue.MulMod(multiplier, _modulusWide);
+        PreviousExponent = exponent;
+    }
+}

PerfectNumbers.Core/Gpu/GpuDivisorPartialData.cs

@@ -0,0 +1,23 @@
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core.Gpu;
+
+internal readonly struct GpuDivisorPartialData
+{
+    public readonly ulong Modulus;
+    public readonly ReadOnlyGpuUInt128 ModulusWide;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuDivisorPartialData(ulong modulus)
+    {
+        Modulus = modulus;
+        ModulusWide = new ReadOnlyGpuUInt128(modulus);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuDivisorPartialData(ulong modulus, in ReadOnlyGpuUInt128 modulusWide)
+    {
+        Modulus = modulus;
+        ModulusWide = modulusWide;
+    }
+}

PerfectNumbers.Core/Gpu/GpuUInt128.cs

@@ -11,9 +11,7 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
 
     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();
 
     private const int Pow2WindowSizeBits = 8;
     private const int Pow2WindowOddPowerCount = 1 << (Pow2WindowSizeBits - 1);
@@ -902,23 +900,6 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public readonly ulong MulMod(ulong value, ulong modulus)
-    {
-        ulong a = Low % modulus;
-        ulong b = value % modulus;
-
-        if (a == 0UL || b == 0UL)
-        {
-            return 0UL;
-        }
-
-        ulong ulongRange = ulong.MaxValue / a;
-        return b <= ulongRange
-            ? (a * b) % modulus
-            : MulMod64(a, b, modulus);
-    }
-
-    [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    public readonly ulong MulModSimplified(ulong value, ulong modulus)
     {
         ulong modulusLocal = modulus;
         ulong a = Low % modulusLocal;
@@ -1028,10 +1009,10 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
     private static ulong MultiplyChunkModulo(ulong value, ulong chunk, ulong modulus)
     {
         ulong result = 0UL;
-        ulong[] nativeModuloBitMasks = NativeModuloBitMasks;
+        ulong bitMask = 1UL;
-        for (int bit = 0; bit < NativeModuloChunkBits; bit++)
+        for (int bit = 0; bit < NativeModuloChunkBits; bit++, bitMask <<= 1)
         {
-            if ((chunk & nativeModuloBitMasks[bit]) != 0UL)
+            if ((chunk & bitMask) != 0UL)
             {
                 result = (result + value) % modulus;
             }
@@ -1548,16 +1529,6 @@ public struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
         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: Check if the TODO below is still relevant.

PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs

@@ -7,280 +7,82 @@ namespace PerfectNumbers.Core.Gpu;
 
 internal static class DivisorByDivisorKernels
 {
-    public static void CheckKernel(Index1D index, ArrayView<MontgomeryDivisorData> divisors, ArrayView<ulong> exponents, ArrayView<byte> hits)
+    public static void CheckKernel(
-    {
-        MontgomeryDivisorData divisor = divisors[index];
-        ulong modulus = divisor.Modulus;
-        // The GPU by-divisor pipeline only materializes odd moduli greater than one (q = 2kp + 1),
-        // so the defensive guard for invalid values stays disabled here to keep the kernel branch-free.
-        ulong exponent = exponents[index];
-        ulong montgomeryResult = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent);
-        hits[index] = montgomeryResult == 1UL ? (byte)1 : (byte)0;
-    }
-
-    public static void GenerateCandidatesAndGapsKernel(
-        Index1D index,
-        ArrayView<ulong> candidates,
-        ArrayView<ulong> gaps,
-        ArrayView<GpuUInt128> nextStarts,
-        GpuUInt128 startValue,
-        GpuUInt128 stride)
-    {
-        int globalIndex = index;
-        // EvenPerfectBitScanner always launches this kernel with an extent matching the candidate view length,
-        // so the bounds guard stays commented out to keep the hot path branch-free.
-        // int candidateLength = (int)candidates.Length;
-        // if (globalIndex >= candidateLength)
-        // {
-        //     return;
-        // }
-
-        ulong offset = (ulong)globalIndex;
-        GpuUInt128 candidateValue = stride;
-        candidateValue.Mul(offset);
-        candidateValue.Add(startValue);
-        candidates[globalIndex] = candidateValue.Low;
-
-        // The scanner supplies gap buffers that mirror the candidate view length,
-        // so the additional bounds guard remains commented out for the same reason as above.
-        // int gapLength = (int)gaps.Length;
-        // if (globalIndex >= gapLength)
-        // {
-        //     return;
-        // }
-
-        ulong previousOffset = offset > 0UL ? offset - 1UL : 0UL;
-        GpuUInt128 previousValue = stride;
-        previousValue.Mul(previousOffset);
-        previousValue.Add(startValue);
-
-        GpuUInt128 gapValue = candidateValue;
-        gapValue.Sub(previousValue);
-        gaps[globalIndex] = gapValue.Low;
-
-        if (globalIndex == 0)
-        {
-            // EvenPerfectBitScanner provides a single-element nextStarts view per launch, so the guard stays commented out.
-            // int nextStartLength = (int)nextStarts.Length;
-            // if (nextStartLength <= 0)
-            // {
-            //     return;
-            // }
-
-            ulong candidateCount = (ulong)candidates.Length;
-            GpuUInt128 nextValue = stride;
-            nextValue.Mul(candidateCount);
-            nextValue.Add(startValue);
-            nextStarts[0] = nextValue;
-        }
-    }
-
-    public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView<ulong> gaps, ArrayView<byte> deltas, byte modulus)
-    {
-        deltas[index] = (byte)(gaps[index] % modulus);
-    }
-
-    public static void AccumulateRemaindersKernel(
-        Index1D index,
-        ArrayView<byte> deltas,
-        ArrayView<byte> remainders,
-        byte baseRemainder,
-        byte modulus)
-    {
-        int globalIndex = index;
-        // EvenPerfectBitScanner launches this stream kernel with an extent matching the remainder span length,
-        // so the defensive guard stays commented out to keep the kernel branch-free.
-        // int length = (int)remainders.Length;
-        // if (globalIndex >= length)
-        // {
-        //     return;
-        // }
-
-        int localIndex = Group.IdxX;
-        int groupSize = Group.Dimension.X;
-        var shared = SharedMemory.GetDynamic<int>();
-
-        shared[localIndex] = deltas[globalIndex];
-        Group.Barrier();
-
-        int offset = 1;
-        int mod = modulus;
-        while (offset < groupSize)
-        {
-            bool apply = localIndex >= offset;
-            int addend = apply ? shared[localIndex - offset] : 0;
-
-            Group.Barrier();
-
-            int current = shared[localIndex];
-            int sum = apply ? current + addend : current;
-            int wrap = sum >= mod ? 1 : 0;
-            sum -= wrap * mod;
-            shared[localIndex] = apply ? sum : current;
-
-            Group.Barrier();
-
-            offset <<= 1;
-        }
-
-        int remainder = baseRemainder + shared[localIndex];
-        int remainderWrap = remainder >= mod ? 1 : 0;
-        remainder -= remainderWrap * mod;
-        remainders[globalIndex] = (byte)remainder;
-    }
-
-    public static void EvaluateCandidateMaskKernel(
-        Index1D index,
-        ArrayView<byte> remainder10,
-        ArrayView<byte> remainder8,
-        ArrayView<byte> remainder5,
-        ArrayView<byte> remainder3,
-        byte lastIsSevenFlag,
-        ArrayView<byte> mask)
-    {
-        int globalIndex = index;
-        // EvenPerfectBitScanner launches this auto-grouped kernel with an extent matching the mask length,
-        // so the defensive guard stays commented out to keep the kernel branch-free.
-        // int length = (int)mask.Length;
-        // if (globalIndex >= length)
-        // {
-        //     return;
-        // }
-
-        byte value10 = remainder10[globalIndex];
-        bool lastIsSeven = lastIsSevenFlag != 0;
-        bool acceptSeven = value10 == 3 || value10 == 7 || value10 == 9;
-        bool acceptNonSeven = value10 == 1 || value10 == 3 || value10 == 7 || value10 == 9;
-        bool accept10 = lastIsSeven ? acceptSeven : acceptNonSeven;
-        byte value8 = remainder8[globalIndex];
-        bool accept8 = value8 == 1 || value8 == 7;
-        bool accept5 = remainder5[globalIndex] != 0;
-        bool accept3 = remainder3[globalIndex] != 0;
-        byte accepted = (byte)((accept10 && accept8 && accept3 && accept5) ? 1 : 0);
-        mask[globalIndex] = accepted;
-    }
-
-    public static void CompactCandidateMaskKernel(
-        Index1D index,
-        ArrayView<byte> mask,
-        ArrayView<ulong> sourceCandidates,
-        ArrayView<ulong> compactedCandidates,
-        ArrayView<int> compactedCount)
-    {
-        int globalIndex = index;
-        int length = (int)mask.Length;
-        // EvenPerfectBitScanner launches this stream kernel with an extent matching the mask length,
-        // so the defensive guard stays commented out to keep the kernel branch-free.
-        // if (globalIndex >= length)
-        // {
-        //     return;
-        // }
-
-        int localIndex = Group.IdxX;
-        int groupSize = Group.Dimension.X;
-        var shared = SharedMemory.GetDynamic<int>();
-
-        int prefixBase = 0;
-        int mappingBase = length;
-        int fallbackBase = length * 2;
-        // The dynamic shared-memory slice stores prefix sums, compacted index mappings, and fallback slots
-        // in separate regions so the scatter stage can avoid divergent branches.
-
-        byte maskValue = mask[globalIndex];
-        int accepted = maskValue != 0 ? 1 : 0;
-        shared[prefixBase + localIndex] = accepted;
-        shared[fallbackBase + localIndex] = globalIndex;
-
-        Group.Barrier();
-
-        int offset = 1;
-        while (offset < groupSize)
-        {
-            int addendIndex = prefixBase + localIndex - offset;
-            int addend = localIndex >= offset ? shared[addendIndex] : 0;
-
-            Group.Barrier();
-
-            int currentIndex = prefixBase + localIndex;
-            int current = shared[currentIndex];
-            int sum = current + addend;
-            int shouldUpdate = localIndex >= offset ? 1 : 0;
-            shared[currentIndex] = shouldUpdate != 0 ? sum : current;
-
-            Group.Barrier();
-
-            offset <<= 1;
-        }
-
-        int inclusivePrefix = shared[prefixBase + localIndex];
-        int exclusiveIndex = inclusivePrefix - accepted;
-        int targetBase = accepted != 0 ? mappingBase : fallbackBase;
-        int targetOffset = accepted != 0 ? exclusiveIndex : localIndex;
-        shared[targetBase + targetOffset] = globalIndex;
-
-        Group.Barrier();
-
-        int totalAccepted = shared[prefixBase + length - 1];
-        int writerMask = localIndex < totalAccepted ? 1 : 0;
-        int lookupBase = writerMask != 0 ? mappingBase : fallbackBase;
-        int sourceIndex = shared[lookupBase + localIndex];
-
-        ulong candidateValue = sourceCandidates[sourceIndex];
-        int destinationIndex = localIndex;
-        ulong existingValue = compactedCandidates[destinationIndex];
-        compactedCandidates[destinationIndex] = writerMask != 0 ? candidateValue : existingValue;
-
-        int lastLaneMask = globalIndex == length - 1 ? 1 : 0;
-        int countValue = shared[prefixBase + localIndex];
-        if (lastLaneMask != 0)
-        {
-            compactedCount[0] = countValue;
-        }
-    }
-
-    public static void PublishFirstHitKernel(
         Index1D index,
+        ArrayView<GpuDivisorPartialData> divisors,
+        ArrayView<int> exponentOffsets,
+        ArrayView<int> exponentCounts,
+        ArrayView<ulong> exponents,
+        ArrayView<ulong> divisorCycles,
+        ArrayView<ulong> firstCycleRemainders,
         ArrayView<byte> hits,
         ArrayView<int> firstHit)
     {
         int globalIndex = index;
-        // EvenPerfectBitScanner launches this auto-grouped kernel with an extent matching the hit view length,
+        int count = exponentCounts[globalIndex];
-        // so the defensive guard stays commented out to keep the hot path branch-free.
+        if (count <= 0)
-        // int length = (int)hits.Length;
+        {
-        // if (globalIndex >= length)
+            return;
-        // {
+        }
-        //     return;
-        // }
 
-        byte hitValue = hits[globalIndex];
+        int offset = exponentOffsets[globalIndex];
-        int candidate = hitValue != 0 ? globalIndex : int.MaxValue;
+        GpuDivisorPartialData divisor = divisors[globalIndex];
-        Atomic.Min(ref firstHit[0], candidate);
+        ulong cycleLength = divisorCycles[globalIndex];
-    }
+        ulong firstCycleRemainder = firstCycleRemainders[globalIndex];
 
-    public static void ConvertMontgomeryResultsToHitsKernel(
+        var stepper = new ExponentRemainderStepperGpu(divisor);
-        Index1D index,
-        ArrayView<ulong> results,
-        ArrayView<byte> hits)
-    {
-        int globalIndex = index;
-        // EvenPerfectBitScanner launches this kernel with an extent matching the result and hit spans,
-        // so the defensive bounds guard stays commented out to keep the conversion branch-free.
-        // int length = (int)results.Length;
-        // if (globalIndex >= length)
-        // {
-        //     return;
-        // }
 
-        ulong montgomeryResult = results[globalIndex];
+        ulong firstExponent = exponents[offset];
-        hits[globalIndex] = montgomeryResult == 1UL ? (byte)1 : (byte)0;
+        bool firstUnity = stepper.InitializeIsUnityGpu(firstExponent);
-    }
 
-    public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
+        byte firstHitValue;
-    {
+        ulong previousExponent = firstExponent;
-        ulong modulus = divisor.Modulus;
+        ulong cycleRemainder = 0UL;
-        // Each divisor flowing through this kernel follows q = 2kp + 1 with k >= 1, so modulus is always odd and exceeds one.
+        if (cycleLength == 0UL)
-        ulong exponent = exponents[index];
+        {
-        ulong montgomeryResult = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponent);
+            firstHitValue = firstUnity ? (byte)1 : (byte)0;
-        results[index] = montgomeryResult;
+        }
+        else
+        {
+            cycleRemainder = firstCycleRemainder;
+            firstHitValue = cycleRemainder == 0UL && firstUnity ? (byte)1 : (byte)0;
+        }
+
+        hits[offset] = firstHitValue;
+        byte anyHit = firstHitValue;
+
+        for (int i = 1; i < count; i++)
+        {
+            ulong exponent = exponents[offset + i];
+            ulong delta = exponent - previousExponent;
+            previousExponent = exponent;
+
+            if (cycleLength != 0UL)
+            {
+                var remainderValue = new GpuUInt128(cycleRemainder);
+                remainderValue.AddMod(delta, cycleLength);
+                cycleRemainder = remainderValue.Low;
+
+                byte hit = 0;
+                if (cycleRemainder == 0UL)
+                {
+                    hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0;
+                }
+
+                hits[offset + i] = hit;
+                anyHit |= hit;
+            }
+            else
+            {
+                byte hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0;
+                hits[offset + i] = hit;
+                anyHit |= hit;
+            }
+        }
+
+        if (anyHit != 0 && firstHit.Length > 0)
+        {
+            Atomic.Min(ref firstHit[0], globalIndex);
+        }
     }
 }

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -81,7 +81,7 @@ public class MersenneDivisorCycles
             return heuristicCycle == exponent;
         }
 
-        if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL)
+        if (TryCalculateCycleLengthForExponentCpu(divisor, exponent, divisorData, out ulong cycleLength, out bool _) && cycleLength != 0UL)
         {
             return cycleLength == exponent;
         }
@@ -91,10 +91,12 @@ public class MersenneDivisorCycles
 
     public static bool CycleEqualsExponentForMersenneCandidate(ulong divisor, in MontgomeryDivisorData divisorData, ulong exponent)
     {
-        if (!IsValidMersenneDivisorCandidate(divisor, exponent))
+        // EvenPerfectBitScanner only calls this helper with divisors generated from 2 * k * p + 1 and primes greater than one,
-        {
+        // so the validation guard would never trigger in production. Keep it commented out to leave the hot path branch-free.
-            return false;
+        // if (!IsValidMersenneDivisorCandidate(divisor, exponent))
-        }
+        // {
+        //     return false;
+        // }
 
         if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles)
         {
@@ -245,7 +247,7 @@ public class MersenneDivisorCycles
         return order;
     }
 
-    public static bool TryCalculateCycleLengthForExponent(
+    public static bool TryCalculateCycleLengthForExponentCpu(
         ulong divisor,
         ulong exponent,
         in MontgomeryDivisorData divisorData,
@@ -255,17 +257,19 @@ public class MersenneDivisorCycles
         cycleLength = 0UL;
         primeOrderFailed = false;
 
-        if ((divisor & (divisor - 1UL)) == 0UL)
+        // EvenPerfectBitScanner never feeds powers of two into this helper, so leave the branch commented out to avoid redundant work.
-        {
+        // if ((divisor & (divisor - 1UL)) == 0UL)
-            cycleLength = 1UL;
+        // {
-            return true;
+        //     cycleLength = 1UL;
-        }
+        //     return true;
+        // }
 
-        if (divisor <= 3UL || exponent <= 1UL)
+        // The scanner only calls this path with prime exponents greater than one and odd divisors above three, keeping the small-parameter guard unreachable.
-        {
+        // if (divisor <= 3UL || exponent <= 1UL)
-            cycleLength = CalculateCycleLength(divisor, divisorData);
+        // {
-            return true;
+        //     cycleLength = CalculateCycleLength(divisor, divisorData);
-        }
+        //     return true;
+        // }
 
         ulong computedOrder = PrimeOrderCalculator.Calculate(
             divisor,
@@ -482,7 +486,7 @@ public class MersenneDivisorCycles
         ulong[]? heapCandidateArray = null;
         bool[]? heapEvaluationArray = null;
 
-        ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+        ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
 
         ulong order = initialOrder;
         for (int i = 0; i < count; i++)
@@ -513,7 +517,7 @@ public class MersenneDivisorCycles
             ProcessReduceOrderPrime(heapCandidateArray.AsSpan(0, multiplicity), heapEvaluationArray.AsSpan(0, multiplicity), ref order, prime, multiplicity, ref stepper);
         }
 
-        ThreadStaticPools.ReturnExponentStepper(stepper);
+        ThreadStaticPools.ReturnExponentStepperCpu(stepper);
         ThreadStaticPools.UlongPool.Return(primes);
 
         if (heapCandidateArray is not null)
@@ -525,7 +529,7 @@ public class MersenneDivisorCycles
         return order;
     }
 
-    private static void ProcessReduceOrderPrime(Span<ulong> candidateBuffer, Span<bool> evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepper stepper)
+    private static void ProcessReduceOrderPrime(Span<ulong> candidateBuffer, Span<bool> evaluationBuffer, ref ulong order, ulong prime, int multiplicity, ref ExponentRemainderStepperCpu stepper)
     {
         ulong working = order;
         int actual = 0;

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs

@@ -163,12 +163,12 @@ internal static partial class PrimeOrderCalculator
 		Span<ulong> candidates = buffer[..actual];
 		candidates.Sort();
 
-		ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+		ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
 		int candidateCount = candidates.Length;
 
 		if (stepper.InitializeCpuIsUnity(candidates[0]))
 		{
-			ThreadStaticPools.ReturnExponentStepper(stepper);
+			ThreadStaticPools.ReturnExponentStepperCpu(stepper);
 			return false;
 		}
 
@@ -176,12 +176,12 @@ internal static partial class PrimeOrderCalculator
 		{
 			if (stepper.ComputeNextIsUnity(candidates[i]))
 			{
-				ThreadStaticPools.ReturnExponentStepper(stepper);
+				ThreadStaticPools.ReturnExponentStepperCpu(stepper);
 				return false;
 			}
 		}
 
-		ThreadStaticPools.ReturnExponentStepper(stepper);
+		ThreadStaticPools.ReturnExponentStepperCpu(stepper);
 		return true;
 	}
 
@@ -219,7 +219,7 @@ internal static partial class PrimeOrderCalculator
 
 		buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value));
 
-		ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+		ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
 
 		const int StackExponentCapacity = 16;
 		const int ExponentHardLimit = 256;
@@ -256,7 +256,7 @@ internal static partial class PrimeOrderCalculator
 			ProcessExponentLoweringPrime(heapCandidateArray.AsSpan(0, exponent), heapEvaluationArray.AsSpan(0, exponent), ref order, primeFactor, exponent, ref stepper);
 		}
 
-		ThreadStaticPools.ReturnExponentStepper(stepper);
+		ThreadStaticPools.ReturnExponentStepperCpu(stepper);
 
 		if (heapCandidateArray is not null)
 		{
@@ -268,7 +268,7 @@ internal static partial class PrimeOrderCalculator
 		return order;
 	}
 
-	private static void ProcessExponentLoweringPrime(Span<ulong> candidateBuffer, Span<bool> evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepper stepper)
+	private static void ProcessExponentLoweringPrime(Span<ulong> candidateBuffer, Span<bool> evaluationBuffer, ref ulong order, ulong primeFactor, int exponent, ref ExponentRemainderStepperCpu stepper)
 	{
 		ulong working = order;
 		int actual = 0;
@@ -388,7 +388,7 @@ internal static partial class PrimeOrderCalculator
 		Span<ulong> stackBuffer = stackalloc ulong[StackExponentCapacity];
 
 		ulong[]? heapCandidateArray = null;
-		ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+		ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
 		bool violates = false;
 
 		for (int i = 0; i < length; i++)
@@ -425,7 +425,7 @@ internal static partial class PrimeOrderCalculator
 			}
 		}
 
-		ThreadStaticPools.ReturnExponentStepper(stepper);
+		ThreadStaticPools.ReturnExponentStepperCpu(stepper);
 
 		if (heapCandidateArray is not null)
 		{
@@ -435,7 +435,7 @@ internal static partial class PrimeOrderCalculator
 		return !violates;
 	}
 
-	private static bool ValidateOrderForFactor(Span<ulong> buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepper stepper)
+	private static bool ValidateOrderForFactor(Span<ulong> buffer, ulong primeFactor, int exponent, ulong order, ref ExponentRemainderStepperCpu stepper)
 	{
 		ulong working = order;
 		int actual = 0;
@@ -547,7 +547,7 @@ internal static partial class PrimeOrderCalculator
 			int candidateCount = candidates.Count;
 			bool allowGpuBatch = true;
 			Span<ulong> candidateSpan = CollectionsMarshal.AsSpan(candidates);
-			ExponentRemainderStepper powStepper = ThreadStaticPools.RentExponentStepper(divisorData);
+			ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
 			bool powStepperInitialized = false;
 
 			// DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)");
@@ -626,7 +626,7 @@ internal static partial class PrimeOrderCalculator
 
 						candidates.Clear();
 						ThreadStaticPools.ReturnUlongList(candidates);
-						ThreadStaticPools.ReturnExponentStepper(powStepper);
+						ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
 						result = candidate;
 						return true;
 					}
@@ -666,7 +666,7 @@ internal static partial class PrimeOrderCalculator
 
 					candidates.Clear();
 					ThreadStaticPools.ReturnUlongList(candidates);
-					ThreadStaticPools.ReturnExponentStepper(powStepper);
+					ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
 					result = candidate;
 					return true;
 				}
@@ -681,7 +681,7 @@ internal static partial class PrimeOrderCalculator
 
 			candidates.Clear();
 			ThreadStaticPools.ReturnUlongList(candidates);
-			ThreadStaticPools.ReturnExponentStepper(powStepper);
+			ThreadStaticPools.ReturnExponentStepperCpu(powStepper);
 			// DebugLog("No candidate confirmed");
 			return false;
 		}
@@ -892,7 +892,7 @@ internal static partial class PrimeOrderCalculator
 
 			Span<ulong> stackBuffer = stackalloc ulong[StackExponentCapacity];
 			ulong[]? heapCandidateArray = null;
-			ExponentRemainderStepper stepper = ThreadStaticPools.RentExponentStepper(divisorData);
+			ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData);
 			bool violates = false;
 
 			for (int i = 0; i < length; i++)
@@ -929,7 +929,7 @@ internal static partial class PrimeOrderCalculator
 				}
 			}
 
-			ThreadStaticPools.ReturnExponentStepper(stepper);
+			ThreadStaticPools.ReturnExponentStepperCpu(stepper);
 
 			if (heapCandidateArray is not null)
 			{
@@ -944,7 +944,7 @@ internal static partial class PrimeOrderCalculator
 		}
 	}
 
-	private static bool CheckCandidateViolation(Span<ulong> buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepper stepper)
+	private static bool CheckCandidateViolation(Span<ulong> buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepperCpu stepper)
 	{
 		ulong working = candidate;
 		int actual = 0;

PerfectNumbers.Core/ThreadStaticPools.cs

@@ -2,6 +2,7 @@ using System.Buffers;
 using System.Collections.Generic;
 using System.Runtime.CompilerServices;
 using PerfectNumbers.Core.Cpu;
+using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core
 {
@@ -102,6 +103,18 @@ namespace PerfectNumbers.Core
             }
         }
 
+        [ThreadStatic]
+        private static ArrayPool<GpuDivisorPartialData>? _gpuDivisorPartialDataPool;
+
+        internal static ArrayPool<GpuDivisorPartialData> GpuDivisorPartialDataPool
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get
+            {
+                return _gpuDivisorPartialDataPool ??= ArrayPool<GpuDivisorPartialData>.Create();
+            }
+        }
+
         [ThreadStatic]
         private static ArrayPool<int>? _intPool;
 
@@ -329,31 +342,31 @@ namespace PerfectNumbers.Core
         }
 
         [ThreadStatic]
-        private static ExponentRemainderStepper _exponentRemainderStepper;
+        private static ExponentRemainderStepperCpu _exponentRemainderStepperCpu;
 
         [ThreadStatic]
-        private static bool _hasExponentRemainderStepper;
+        private static bool _hasExponentRemainderStepperCpu;
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static ExponentRemainderStepper RentExponentStepper(in MontgomeryDivisorData divisorData)
+        internal static ExponentRemainderStepperCpu RentExponentStepperCpu(in MontgomeryDivisorData divisorData)
         {
-            if (_hasExponentRemainderStepper && _exponentRemainderStepper.MatchesDivisor(divisorData))
+            if (_hasExponentRemainderStepperCpu && _exponentRemainderStepperCpu.MatchesDivisor(divisorData))
             {
-                ExponentRemainderStepper stepper = _exponentRemainderStepper;
+                ExponentRemainderStepperCpu stepper = _exponentRemainderStepperCpu;
-                _hasExponentRemainderStepper = false;
+                _hasExponentRemainderStepperCpu = false;
                 stepper.Reset();
                 return stepper;
             }
 
-            return new ExponentRemainderStepper(divisorData);
+            return new ExponentRemainderStepperCpu(divisorData);
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static void ReturnExponentStepper(ExponentRemainderStepper stepper)
+        internal static void ReturnExponentStepperCpu(ExponentRemainderStepperCpu stepper)
         {
             stepper.Reset();
-            _exponentRemainderStepper = stepper;
+            _exponentRemainderStepperCpu = stepper;
-            _hasExponentRemainderStepper = true;
+            _hasExponentRemainderStepperCpu = true;
         }
     }
 }

PerfectNumbers.Core/ULongExtensions.Gpu.cs

@@ -1,3 +1,4 @@
+using System;
 using System.Collections.Concurrent;
 using System.Runtime.CompilerServices;
 using ILGPU;
@@ -18,6 +19,7 @@ public static partial class ULongExtensions
 	}
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	[Obsolete("Use ULongExtensions.MulModGpu for GPU-compatible host code or GpuUInt128.MulMod inside kernels.")]
 	public static ulong MulMod64Gpu(this ulong a, ulong b, ulong modulus)
 	{
 		// TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64,
@@ -27,6 +29,7 @@ public static partial class ULongExtensions
 	}
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	[Obsolete("Use ULongExtensions.MulModGpu for GPU-compatible host code or GpuUInt128.MulMod inside kernels.")]
 	public static ulong MulMod64GpuDeferred(this ulong a, ulong b, ulong modulus)
 	{
 		// TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the
@@ -35,6 +38,13 @@ public static partial class ULongExtensions
 		return state.MulModWithNativeModulo(b, modulus);
 	}
 
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong MulModGpu(this ulong left, ulong right, ulong modulus)
+	{
+		GpuUInt128 state = new(left);
+		return state.MulMod(right, modulus);
+	}
+
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	public static ulong Pow2ModWindowedGpu(this ulong exponent, ulong modulus)
 	{
@@ -71,7 +81,7 @@ public static partial class ULongExtensions
 		while (index >= 0)
 		{
 			ulong currentBit = (exponent >> index) & 1UL;
-			ulong squared = result.MulMod64Gpu(result, modulus);
+			ulong squared = result.MulModGpu(result, modulus);
 			result = currentBit == 0UL ? squared : result;
 			index = currentBit == 0UL ? index - 1 : index;
 			if (currentBit == 0UL)
@@ -91,14 +101,14 @@ public static partial class ULongExtensions
 			int windowLength = index - windowStart + 1;
 			for (int square = 0; square < windowLength; square++)
 			{
-				result = result.MulMod64Gpu(result, modulus);
+				result = result.MulModGpu(result, modulus);
 			}
 
 			ulong mask = (1UL << windowLength) - 1UL;
 			ulong windowValue = (exponent >> windowStart) & mask;
 			int tableIndex = (int)((windowValue - 1UL) >> 1);
 			ulong multiplier = oddPowers[tableIndex];
-			result = result.MulMod64Gpu(multiplier, modulus);
+			result = result.MulModGpu(multiplier, modulus);
 
 			index = windowStart - 1;
 		}
@@ -108,6 +118,79 @@ public static partial class ULongExtensions
 		return result;
 	}
 
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	internal static ulong Pow2ModWindowedGpuKernel(ulong exponent, ulong modulus)
+	{
+		int bitLength = GetPortableBitLengthGpu(exponent);
+		int windowSize = GetWindowSizeGpu(bitLength);
+		ReadOnlyGpuUInt128 modulusWide = new(modulus);
+		GpuUInt128 result = new(1UL);
+		int index = bitLength - 1;
+
+		while (index >= 0)
+		{
+			ulong currentBit = (exponent >> index) & 1UL;
+			if (currentBit == 0UL)
+			{
+				result.MulMod(result.AsReadOnly(), modulusWide);
+				index--;
+				continue;
+			}
+
+			int windowStartCandidate = index - windowSize + 1;
+			int negativeMask = windowStartCandidate >> 31;
+			windowStartCandidate &= ~negativeMask;
+			int windowStart = GetNextSetBitIndexGpu(exponent, windowStartCandidate);
+			int windowLength = index - windowStart + 1;
+
+			for (int square = 0; square < windowLength; square++)
+			{
+				result.MulMod(result.AsReadOnly(), modulusWide);
+			}
+
+			ulong mask = (1UL << windowLength) - 1UL;
+			ulong windowValue = (exponent >> windowStart) & mask;
+			ulong multiplier = ComputeWindowedOddPowerGpuKernel(windowValue, modulusWide);
+			result.MulMod(multiplier, modulusWide);
+
+			index = windowStart - 1;
+		}
+
+		return result.Low;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static ulong ComputeWindowedOddPowerGpuKernel(ulong windowValue, in ReadOnlyGpuUInt128 modulus)
+	{
+		GpuUInt128 result = new(2UL);
+		if (windowValue == 1UL)
+		{
+			// Single-bit windows show up whenever the bit scan finds an isolated exponent bit,
+			// which happens routinely on the GPU path, so keep the fast return enabled.
+			return 2UL;
+		}
+
+		ulong remaining = (windowValue - 1UL) >> 1;
+		GpuUInt128 squareBase = new(2UL);
+		squareBase.MulMod(squareBase.AsReadOnly(), modulus);
+
+		while (remaining != 0UL)
+		{
+			if ((remaining & 1UL) != 0UL)
+			{
+				result.MulMod(squareBase.AsReadOnly(), modulus);
+			}
+
+			remaining >>= 1;
+			if (remaining != 0UL)
+			{
+				squareBase.MulMod(squareBase.AsReadOnly(), modulus);
+			}
+		}
+
+		return result.Low;
+	}
+
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	internal static ulong Pow2MontgomeryModWindowedGpuKeepMontgomery(in MontgomeryDivisorData divisor, ulong exponent)
 	{
@@ -228,11 +311,11 @@ public static partial class ULongExtensions
 			return;
 		}
 
-		ulong square = baseValue.MulMod64Gpu(baseValue, modulus);
+		ulong square = baseValue.MulModGpu(baseValue, modulus);
 		for (int i = 1; i < oddPowers.Length; i++)
 		{
 			ulong previous = oddPowers[i - 1];
-			oddPowers[i] = previous.MulMod64Gpu(square, modulus);
+			oddPowers[i] = previous.MulModGpu(square, modulus);
 		}
 	}