diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs new file mode 100644 index 0000000..d4af20d --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedBenchmarkData.cs @@ -0,0 +1,94 @@ +using System; +using PerfectNumbers.Core; + +namespace EvenPerfectBitScanner.Benchmarks; + +internal sealed class Pow2MontgomeryModWindowedBenchmarkInputs +{ + public Pow2MontgomeryModWindowedBenchmarkInputs(Pow2MontgomeryModWindowedBenchmarkCase[] cases) + { + Cases = cases; + } + + public Pow2MontgomeryModWindowedBenchmarkCase[] Cases { get; } +} + +internal readonly struct Pow2MontgomeryModWindowedBenchmarkCase +{ + public Pow2MontgomeryModWindowedBenchmarkCase( + ulong exponent, + ulong modulus, + in MontgomeryDivisorData divisor, + ulong cycleLength, + ulong reducedExponent) + { + Exponent = exponent; + Modulus = modulus; + Divisor = divisor; + CycleLength = cycleLength; + ReducedExponent = reducedExponent; + } + + public ulong Exponent { get; } + + public ulong Modulus { get; } + + public MontgomeryDivisorData Divisor { get; } + + public ulong CycleLength { get; } + + public ulong ReducedExponent { get; } +} + +internal static class Pow2MontgomeryModWindowedBenchmarkInputsProvider +{ + private static readonly Lazy Cache = new(Create); + + public static Pow2MontgomeryModWindowedBenchmarkInputs Instance => Cache.Value; + + private static Pow2MontgomeryModWindowedBenchmarkInputs Create() + { + Pow2MontgomeryModWindowedBenchmarkSeed[] seeds = + { + new Pow2MontgomeryModWindowedBenchmarkSeed(138_000_001UL, 1UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(150_000_013UL, 2UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(175_000_019UL, 3UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(190_000_151UL, 4UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(210_000_089UL, 5UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(230_000_039UL, 6UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(260_000_111UL, 7UL), + new Pow2MontgomeryModWindowedBenchmarkSeed(300_000_007UL, 8UL) + }; + + var cases = new Pow2MontgomeryModWindowedBenchmarkCase[seeds.Length]; + for (int i = 0; i < seeds.Length; i++) + { + Pow2MontgomeryModWindowedBenchmarkSeed seed = seeds[i]; + ulong modulus = 2UL * seed.Multiplier * seed.Exponent + 1UL; + MontgomeryDivisorData divisor = MontgomeryDivisorData.FromModulus(modulus); + ulong cycleLength = MersenneDivisorCycles.CalculateCycleLength(modulus, divisor); + ulong reducedExponent = seed.Exponent % cycleLength; + cases[i] = new Pow2MontgomeryModWindowedBenchmarkCase( + seed.Exponent, + modulus, + divisor, + cycleLength, + reducedExponent); + } + + return new Pow2MontgomeryModWindowedBenchmarkInputs(cases); + } + + private readonly struct Pow2MontgomeryModWindowedBenchmarkSeed + { + public Pow2MontgomeryModWindowedBenchmarkSeed(ulong exponent, ulong multiplier) + { + Exponent = exponent; + Multiplier = multiplier; + } + + public ulong Exponent { get; } + + public ulong Multiplier { get; } + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks.cs new file mode 100644 index 0000000..99ad30a --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks.cs @@ -0,0 +1,94 @@ +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core; + +namespace EvenPerfectBitScanner.Benchmarks; + +/// +/// Benchmarks CPU Montgomery exponentiation helpers for exponents at or above 138 million. +/// Reuses the shared benchmark input cache so every method processes identical samples. +/// +[MemoryDiagnoser] +[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] +public class Pow2MontgomeryModWindowedKeepMontgomeryCpuBenchmarks +{ + private Pow2MontgomeryModWindowedBenchmarkInputs _inputs = null!; + + [GlobalSetup] + public void Setup() + { + _inputs = Pow2MontgomeryModWindowedBenchmarkInputsProvider.Instance; + } + + /// + /// Baseline CPU windowed Montgomery exponentiation that converts the residue back to standard form. + /// + [Benchmark(Baseline = true)] + public ulong WindowedMontgomeryConvert() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedCpu(current.Divisor, keepMontgomery: false); + } + + return checksum; + } + + /// + /// CPU windowed Montgomery exponentiation that keeps the residue in the Montgomery domain. + /// + [Benchmark] + public ulong WindowedMontgomeryKeep() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedCpu(current.Divisor, keepMontgomery: true); + } + + return checksum; + } + + /// + /// Uses the known cycle length to fold the exponent before running the windowed Montgomery ladder. + /// + [Benchmark] + public ulong WindowedMontgomeryWithCycle() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWithCycleCpu(current.CycleLength, current.Divisor); + } + + return checksum; + } + + /// + /// Runs the windowed Montgomery ladder on the reduced exponent that the divisor cycle exposes. + /// + [Benchmark] + public ulong WindowedMontgomeryFromCycleRemainder() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.ReducedExponent.Pow2MontgomeryModFromCycleRemainderCpu(current.Divisor); + } + + return checksum; + } +} diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs new file mode 100644 index 0000000..781c33b --- /dev/null +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks.cs @@ -0,0 +1,78 @@ +using System; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using PerfectNumbers.Core; + +namespace EvenPerfectBitScanner.Benchmarks; + +/// +/// Benchmarks GPU Montgomery exponentiation helpers for exponents at or above 138 million. +/// The sample set is generated once and reused across all benchmark methods so every variant +/// processes the same data without incurring repeated initialization costs. +/// +[MemoryDiagnoser] +[SimpleJob(RuntimeMoniker.Net80, launchCount: 1, warmupCount: 1, iterationCount: 5)] +public class Pow2MontgomeryModWindowedKeepMontgomeryGpuBenchmarks +{ + private Pow2MontgomeryModWindowedBenchmarkInputs _inputs = null!; + + [GlobalSetup] + public void Setup() + { + _inputs = Pow2MontgomeryModWindowedBenchmarkInputsProvider.Instance; + } + + /// + /// Baseline GPU windowed Montgomery exponentiation that converts back to the standard residue. + /// + [Benchmark(Baseline = true)] + public ulong WindowedMontgomeryConvert() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedGpu(current.Divisor, keepMontgomery: false); + } + + return checksum; + } + + /// + /// GPU windowed Montgomery exponentiation that keeps the result in the Montgomery domain. + /// + [Benchmark] + public ulong WindowedMontgomeryKeep() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2MontgomeryModWindowedGpu(current.Divisor, keepMontgomery: true); + } + + return checksum; + } + + /// + /// GPU windowed modular exponentiation without Montgomery reduction for comparison. + /// + [Benchmark] + public ulong WindowedStandardMod() + { + ulong checksum = 0UL; + Pow2MontgomeryModWindowedBenchmarkCase[] cases = _inputs.Cases; + + for (int i = 0; i < cases.Length; i++) + { + ref readonly Pow2MontgomeryModWindowedBenchmarkCase current = ref cases[i]; + checksum ^= current.Exponent.Pow2ModWindowedGpu(current.Modulus); + } + + return checksum; + } +} diff --git a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs index 1ef8592..469ce8e 100644 --- a/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/GpuKernelCompilationTests.cs @@ -75,16 +75,16 @@ public class GpuKernelCompilationTests ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetKernel not found."); private static readonly MethodInfo DivisorByDivisorGetExponentKernelMethod = - typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetKernelByPrimeExponent", BindingFlags.NonPublic | BindingFlags.Instance) + typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetStepperKernel", BindingFlags.NonPublic | BindingFlags.Instance) ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetKernelByPrimeExponent not found."); - private static readonly MethodInfo DivisorByDivisorGetRemainderDeltaKernelMethod = - typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetRemainderDeltaKernel", BindingFlags.NonPublic | BindingFlags.Instance) - ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetRemainderDeltaKernel not found."); + private static readonly MethodInfo DivisorByDivisorGetCandidateMaskKernelMethod = + typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetCandidateMaskKernel", BindingFlags.NonPublic | BindingFlags.Instance) + ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetCandidateMaskKernel not found."); - private static readonly MethodInfo DivisorByDivisorGetRemainderScanKernelMethod = - typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetRemainderScanKernel", BindingFlags.NonPublic | BindingFlags.Instance) - ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetRemainderScanKernel not found."); + private static readonly MethodInfo DivisorByDivisorGetCandidateGenerationKernelMethod = + typeof(MersenneNumberDivisorByDivisorGpuTester).GetMethod("GetCandidateGenerationKernel", BindingFlags.NonPublic | BindingFlags.Instance) + ?? throw new InvalidOperationException("MersenneNumberDivisorByDivisorGpuTester.GetCandidateGenerationKernel not found."); private static readonly MethodInfo DivisorTesterGetKernelMethod = typeof(MersenneNumberDivisorGpuTester).GetMethod("GetKernel", BindingFlags.NonPublic | BindingFlags.Instance) @@ -157,9 +157,9 @@ public class GpuKernelCompilationTests yield return Loader("NttGpuMath.ForwardKernel", CompileNttForwardKernel); yield return Loader("NttGpuMath.InverseKernel", CompileNttInverseKernel); yield return Loader("MersenneNumberDivisorByDivisorGpuTester.CheckKernel", CompileDivisorByDivisorCheckKernel); - yield return Loader("MersenneNumberDivisorByDivisorGpuTester.ComputeMontgomeryExponentKernel", CompileDivisorByDivisorExponentKernel); - yield return Loader("MersenneNumberDivisorByDivisorGpuTester.ComputeRemainderDeltasKernel", CompileDivisorByDivisorRemainderDeltaKernel); - yield return Loader("MersenneNumberDivisorByDivisorGpuTester.AccumulateRemaindersKernel", CompileDivisorByDivisorRemainderScanKernel); + yield return Loader("MersenneNumberDivisorByDivisorGpuTester.EvaluateDivisorWithStepperKernel", CompileDivisorByDivisorExponentKernel); + yield return Loader("MersenneNumberDivisorByDivisorGpuTester.EvaluateCandidateMaskKernel", CompileDivisorByDivisorCandidateMaskKernel); + yield return Loader("MersenneNumberDivisorByDivisorGpuTester.GenerateCandidatesKernel", CompileDivisorByDivisorCandidateGenerationKernel); yield return Loader("MersenneNumberDivisorGpuTester.Kernel", CompileDivisorTesterKernel); yield return Loader("MersenneNumberLucasLehmerGpuTester.Kernel", CompileLucasLehmerKernel); yield return Loader("MersenneNumberLucasLehmerGpuTester.AddSmallKernel", CompileLucasLehmerAddSmallKernel); @@ -287,16 +287,16 @@ public class GpuKernelCompilationTests _ = DivisorByDivisorGetExponentKernelMethod.Invoke(tester, new object[] { accelerator }); } - private static void CompileDivisorByDivisorRemainderDeltaKernel(Accelerator accelerator) + private static void CompileDivisorByDivisorCandidateMaskKernel(Accelerator accelerator) { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - _ = DivisorByDivisorGetRemainderDeltaKernelMethod.Invoke(tester, new object[] { accelerator }); + _ = DivisorByDivisorGetCandidateMaskKernelMethod.Invoke(tester, new object[] { accelerator }); } - private static void CompileDivisorByDivisorRemainderScanKernel(Accelerator accelerator) + private static void CompileDivisorByDivisorCandidateGenerationKernel(Accelerator accelerator) { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - _ = DivisorByDivisorGetRemainderScanKernelMethod.Invoke(tester, new object[] { accelerator }); + _ = DivisorByDivisorGetCandidateGenerationKernelMethod.Invoke(tester, new object[] { accelerator }); } private static void CompileDivisorTesterKernel(Accelerator accelerator) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 8f67d06..a382431 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -85,11 +85,17 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_tester_requires_configuration() + public void ByDivisor_tester_defaults_to_zero_limit_before_configuration() { var tester = new MersenneNumberDivisorByDivisorGpuTester(); - Action act = () => tester.IsPrime(31UL, out _); - act.Should().Throw(); + + tester.DivisorLimit.Should().Be(0UL); + + ulong[] primes = { 31UL, 37UL }; + ulong[] allowed = new ulong[primes.Length]; + + tester.PrepareCandidates(primes, allowed); + allowed.Should().OnlyContain(value => value == 0UL); } [Fact] diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepper.cs index 0722880..9d6ad3f 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepper.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepper.cs @@ -121,6 +121,30 @@ internal struct ExponentRemainderStepper return _currentMontgomery == _montgomeryOne; } + /// + /// Advances the GPU stepping state to the next exponent in a strictly ascending sequence that reuses the same divisor. + /// Callers must invoke beforehand. + /// + /// Next exponent in the ascending sequence. + /// The canonical residue produced by 2^exponent mod divisor. + public ulong ComputeNextGpu(ulong exponent) + { + AdvanceGpu(exponent); + return ReduceCurrent(); + } + + /// + /// Advances the GPU stepping state to the next exponent and reports whether it evaluates to unity. + /// Callers must seed the stepper through or . + /// + /// Next exponent in the ascending sequence. + /// when the canonical residue equals one. + public bool ComputeNextGpuIsUnity(ulong exponent) + { + AdvanceGpu(exponent); + return _currentMontgomery == _montgomeryOne; + } + /// /// Seeds the stepper with a precomputed Montgomery residue produced by the caller. /// The provided exponent becomes the new baseline and all subsequent calls must continue with greater exponents. @@ -145,6 +169,14 @@ internal struct ExponentRemainderStepper /// Montgomery residue delta corresponding to the exponent gap. /// Outputs whether the resulting residue equals Montgomery one. /// when the advance succeeds. + private void AdvanceGpu(ulong exponent) + { + ulong delta = exponent - _previousExponent; + ulong multiplier = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, delta); + _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); + _previousExponent = exponent; + } + public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity) { // TODO: Once the divisor-cycle cache exposes a direct Montgomery delta, multiply it here instead @@ -164,7 +196,7 @@ internal struct ExponentRemainderStepper private void InitializeGpuState(ulong exponent) { - _currentMontgomery = exponent.Pow2MontgomeryModWindowedCpu(_divisor, keepMontgomery: true); + _currentMontgomery = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(_divisor, exponent); _previousExponent = exponent; } diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index eece679..d2ec9bf 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -2,6 +2,7 @@ using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Gpu; @@ -18,10 +19,64 @@ internal static class DivisorByDivisorKernels hits[index] = montgomeryResult == 1UL ? (byte)1 : (byte)0; } - public static void GenerateCandidatesAndGapsKernel( + public static void EvaluateDivisorWithStepperKernel( + Index1D _, + MontgomeryDivisorData divisor, + ulong divisorCycle, + ulong firstCycleRemainder, + ArrayView exponents, + ArrayView hits) + { + int length = (int)exponents.Length; + // EvenPerfectBitScanner always supplies at least one exponent per GPU divisor check, + // so the guard stays commented out to keep the kernel branch-free. + // if (length == 0) + // { + // return; + // } + + var exponentStepper = new ExponentRemainderStepper(divisor); + ulong firstPrime = exponents[0]; + bool firstUnity = exponentStepper.InitializeGpuIsUnity(firstPrime); + // EvenPerfectBitScanner only schedules divisors with non-zero cycle lengths (equal to the tested prime), + // so the no-cycle path remains commented out to keep the kernel branch-free. + // if (divisorCycle == 0UL) + // { + // hits[0] = firstUnity ? (byte)1 : (byte)0; + // for (int i = 1; i < length; i++) + // { + // ulong primeNoCycle = exponents[i]; + // bool unityNoCycle = exponentStepper.ComputeNextGpuIsUnity(primeNoCycle); + // hits[i] = unityNoCycle ? (byte)1 : (byte)0; + // } + // return; + // } + + ulong cycleLength = divisorCycle; + ulong previousPrime = firstPrime; + // The host precomputes the initial cycle remainder to keep the kernel free of modulo operations. + ulong cycleRemainder = firstCycleRemainder; + hits[0] = cycleRemainder == 0UL && firstUnity ? (byte)1 : (byte)0; + + for (int i = 1; i < length; i++) + { + ulong prime = exponents[i]; + ulong delta = prime - previousPrime; + previousPrime = prime; + + var remainderValue = new GpuUInt128(cycleRemainder); + remainderValue.AddMod(delta, cycleLength); + cycleRemainder = remainderValue.Low; + + hits[i] = cycleRemainder == 0UL + ? (exponentStepper.ComputeNextGpuIsUnity(prime) ? (byte)1 : (byte)0) + : (byte)0; + } + } + + public static void GenerateCandidatesKernel( Index1D index, ArrayView candidates, - ArrayView gaps, ArrayView nextStarts, GpuUInt128 startValue, GpuUInt128 stride) @@ -41,23 +96,6 @@ internal static class DivisorByDivisorKernels candidateValue.Add(startValue); candidates[globalIndex] = candidateValue.Low; - // The scanner supplies gap buffers that mirror the candidate view length, - // so the additional bounds guard remains commented out for the same reason as above. - // int gapLength = (int)gaps.Length; - // if (globalIndex >= gapLength) - // { - // return; - // } - - ulong previousOffset = offset > 0UL ? offset - 1UL : 0UL; - GpuUInt128 previousValue = stride; - previousValue.Mul(previousOffset); - previousValue.Add(startValue); - - GpuUInt128 gapValue = candidateValue; - gapValue.Sub(previousValue); - gaps[globalIndex] = gapValue.Low; - if (globalIndex == 0) { // EvenPerfectBitScanner provides a single-element nextStarts view per launch, so the guard stays commented out. @@ -75,68 +113,11 @@ internal static class DivisorByDivisorKernels } } - public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) - { - deltas[index] = (byte)(gaps[index] % modulus); - } - - public static void AccumulateRemaindersKernel( - Index1D index, - ArrayView deltas, - ArrayView remainders, - byte baseRemainder, - byte modulus) - { - int globalIndex = index; - // EvenPerfectBitScanner launches this stream kernel with an extent matching the remainder span length, - // so the defensive guard stays commented out to keep the kernel branch-free. - // int length = (int)remainders.Length; - // if (globalIndex >= length) - // { - // return; - // } - - int localIndex = Group.IdxX; - int groupSize = Group.Dimension.X; - var shared = SharedMemory.GetDynamic(); - - shared[localIndex] = deltas[globalIndex]; - Group.Barrier(); - - int offset = 1; - int mod = modulus; - while (offset < groupSize) - { - bool apply = localIndex >= offset; - int addend = apply ? shared[localIndex - offset] : 0; - - Group.Barrier(); - - int current = shared[localIndex]; - int sum = apply ? current + addend : current; - int wrap = sum >= mod ? 1 : 0; - sum -= wrap * mod; - shared[localIndex] = apply ? sum : current; - - Group.Barrier(); - - offset <<= 1; - } - - int remainder = baseRemainder + shared[localIndex]; - int remainderWrap = remainder >= mod ? 1 : 0; - remainder -= remainderWrap * mod; - remainders[globalIndex] = (byte)remainder; - } - public static void EvaluateCandidateMaskKernel( Index1D index, - ArrayView remainder10, - ArrayView remainder8, - ArrayView remainder5, - ArrayView remainder3, - byte lastIsSevenFlag, - ArrayView mask) + ArrayView candidates, + ArrayView mask, + byte lastIsSevenFlag) { int globalIndex = index; // EvenPerfectBitScanner launches this auto-grouped kernel with an extent matching the mask length, @@ -147,15 +128,18 @@ internal static class DivisorByDivisorKernels // return; // } - byte value10 = remainder10[globalIndex]; + ulong candidate = candidates[globalIndex]; + byte remainder10 = (byte)(candidate % 10UL); bool lastIsSeven = lastIsSevenFlag != 0; - bool acceptSeven = value10 == 3 || value10 == 7 || value10 == 9; - bool acceptNonSeven = value10 == 1 || value10 == 3 || value10 == 7 || value10 == 9; + bool acceptSeven = remainder10 == 3 || remainder10 == 7 || remainder10 == 9; + bool acceptNonSeven = remainder10 == 1 || remainder10 == 3 || remainder10 == 7 || remainder10 == 9; bool accept10 = lastIsSeven ? acceptSeven : acceptNonSeven; - byte value8 = remainder8[globalIndex]; - bool accept8 = value8 == 1 || value8 == 7; - bool accept5 = remainder5[globalIndex] != 0; - bool accept3 = remainder3[globalIndex] != 0; + byte remainder8 = (byte)(candidate & 7UL); + bool accept8 = remainder8 == 1 || remainder8 == 7; + byte remainder5 = (byte)(candidate % 5UL); + bool accept5 = remainder5 != 0; + byte remainder3 = (byte)(candidate % 3UL); + bool accept3 = remainder3 != 0; byte accepted = (byte)((accept10 && accept8 && accept3 && accept5) ? 1 : 0); mask[globalIndex] = accepted; } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 33ed320..4305f0e 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -1,3 +1,4 @@ +using System; using System.Buffers; using System.Collections.Concurrent; using System.Runtime.CompilerServices; @@ -11,11 +12,10 @@ namespace PerfectNumbers.Core.Gpu; /// /// Scans Mersenne divisors on the GPU for prime exponents p >= 31 using per-divisor Montgomery data. +/// Consumers must call before invoking other members so the divisor limits are populated. /// public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDivisorByDivisorTester { - private const int RemainderSlotCount = 4; - private int _gpuBatchSize = GpuConstants.ScanBatchSize; // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, // so the synchronization fields from the previous implementation remain commented out here. @@ -24,14 +24,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // private bool _isConfigured; private readonly ConcurrentDictionary, ArrayView, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); - private readonly ConcurrentDictionary, ArrayView, byte>> _remainderDeltaKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView, byte, byte>> _remainderScanKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView, byte, ArrayView>> _candidateMaskKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView>> _stepperKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView, byte>> _candidateMaskKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView>> _candidateCompactionKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView>> _hitDetectionKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView>> _hitConversionKernelCache = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentBag _sessionPool = new(); @@ -39,21 +36,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private Action, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(DivisorByDivisorKernels.CheckKernel)); - private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.ComputeMontgomeryExponentKernel)); + private Action, ArrayView> GetStepperKernel(Accelerator accelerator) => + _stepperKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.EvaluateDivisorWithStepperKernel)); - private Action, ArrayView, byte> GetRemainderDeltaKernel(Accelerator accelerator) => - _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, byte>(DivisorByDivisorKernels.ComputeRemainderDeltasKernel)); - - private Action, ArrayView, byte, byte> GetRemainderScanKernel(Accelerator accelerator) => - _remainderScanKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadStreamKernel, ArrayView, byte, byte>(DivisorByDivisorKernels.AccumulateRemaindersKernel)); - - private Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> GetCandidateMaskKernel(Accelerator accelerator) => + private Action, ArrayView, byte> GetCandidateMaskKernel(Accelerator accelerator) => _candidateMaskKernelCache.GetOrAdd( accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, byte, ArrayView>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); + acc => acc.LoadAutoGroupedStreamKernel, ArrayView, byte>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); private Action, ArrayView, ArrayView, ArrayView> GetCandidateCompactionKernel(Accelerator accelerator) => _candidateCompactionKernelCache.GetOrAdd( @@ -65,15 +54,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.PublishFirstHitKernel)); - private Action, ArrayView> GetResultConversionKernel(Accelerator accelerator) => - _hitConversionKernelCache.GetOrAdd( - accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.ConvertMontgomeryResultsToHitsKernel)); - - private Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => + private Action, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => _candidateGenerationKernelCache.GetOrAdd( accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, GpuUInt128, GpuUInt128>(DivisorByDivisorKernels.GenerateCandidatesAndGapsKernel)); + acc => acc.LoadAutoGroupedStreamKernel, ArrayView, GpuUInt128, GpuUInt128>(DivisorByDivisorKernels.GenerateCandidatesKernel)); public int GpuBatchSize { @@ -152,12 +136,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.HitIndexBuffer, resources.CandidateBuffer, resources.NextStartBuffer, - resources.DivisorDeltaBuffer, - resources.RemainderDeltaBuffer, - resources.RemainderBuffer, resources.FilterCountBuffer, - resources.RemainderDeltaKernel, - resources.RemainderScanKernel, resources.CandidateMaskKernel, resources.CandidateCompactionKernel, resources.CandidateGenerationKernel, @@ -226,15 +205,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D hitIndexBuffer, MemoryBuffer1D candidateBuffer, MemoryBuffer1D nextStartBuffer, - MemoryBuffer1D divisorDeltaBuffer, - MemoryBuffer1D remainderDeltaBuffer, - MemoryBuffer1D remainderBuffer, MemoryBuffer1D filterCountBuffer, - Action, ArrayView, byte> remainderDeltaKernel, - Action, ArrayView, byte, byte> remainderScanKernel, - Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> candidateMaskKernel, + Action, ArrayView, byte> candidateMaskKernel, Action, ArrayView, ArrayView, ArrayView> candidateCompactionKernel, - Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, + Action, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, Action, ArrayView> hitDetectionKernel, in ulong[] divisors, ulong[] exponents, @@ -276,31 +250,20 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; ulong currentK = 1UL; - byte step10 = (byte)(twoP128 % 10UL); - byte step8 = (byte)(twoP128 % 8UL); - byte step5 = (byte)(twoP128 % 5UL); - byte step3 = (byte)(twoP128 % 3UL); - UInt128 currentDivisor128 = firstDivisor128; - byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL)); - byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL)); - byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL)); - byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL)); bool lastIsSeven = (prime & 3UL) == 3UL; byte lastIsSevenFlag = lastIsSeven ? (byte)1 : (byte)0; DivisorCycleCache cycleCache = DivisorCycleCache.Shared; - int remainderStride = (int)divisorDeltaBuffer.Length; - remainderStride = remainderStride <= 0 ? chunkCapacity : remainderStride; ArrayPool ulongPool = ThreadStaticPools.UlongPool; ulong[] filteredDivisors = ulongPool.Rent(chunkCapacity); - var viewCache = new BufferSliceCache(divisorDeltaBuffer, remainderDeltaBuffer, remainderBuffer, hitsBuffer, remainderStride); Span candidateStorage = cycleCandidates.AsSpan(); Span filteredStorage = filteredDivisors.AsSpan(); Span cycleLengthStorage = cycleLengths.AsSpan(); Span divisorDataStorage = divisorData.AsSpan(); Span divisorStorage = divisors.AsSpan(); Span exponentStorage = exponents.AsSpan(); + ArrayView1D filteredCountView = filterCountBuffer.View; while (currentK <= maxK && !composite) { @@ -317,29 +280,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Span candidateSpan = candidateStorage[..effectiveCount]; - viewCache.Ensure(effectiveCount); - - ArrayView1D gapView = viewCache.DivisorDelta; - ArrayView1D deltaView = viewCache.RemainderDelta; - ArrayView1D remainder10View = viewCache.Remainder10; - ArrayView1D remainder8View = viewCache.Remainder8; - ArrayView1D remainder5View = viewCache.Remainder5; - ArrayView1D remainder3View = viewCache.Remainder3; - ArrayView1D maskView = viewCache.Hits; - ArrayView1D filteredCountView = filterCountBuffer.View; - - UInt128 localDivisor = currentDivisor128; - ArrayView1D candidateView = candidateBuffer.View.SubView(0, effectiveCount); ArrayView1D nextStartView = nextStartBuffer.View.SubView(0, 1); - GpuUInt128 gpuStartValue = new GpuUInt128(localDivisor); + GpuUInt128 gpuStartValue = new GpuUInt128(currentDivisor128); GpuUInt128 gpuStride = new GpuUInt128(twoP128); GpuUInt128 nextDivisorGpu = GenerateCandidatesOnGpu( effectiveCount, gpuStartValue, gpuStride, candidateView, - gapView, nextStartView, candidateGenerationKernel); if (effectiveCount > 0) @@ -350,22 +299,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv currentDivisor128 = (UInt128)nextDivisorGpu; } - ComputeRemaindersOnGpu(effectiveCount, remainder10, 10, gapView, deltaView, remainder10View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(effectiveCount, remainder8, 8, gapView, deltaView, remainder8View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(effectiveCount, remainder5, 5, gapView, deltaView, remainder5View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(effectiveCount, remainder3, 3, gapView, deltaView, remainder3View, remainderDeltaKernel, remainderScanKernel); - - BuildCandidateMaskOnGpu(effectiveCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, candidateMaskKernel); - - byte lastRemainder10 = effectiveCount > 0 ? FetchLastRemainder(remainder10View, effectiveCount) : (byte)0; - byte lastRemainder8 = effectiveCount > 0 ? FetchLastRemainder(remainder8View, effectiveCount) : (byte)0; - byte lastRemainder5 = effectiveCount > 0 ? FetchLastRemainder(remainder5View, effectiveCount) : (byte)0; - byte lastRemainder3 = effectiveCount > 0 ? FetchLastRemainder(remainder3View, effectiveCount) : (byte)0; - - remainder10 = AddModGpu(lastRemainder10, step10, 10); - remainder8 = AddModGpu(lastRemainder8, step8, 8); - remainder5 = AddModGpu(lastRemainder5, step5, 5); - remainder3 = AddModGpu(lastRemainder3, step3, 3); + ArrayView1D maskView = hitsBuffer.View.SubView(0, effectiveCount); + EvaluateCandidateMaskOnGpu(effectiveCount, candidateView, maskView, lastIsSevenFlag, candidateMaskKernel); int filteredCount = CompactCandidatesOnGpu( effectiveCount, @@ -379,8 +314,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { Span filteredSlice = filteredStorage[..filteredCount]; ArrayView1D filteredView = candidateView.SubView(0, filteredCount); - // The compaction kernel rewrites the accepted divisors into the front of the candidate buffer, - // so copying the compacted prefix is enough to build the CPU span for cycle checks. filteredView.CopyToCPU(ref MemoryMarshal.GetReference(filteredSlice), filteredCount); } @@ -442,76 +375,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return composite; } - private ref struct BufferSliceCache - { - private readonly MemoryBuffer1D _divisorDeltaBuffer; - private readonly MemoryBuffer1D _remainderDeltaBuffer; - private readonly MemoryBuffer1D _remainderBuffer; - private readonly MemoryBuffer1D _hitsBuffer; - private readonly int _remainderStride; - private int _length; - private ArrayView1D _divisorDelta; - private ArrayView1D _remainderDelta; - private ArrayView1D _remainder10; - private ArrayView1D _remainder8; - private ArrayView1D _remainder5; - private ArrayView1D _remainder3; - private ArrayView1D _hits; - internal BufferSliceCache( - MemoryBuffer1D divisorDeltaBuffer, - MemoryBuffer1D remainderDeltaBuffer, - MemoryBuffer1D remainderBuffer, - MemoryBuffer1D hitsBuffer, - int remainderStride) - { - _divisorDeltaBuffer = divisorDeltaBuffer; - _remainderDeltaBuffer = remainderDeltaBuffer; - _remainderBuffer = remainderBuffer; - _hitsBuffer = hitsBuffer; - _remainderStride = remainderStride; - _length = -1; - _divisorDelta = default; - _remainderDelta = default; - _remainder10 = default; - _remainder8 = default; - _remainder5 = default; - _remainder3 = default; - _hits = default; - } - internal void Ensure(int length) - { - if (_length == length) - { - return; - } - _length = length; - ArrayView1D remainderBaseView = _remainderBuffer.View; - _divisorDelta = _divisorDeltaBuffer.View.SubView(0, length); - _remainderDelta = _remainderDeltaBuffer.View.SubView(0, length); - _remainder10 = remainderBaseView.SubView(0, length); - _remainder8 = remainderBaseView.SubView(_remainderStride, length); - _remainder5 = remainderBaseView.SubView(_remainderStride * 2, length); - _remainder3 = remainderBaseView.SubView(_remainderStride * 3, length); - _hits = _hitsBuffer.View.SubView(0, length); - } - - internal ArrayView1D DivisorDelta => _divisorDelta; - - internal ArrayView1D RemainderDelta => _remainderDelta; - - internal ArrayView1D Remainder10 => _remainder10; - - internal ArrayView1D Remainder8 => _remainder8; - - internal ArrayView1D Remainder5 => _remainder5; - - internal ArrayView1D Remainder3 => _remainder3; - - internal ArrayView1D Hits => _hits; - } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static GpuUInt128 GenerateCandidatesOnGpu( @@ -519,9 +385,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv GpuUInt128 startValue, GpuUInt128 stride, ArrayView1D candidateView, - ArrayView1D gapView, ArrayView1D nextStartView, - Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> kernel) + Action, ArrayView, GpuUInt128, GpuUInt128> kernel) { // EvenPerfectBitScanner always launches this helper with a positive count, // so the guard stays commented out to keep the hot path branch-free. @@ -530,7 +395,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // return GpuUInt128.Zero; // } - kernel(count, candidateView, gapView, nextStartView, startValue, stride); + kernel(count, candidateView, nextStartView, startValue, stride); GpuUInt128 nextStart = default; nextStartView.SubView(0, 1).CopyToCPU(ref nextStart, 1); @@ -538,15 +403,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void ComputeRemaindersOnGpu( + private static void EvaluateCandidateMaskOnGpu( int count, - byte baseRemainder, - byte modulus, - ArrayView1D gapView, - ArrayView1D deltaView, - ArrayView1D remainderView, - Action, ArrayView, byte> deltaKernel, - Action, ArrayView, byte, byte> scanKernel) + ArrayView1D candidateView, + ArrayView1D maskView, + byte lastIsSevenFlag, + Action, ArrayView, byte> maskKernel) { // EvenPerfectBitScanner always supplies a positive count here, // so the guard stays commented out to keep the hot path branch-free. @@ -555,35 +417,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // return; // } - deltaKernel(count, gapView, deltaView, modulus); - SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count); - Index1D gridDim = new Index1D(1); - Index1D groupDim = new Index1D(count); - KernelConfig kernelConfig = new KernelConfig(gridDim, groupDim, sharedMemoryConfig); - - scanKernel(kernelConfig, new Index1D(count), deltaView, remainderView, baseRemainder, modulus); + maskKernel(count, candidateView, maskView, lastIsSevenFlag); } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void BuildCandidateMaskOnGpu( - int count, - ArrayView1D remainder10View, - ArrayView1D remainder8View, - ArrayView1D remainder5View, - ArrayView1D remainder3View, - byte lastIsSevenFlag, - ArrayView1D maskView, - Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> maskKernel) - { - // EvenPerfectBitScanner always supplies a positive count here, - // so the guard stays commented out to keep the hot path branch-free. - // if (count <= 0) - // { - // return; - // } - maskKernel(count, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView); - } private static int CompactCandidatesOnGpu( int count, @@ -636,23 +476,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return firstHit >= count ? -1 : firstHit; } - private static byte FetchLastRemainder(ArrayView1D remainderView, int count) - { - // Callers guard against zero counts before invoking this helper to keep the hot path branch-free. - byte value = 0; - remainderView.SubView(count - 1, 1).CopyToCPU(ref value, 1); - return value; - } - private static byte AddModGpu(byte value, byte delta, byte modulus) - { - int result = value + delta; - int firstWrap = result >= modulus ? 1 : 0; - result -= firstWrap * modulus; - int secondWrap = result >= modulus ? 1 : 0; - result -= secondWrap * modulus; - return (byte)result; - } + + private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime) { @@ -691,10 +517,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; - private Action, ArrayView> _kernel = null!; - private Action, ArrayView> _hitConversionKernel = null!; + private Action, ArrayView> _stepperKernel = null!; private MemoryBuffer1D _exponentsBuffer = null!; - private MemoryBuffer1D _resultsBuffer = null!; private MemoryBuffer1D _hitBuffer = null!; private ulong[] _hostBuffer = null!; private int _capacity; @@ -713,14 +537,37 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _disposed = false; } - private void EnsureExecutionResourcesLocked() + private void EnsureExecutionResourcesLocked(int requiredCapacity) { - _kernel = _owner.GetKernelByPrimeExponent(_accelerator); - _hitConversionKernel = _owner.GetResultConversionKernel(_accelerator); - _exponentsBuffer = _accelerator.Allocate1D(_capacity); - _resultsBuffer = _accelerator.Allocate1D(_capacity); - _hitBuffer = _accelerator.Allocate1D(_capacity); - _hostBuffer = ThreadStaticPools.UlongPool.Rent(_capacity); + if (_stepperKernel == null) + { + _stepperKernel = _owner.GetStepperKernel(_accelerator); + } + + bool allocateBuffers = _exponentsBuffer == null; + if (!allocateBuffers && _capacity < requiredCapacity) + { + _exponentsBuffer?.Dispose(); + _hitBuffer?.Dispose(); + if (_hostBuffer != null) + { + ThreadStaticPools.UlongPool.Return(_hostBuffer, clearArray: false); + _hostBuffer = null!; + } + + _exponentsBuffer = null!; + _hitBuffer = null!; + allocateBuffers = true; + } + + if (allocateBuffers) + { + int desiredCapacity = requiredCapacity > _capacity ? requiredCapacity : _capacity; + _capacity = desiredCapacity; + _exponentsBuffer = _accelerator.Allocate1D(_capacity); + _hitBuffer = _accelerator.Allocate1D(_capacity); + _hostBuffer = ThreadStaticPools.UlongPool.Rent(_capacity); + } } public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes, Span hits) @@ -731,7 +578,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } int length = primes.Length; - // EvenPerfectBitScanner always supplies at least one exponent per divisor check, so the previous zero-length guard stays commented out. + // EvenPerfectBitScanner always supplies at least one exponent per divisor check, so the guard stays commented out. // if (length == 0) // { // return; @@ -745,42 +592,51 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // return; // } + MontgomeryDivisorData cachedData = divisorData; + if (cachedData.Modulus != divisor) + { + cachedData = MontgomeryDivisorData.FromModulus(divisor); + } + + ulong cycle = divisorCycle; + // EvenPerfectBitScanner resolves divisor cycles before scheduling GPU work, so the zero-cycle fallback stays disabled. + // if (cycle == 0UL) + // { + // throw new InvalidOperationException("GPU divisor cycle length must be non-zero."); + // } + + ulong firstPrime = primes[0]; + // Precompute the initial cycle remainder on the host to keep the GPU kernel branch-free. + // EvenPerfectBitScanner resolves divisor cycles before scheduling GPU work, so cycle never equals zero here. + ulong firstRemainder = firstPrime % cycle; + Monitor.Enter(_lease.ExecutionLock); - if (_kernel == null || _hitConversionKernel == null) - { - EnsureExecutionResourcesLocked(); - } + EnsureExecutionResourcesLocked(length); - var kernel = _kernel!; - var hitConversionKernel = _hitConversionKernel!; - - ArrayView1D exponentsView = _exponentsBuffer.View; - ArrayView1D resultsView = _resultsBuffer.View; + ArrayView1D exponentView = _exponentsBuffer.View; ArrayView1D hitView = _hitBuffer.View; - Span hostSpan = _hostBuffer.AsSpan(0, _capacity); - int offset = 0; - while (offset < length) - { - int batchSize = Math.Min(_capacity, length - offset); - ReadOnlySpan batch = primes.Slice(offset, batchSize); - batch.CopyTo(hostSpan); + Span hostSpan = _hostBuffer.AsSpan(0, length); + primes.CopyTo(hostSpan); - ArrayView1D exponentSlice = exponentsView.SubView(0, batchSize); - exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), batchSize); + ArrayView1D exponentSlice = exponentView.SubView(0, length); + exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length); - ArrayView1D resultSlice = resultsView.SubView(0, batchSize); - kernel(batchSize, divisorData, exponentSlice, resultSlice); + ArrayView1D hitSliceView = hitView.SubView(0, length); + var stepperKernel = _stepperKernel!; + // EnsureExecutionResourcesLocked binds the GPU stepper kernel during the first session launch, so the null check + // remains commented out to avoid a redundant branch on the production path. + // if (stepperKernel is null) + // { + // Monitor.Exit(_lease.ExecutionLock); + // throw new InvalidOperationException("GPU stepper kernel was not initialized."); + // } - ArrayView1D hitSliceView = hitView.SubView(0, batchSize); - hitConversionKernel(batchSize, resultSlice, hitSliceView); + stepperKernel(1, cachedData, cycle, firstRemainder, exponentSlice, hitSliceView); - Span hitSlice = hits.Slice(offset, batchSize); - hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), batchSize); - - offset += batchSize; - } + Span hitSlice = hits.Slice(0, length); + hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), length); Monitor.Exit(_lease.ExecutionLock); } @@ -859,9 +715,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv HitIndexBuffer = accelerator.Allocate1D(1); CandidateBuffer = accelerator.Allocate1D(actualCapacity); NextStartBuffer = accelerator.Allocate1D(1); - DivisorDeltaBuffer = accelerator.Allocate1D(actualCapacity); - RemainderDeltaBuffer = accelerator.Allocate1D(actualCapacity); - RemainderBuffer = accelerator.Allocate1D(actualCapacity * RemainderSlotCount); FilterCountBuffer = accelerator.Allocate1D(1); ArrayPool ulongPool = ThreadStaticPools.UlongPool; Divisors = ulongPool.Rent(actualCapacity); @@ -869,8 +722,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorData = ThreadStaticPools.MontgomeryDivisorDataPool.Rent(actualCapacity); CycleCandidates = ulongPool.Rent(cycleCapacity); CycleLengths = ulongPool.Rent(cycleCapacity); - RemainderDeltaKernel = owner.GetRemainderDeltaKernel(accelerator); - RemainderScanKernel = owner.GetRemainderScanKernel(accelerator); CandidateMaskKernel = owner.GetCandidateMaskKernel(accelerator); CandidateCompactionKernel = owner.GetCandidateCompactionKernel(accelerator); CandidateGenerationKernel = owner.GetCandidateGenerationKernel(accelerator); @@ -891,11 +742,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly MemoryBuffer1D NextStartBuffer; - internal readonly MemoryBuffer1D DivisorDeltaBuffer; - - internal readonly MemoryBuffer1D RemainderDeltaBuffer; - - internal readonly MemoryBuffer1D RemainderBuffer; internal readonly MemoryBuffer1D FilterCountBuffer; internal readonly ulong[] Divisors; @@ -908,12 +754,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly ulong[] CycleLengths; - internal readonly Action, ArrayView, byte> RemainderDeltaKernel; + internal readonly Action, ArrayView, byte> CandidateMaskKernel; - internal readonly Action, ArrayView, byte, byte> RemainderScanKernel; - internal readonly Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> CandidateMaskKernel; internal readonly Action, ArrayView, ArrayView, ArrayView> CandidateCompactionKernel; - internal readonly Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> CandidateGenerationKernel; + + internal readonly Action, ArrayView, GpuUInt128, GpuUInt128> CandidateGenerationKernel; + internal readonly Action, ArrayView> HitDetectionKernel; internal readonly int Capacity; @@ -928,9 +774,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv HitIndexBuffer.Dispose(); CandidateBuffer.Dispose(); NextStartBuffer.Dispose(); - DivisorDeltaBuffer.Dispose(); - RemainderDeltaBuffer.Dispose(); - RemainderBuffer.Dispose(); FilterCountBuffer.Dispose(); ArrayPool ulongPool = ThreadStaticPools.UlongPool; ArrayPool divisorPool = ThreadStaticPools.MontgomeryDivisorDataPool; @@ -942,6 +785,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } + private sealed class AcceleratorReferenceComparer : IEqualityComparer { public static AcceleratorReferenceComparer Instance { get; } = new();