Align GPU pow2 heuristics with Montgomery helper

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -4,8 +4,6 @@ using ILGPU;
 using ILGPU.Runtime;
 using PerfectNumbers.Core;
 using System.Runtime.InteropServices;
-using System.Numerics;
-
 namespace PerfectNumbers.Core.Gpu;
 
 internal enum GpuPow2ModStatus
@@ -18,7 +16,7 @@ internal enum GpuPow2ModStatus
 internal static class PrimeOrderGpuHeuristics
 {
     private static readonly ConcurrentDictionary<ulong, byte> OverflowedPrimes = new();
-    private static readonly ConcurrentDictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, GpuUInt128, ArrayView1D<ulong, Stride1D.Dense>>> Pow2ModKernelCache = new();
+    private static readonly ConcurrentDictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>>> Pow2ModKernelCache = new();
     private static PrimeOrderGpuCapability s_capability = PrimeOrderGpuCapability.Default;
 
     internal static ConcurrentDictionary<ulong, byte> OverflowRegistry => OverflowedPrimes;
@@ -108,7 +106,8 @@ internal static class PrimeOrderGpuHeuristics
                 exponentBuffer.View.CopyFromCPU(ref MemoryMarshal.GetReference(exponents), exponents.Length);
                 remainderBuffer.MemSetToZero();
 
-                kernel(stream, exponents.Length, exponentBuffer.View, new GpuUInt128(prime), remainderBuffer.View);
+                MontgomeryDivisorData divisor = MontgomeryDivisorDataCache.Get(prime);
+                kernel(stream, exponents.Length, exponentBuffer.View, divisor, remainderBuffer.View);
 
                 stream.Synchronize();
 
@@ -142,44 +141,24 @@ internal static class PrimeOrderGpuHeuristics
             return 0UL;
         }
 
-        ulong result = 1UL % modulus;
+        MontgomeryDivisorData divisor = MontgomeryDivisorDataCache.Get(modulus);
-        ulong baseValue = 2UL % modulus;
+        return exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
-        ulong remaining = exponent;
-
-        while (remaining > 0UL)
-        {
-            if ((remaining & 1UL) != 0UL)
-            {
-                result = (ulong)(((UInt128)result * baseValue) % modulus);
     }
 
-            remaining >>= 1;
+    private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
-            if (remaining == 0UL)
-            {
-                break;
-            }
-
-            baseValue = (ulong)(((UInt128)baseValue * baseValue) % modulus);
-        }
-
-        return result;
-    }
-
-    private static Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, GpuUInt128, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator)
     {
         return Pow2ModKernelCache.GetOrAdd(accelerator, static accel =>
         {
-            var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<ulong, Stride1D.Dense>, GpuUInt128, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernel);
+            var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernel);
             var kernel = KernelUtil.GetKernel(loaded);
-            return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, GpuUInt128, ArrayView1D<ulong, Stride1D.Dense>>>();
+            return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>>>();
         });
     }
 
-    private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, GpuUInt128 modulus, ArrayView1D<ulong, Stride1D.Dense> remainders)
+    private static void Pow2ModKernel(Index1D index, ArrayView1D<ulong, Stride1D.Dense> exponents, MontgomeryDivisorData divisor, ArrayView1D<ulong, Stride1D.Dense> remainders)
     {
         ulong exponent = exponents[index];
-        GpuUInt128 result = GpuUInt128.Pow2Mod(exponent, modulus);
+        remainders[index] = exponent.Pow2MontgomeryModWindowed(divisor, keepMontgomery: false);
-        remainders[index] = result.Low;
     }
 
     internal readonly record struct PrimeOrderGpuCapability(int ModulusBits, int ExponentBits)

todo.txt

@@ -1,2 +1 @@
-- Swap PrimeOrderGpuHeuristics.Pow2ModCpu for the existing Pow2MontgomeryModWindowed helper and plan the same windowed ladder for the GPU kernel implementation.
 - Extend the GPU numeric type probe to cover proposed custom types such as GpuUInt256 and GpuUInt512 once their definitions are available so we can plan the wide-order heuristics on the device.