Switch single-use GPU exponent checks to standard modular math

PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs

@@ -652,7 +652,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
         private bool _disposed;
         private ulong[]? _primeDeltas;
         private ulong[]? _cycleRemainders;
-        private ulong[]? _montgomeryDeltas;
+        private ulong[]? _residues;
         private int _bufferCapacity;
         private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator;
 
@@ -765,13 +765,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
             ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan);
 
             EnsureGpuEvaluator();
-            Span<ulong> residueSpan = _montgomeryDeltas!.AsSpan(0, length);
+            Span<ulong> residueSpan = _residues!.AsSpan(0, length);
             _gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan);
 
-            ulong montgomeryOne = divisorData.MontgomeryOne;
             for (int i = 0; i < length; i++)
             {
-                hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0;
+                hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == 1UL ? (byte)1 : (byte)0;
             }
         }
 
@@ -801,7 +800,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 return;
             }
 
-            if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null)
+            if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null)
             {
                 return;
             }
@@ -821,7 +820,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
 
             _primeDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
             _cycleRemainders = ArrayPool<ulong>.Shared.Rent(newCapacity);
-            _montgomeryDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
+            _residues = ArrayPool<ulong>.Shared.Rent(newCapacity);
             _bufferCapacity = newCapacity;
         }
 
@@ -839,10 +838,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
                 _cycleRemainders = null;
             }
 
-            if (_montgomeryDeltas is not null)
+            if (_residues is not null)
             {
-                ArrayPool<ulong>.Shared.Return(_montgomeryDeltas, clearArray: false);
-                _montgomeryDeltas = null;
+                ArrayPool<ulong>.Shared.Return(_residues, clearArray: false);
+                _residues = null;
             }
 
             _bufferCapacity = 0;

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -35,7 +35,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<MontgomeryDivisorData>, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
 
     private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeExponent(Accelerator accelerator) =>
-        _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeMontgomeryExponentKernel));
+        _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeExponentKernel));
 
     private Action<Index1D, ArrayView<ulong>, ArrayView<byte>, byte> GetRemainderDeltaKernel(Accelerator accelerator) =>
         _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<byte>, byte>(ComputeRemainderDeltasKernel));
@@ -569,7 +569,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulong exponent = exponents[index];
-        hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL ? (byte)1 : (byte)0;
+        hits[index] = exponent.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0;
     }
 
     public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession
@@ -777,7 +777,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         mask[globalIndex] = 1;
     }
 
-    private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
+    private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
     {
         ulong modulus = divisor.Modulus;
         if (modulus <= 1UL || (modulus & 1UL) == 0UL)
@@ -788,7 +788,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         ulong exponent = exponents[index];
 
-        results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true);
+        results[index] = exponent.Pow2ModBinaryGpu(modulus);
     }
 
     public ulong DivisorLimit

PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs

@@ -109,6 +109,13 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable
         }
 
         ulong exponent = exponents[idx];
-        results[idx] = exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: true);
+        ulong modulus = divisorData.Modulus;
+        if (modulus <= 1UL || (modulus & 1UL) == 0UL)
+        {
+            results[idx] = 0UL;
+            return;
+        }
+
+        results[idx] = exponent.Pow2ModBinaryGpu(modulus);
     }
 }

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -811,7 +811,7 @@ internal static partial class PrimeOrderGpuHeuristics
             }
 
             ulong reduced = phi / factor;
-            if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL)
+            if (reduced.Pow2ModBinaryGpu(divisor.Modulus) == 1UL)
             {
                 return false;
             }
@@ -1595,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics
 
     private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor)
     {
-        return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL;
+        return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL;
     }
 
     private static ulong CalculateByDoublingKernel(ulong prime)

PerfectNumbers.Core/ULongExtensions.cs

@@ -310,6 +310,37 @@ public static class ULongExtensions
         return state.MulModWithNativeModulo(b, modulus);
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong Pow2ModBinaryGpu(this ulong exponent, ulong modulus)
+    {
+        if (modulus <= 1UL)
+        {
+            return 0UL;
+        }
+
+        ulong remainingExponent = exponent;
+        ulong result = 1UL % modulus;
+        ulong baseValue = 2UL % modulus;
+
+        while (remainingExponent != 0UL)
+        {
+            if ((remainingExponent & 1UL) != 0UL)
+            {
+                result = result.MulMod64GpuCompatible(baseValue, modulus);
+            }
+
+            remainingExponent >>= 1;
+            if (remainingExponent == 0UL)
+            {
+                break;
+            }
+
+            baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus);
+        }
+
+        return result;
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public static ulong MontgomeryMultiply(this ulong a, ulong b, ulong modulus, ulong nPrime)
     {

gpu-montgomery-usage.md

@@ -0,0 +1,18 @@
+# GPU Modular Arithmetic in the `--mersenne=bydivisor` Path
+
+The GPU implementation now routes single-use exponent checks through standard modular arithmetic while keeping Montgomery reduction only where repeated multiplications reuse the same modulus.
+
+## Kernels using standard modular arithmetic
+
+- `CheckKernel` in `MersenneNumberDivisorByDivisorGpuTester` evaluates each divisor candidate with `Pow2ModBinaryGpu` and compares the result directly against 1. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L561-L645】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
+- `ComputeExponentKernel` in the same tester precomputes batched residues with the plain helper so CPU lookups receive standard-domain results. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs†L780-L791】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
+- `ComputeResiduesKernel` inside `MersenneNumberDivisorResidueGpuEvaluator` now mirrors that logic when building residue tables for mixed CPU/GPU scans. 【F:PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs†L102-L114】【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】
+- The order-search helpers in `PrimeOrderGpuHeuristics` (`TrySpecialMaxKernel`, `ExponentLoweringKernel`, and `Pow2EqualsOneKernel`) also rely on `Pow2ModBinaryGpu` for their equality checks. 【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L798-L866】【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L1596-L1598】
+
+## Kernels still using Montgomery arithmetic
+
+- `MontgomeryOddPowerKernelScan` in `GpuKernelPool` continues to build windowed odd-power tables with `MontgomeryMultiply`, which amortizes the Montgomery setup over many multiplications per modulus. 【F:PerfectNumbers.Core/Gpu/GpuKernelPool.cs†L200-L237】【F:PerfectNumbers.Core/ULongExtensions.cs†L532-L559】
+
+## Shared helper
+
+- `Pow2ModBinaryGpu` in `ULongExtensions` implements the GPU-friendly binary exponentiation loop that powers the standard modular path across these kernels. 【F:PerfectNumbers.Core/ULongExtensions.cs†L314-L339】