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..2de1d6e 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) @@ -1810,58 +1810,15 @@ internal static partial class PrimeOrderGpuHeuristics private static void Pow2ModKernel(Index1D index, ArrayView1D exponents, MontgomeryDivisorData divisor, ArrayView1D remainders) { - ulong exponent = exponents[index]; ulong modulus = divisor.Modulus; - // ulong nPrime = divisor.NPrime; - // ulong result = divisor.MontgomeryOne; - // ulong baseValue = divisor.MontgomeryTwo; - // ulong remainingExponent = exponent; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + remainders[index] = 0UL; + return; + } - // if (remainingExponent == 0UL) - // { - // remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime); - // return; - // } - - // while (remainingExponent != 0UL) - // { - // if ((remainingExponent & 1UL) != 0UL) - // { - // result = result.MontgomeryMultiply(baseValue, modulus, nPrime); - // } - - // remainingExponent >>= 1; - // if (remainingExponent == 0UL) - // { - // break; - // } - - // // Reusing baseValue to hold base^2 so we avoid allocating another register for the squared term. - // baseValue = baseValue.MontgomeryMultiply(baseValue, modulus, nPrime); - // } - - // remainders[index] = result.MontgomeryMultiply(1UL, modulus, nPrime); - - // Plain modular exponentiation without Montgomery reduction kept for reference. - ulong plainResult = 1UL;// % modulus; - ulong plainBase = 2UL;// % modulus; - while (exponent != 0UL) - { - if ((exponent & 1UL) != 0UL) - { - plainResult = (plainResult * plainBase) % modulus; - } - - exponent >>= 1; - if (exponent == 0UL) - { - break; - } - - plainBase = (plainBase * plainBase) % modulus; - } - - remainders[index] = plainResult; + ulong exponent = exponents[index]; + remainders[index] = exponent.Pow2ModBinaryGpu(modulus); } private static Action, GpuUInt128, ArrayView1D> GetPow2ModWideKernel(Accelerator accelerator) 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..5589ccc --- /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`, `Pow2EqualsOneKernel`, and the narrow `Pow2ModKernel`) also rely on `Pow2ModBinaryGpu` for their equality checks and residue generation. 【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L798-L866】【F:PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs†L1812-L1823】 + +## 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】 diff --git a/todo_list.txt b/todo_list.txt index 3ed48dc..9a7cdab 100644 --- a/todo_list.txt +++ b/todo_list.txt @@ -122,7 +122,6 @@ PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:203: PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207: // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214: // TODO: Remove this no-cycle branch and require on-demand cycle computation (without PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240: // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the -PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262: : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283: // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356: // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396: // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here @@ -266,7 +265,6 @@ todo_list.txt:121:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTeste todo_list.txt:122:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:207: // TODO: When divisorCycle is zero (cache miss), compute only that single cycle on the device todo_list.txt:123:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:214: // TODO: Remove this no-cycle branch and require on-demand cycle computation (without todo_list.txt:124:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:240: // TODO: Replace this linear increment with the batched divisor-cycle walker validated in the -todo_list.txt:125:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:262: : prime.Pow2MontgomeryMod(divisorData); // TODO: Switch to the upcoming ProcessEightBitWindows helper once the todo_list.txt:126:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:283: // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI todo_list.txt:127:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:356: // TODO: Replace this guard with an on-demand CPU/GPU cycle computation path that skips todo_list.txt:128:PerfectNumbers.Core\Cpu\MersenneNumberDivisorByDivisorCpuTester.cs:396: // TODO: Remove this no-cycle fallback once divisor cycles are mandatory; the per-prime pow2 ladder here @@ -407,17 +405,6 @@ PerfectNumbers.Core\ModResidueTracker.cs:257: // TODO: Us PerfectNumbers.Core\ModResidueTracker.cs:263: // TODO: Consult MersenneDivisorCycles.Shared (or similar caches) here so we reuse precomputed cycle lengths PerfectNumbers.Core\ModResidueTracker.cs:295: // TODO: Route this through the planned UInt128 windowed powmod so the inner MulMod adopts the PerfectNumbers.Core\ModResidueTracker.cs:314: // TODO: Replace this double-and-add fallback with the UInt128 intrinsic-backed multiplier once the -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:49: set => _useDivisorCycles = value; // TODO: Remove this flag once GPU runs unconditionally fetch cycle lengths so every divisor check benefits from the measured speedups. -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:298: // TODO: Stage MontgomeryDivisorData in a GPU-ready cache keyed by divisor so we can copy -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:372: // TODO: Replace the block lease with a direct single-cycle lookup so the GPU path -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:395: // TODO: Drop this block hand-off once the cache exposes a direct -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:409: // TODO: Replace this `%` with the cycle-cache helper that subtracts precomputed -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:425: // TODO: Trigger an immediate GPU computation for the single missing cycle without storing -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:476: // TODO: Preload ProcessEightBitWindows-friendly divisor metadata for GPU batches so this -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:563: hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; // TODO: Replace Pow2MontgomeryMod with the ProcessEightBitWindows helper once available so GPU by-divisor scans reuse the benchmarked fast pow2 ladder. -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:632: // TODO: Replace this guard with an on-demand cycle computation that runs on the device selected -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:787: // TODO: Replace this `%` with the divisor-cycle remainder helper so we reuse cached -PerfectNumbers.Core\Gpu\MersenneNumberDivisorByDivisorGpuTester.cs:858: results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); // TODO: Swap this Montgomery-domain pow2 with the windowed helper slated to replace the slower ladder in Pow2Minus1Mod benchmarks. PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:43: // TODO: Swap this per-run enumerator creation for the shared PrimeIterator so by-divisor scans reuse the PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:219: // TODO: Swap Task.Factory.StartNew for the shared low-overhead work queue once the scanner's PerfectNumbers.Core\MersenneNumberDivisorByDivisorTester.cs:240: bool useDivisorCycles = tester.UseDivisorCycles; // TODO: Remove the conditional path once all testers always enable divisor cycles to keep every divisor scan on the faster cached-length track.