From b07bf769c9feca0f0881c37f63cd87001f30f86c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 04:58:20 +0200 Subject: [PATCH] Switch single-use GPU exponent checks to standard modular math --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 17 +++++----- ...MersenneNumberDivisorByDivisorGpuTester.cs | 8 ++--- ...ersenneNumberDivisorResidueGpuEvaluator.cs | 9 +++++- .../Gpu/PrimeOrderGpuHeuristics.cs | 4 +-- PerfectNumbers.Core/ULongExtensions.cs | 31 +++++++++++++++++++ gpu-montgomery-usage.md | 18 +++++++++++ 6 files changed, 71 insertions(+), 16 deletions(-) create mode 100644 gpu-montgomery-usage.md diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 7269b31..61068f7 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/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 residueSpan = _montgomeryDeltas!.AsSpan(0, length); + Span 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.Shared.Rent(newCapacity); _cycleRemainders = ArrayPool.Shared.Rent(newCapacity); - _montgomeryDeltas = ArrayPool.Shared.Rent(newCapacity); + _residues = ArrayPool.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.Shared.Return(_montgomeryDeltas, clearArray: false); - _montgomeryDeltas = null; + ArrayPool.Shared.Return(_residues, clearArray: false); + _residues = null; } _bufferCapacity = 0; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 21769b7..bd96720 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -35,7 +35,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(CheckKernel)); private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeMontgomeryExponentKernel)); + _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeExponentKernel)); private Action, ArrayView, byte> GetRemainderDeltaKernel(Accelerator accelerator) => _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, 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 exponents, ArrayView results) + private static void ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView 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 diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs index 8fdd163..8e8879a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs +++ b/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); } } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index b403218..4df1437 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/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) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index fca8662..3388983 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/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) { diff --git a/gpu-montgomery-usage.md b/gpu-montgomery-usage.md new file mode 100644 index 0000000..c040d5c --- /dev/null +++ b/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】