diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 61068f7..7269b31 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[]? _residues; + private ulong[]? _montgomeryDeltas; private int _bufferCapacity; private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator; @@ -765,12 +765,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); EnsureGpuEvaluator(); - Span residueSpan = _residues!.AsSpan(0, length); + Span residueSpan = _montgomeryDeltas!.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] == 1UL ? (byte)1 : (byte)0; + hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0; } } @@ -800,7 +801,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } - if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null) + if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null) { return; } @@ -820,7 +821,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _primeDeltas = ArrayPool.Shared.Rent(newCapacity); _cycleRemainders = ArrayPool.Shared.Rent(newCapacity); - _residues = ArrayPool.Shared.Rent(newCapacity); + _montgomeryDeltas = ArrayPool.Shared.Rent(newCapacity); _bufferCapacity = newCapacity; } @@ -838,10 +839,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _cycleRemainders = null; } - if (_residues is not null) + if (_montgomeryDeltas is not null) { - ArrayPool.Shared.Return(_residues, clearArray: false); - _residues = null; + ArrayPool.Shared.Return(_montgomeryDeltas, clearArray: false); + _montgomeryDeltas = null; } _bufferCapacity = 0; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index bd96720..21769b7 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>(ComputeExponentKernel)); + _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeMontgomeryExponentKernel)); 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.Pow2ModBinaryGpu(modulus) == 1UL ? (byte)1 : (byte)0; + hits[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 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 ComputeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) + private static void ComputeMontgomeryExponentKernel(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.Pow2ModBinaryGpu(modulus); + results[index] = exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: true); } public ulong DivisorLimit diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs index 8e8879a..8fdd163 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs @@ -109,13 +109,6 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable } ulong exponent = exponents[idx]; - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[idx] = 0UL; - return; - } - - results[idx] = exponent.Pow2ModBinaryGpu(modulus); + results[idx] = exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: true); } } diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 4df1437..b403218 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.Pow2ModBinaryGpu(divisor.Modulus) == 1UL) + if (reduced.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL) { return false; } @@ -1595,7 +1595,7 @@ internal static partial class PrimeOrderGpuHeuristics private static bool Pow2EqualsOneKernel(ulong exponent, in MontgomeryDivisorData divisor) { - return exponent.Pow2ModBinaryGpu(divisor.Modulus) == 1UL; + return exponent.Pow2MontgomeryModWindowedGpu(divisor, keepMontgomery: false) == 1UL; } private static ulong CalculateByDoublingKernel(ulong prime) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 3388983..fca8662 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -310,37 +310,6 @@ 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 deleted file mode 100644 index c040d5c..0000000 --- a/gpu-montgomery-usage.md +++ /dev/null @@ -1,18 +0,0 @@ -# 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】