From 681e9a4a5e4c62213beca379a35f1b852f141b54 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 12 Oct 2025 22:46:58 +0200 Subject: [PATCH 1/8] Batch CPU divisor scanning across reusable buffers --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 392 ++++++++++++++---- 1 file changed, 308 insertions(+), 84 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index d269982..d64ef23 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,6 +1,6 @@ using System; +using System.Buffers; using System.Collections.Concurrent; -using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; @@ -200,12 +200,25 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - Dictionary? factorCache = null; + int batchCapacity = Math.Max(1, _batchSize); DivisorCycleCache cycleCache = DivisorCycleCache.Shared; + int preferredCycleBatch = cycleCache.PreferredBatchSize; + if (preferredCycleBatch > 0 && preferredCycleBatch < batchCapacity) + { + batchCapacity = preferredCycleBatch; + } - ulong cachedModulus = 0UL; - MontgomeryDivisorData cachedDivisorData = default; - bool hasCachedDivisorData = false; + ulong[] candidateBuffer = ArrayPool.Shared.Rent(batchCapacity); + ulong[] cycleBuffer = ArrayPool.Shared.Rent(batchCapacity); + MontgomeryDivisorData[] divisorDataBuffer = ArrayPool.Shared.Rent(batchCapacity); + byte[] remainderBuffer = ArrayPool.Shared.Rent(batchCapacity * 6); + + Span remainder10Span = remainderBuffer.AsSpan(0, batchCapacity); + Span remainder8Span = remainderBuffer.AsSpan(batchCapacity, batchCapacity); + Span remainder5Span = remainderBuffer.AsSpan(batchCapacity * 2, batchCapacity); + Span remainder3Span = remainderBuffer.AsSpan(batchCapacity * 3, batchCapacity); + Span remainder7Span = remainderBuffer.AsSpan(batchCapacity * 4, batchCapacity); + Span remainder11Span = remainderBuffer.AsSpan(batchCapacity * 5, batchCapacity); byte step10 = (byte)(step % 10UL); byte step8 = (byte)(step % 8UL); @@ -221,83 +234,151 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv byte remainder7 = (byte)(divisor % 7UL); byte remainder11 = (byte)(divisor % 11UL); - // Keep the divisibility filters aligned with the divisor-cycle generator so the - // CPU path never requests cycles that were skipped during cache creation. bool lastIsSeven = (prime & 3UL) == 3UL; - while (divisor <= limit) + try { - if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) + while (divisor <= limit) { - processedAll = false; - return false; - } + Span candidateSpan = candidateBuffer.AsSpan(0, batchCapacity); + int chunkCount = 0; + UInt128 localDivisor = divisor; - iterationCounter++; - ulong candidate = (ulong)divisor; - processedCount++; - lastProcessed = candidate; + byte localRemainder10 = remainder10; + byte localRemainder8 = remainder8; + byte localRemainder5 = remainder5; + byte localRemainder3 = remainder3; + byte localRemainder7 = remainder7; + byte localRemainder11 = remainder11; - bool admissible = lastIsSeven - ? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9) - : (remainder10 == 1 || remainder10 == 3 || remainder10 == 9); - - if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0) - { - if (!hasCachedDivisorData || cachedModulus != candidate) + while (chunkCount < batchCapacity && localDivisor <= limit) { - cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate); - cachedModulus = candidate; - hasCachedDivisorData = true; - } - - MontgomeryDivisorData divisorData = cachedDivisorData; - ulong divisorCycle; - - if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) - { - divisorCycle = cycleCache.GetCycleLength(candidate); - } - else - { - factorCache ??= new Dictionary(8); - // EvenPerfectBitScanner reaches this path only for divisors above the precomputed cycle snapshot. - // Those candidates follow q = 2 * prime * k + 1 with prime >= 5, so the downstream order calculator - // can rely on an odd divisor and a prime exponent when deriving the cycle length. - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) { - divisorCycle = cycleCache.GetCycleLength(candidate); + processedAll = false; + return false; + } + + ulong candidate = (ulong)localDivisor; + candidateSpan[chunkCount] = candidate; + remainder10Span[chunkCount] = localRemainder10; + remainder8Span[chunkCount] = localRemainder8; + remainder5Span[chunkCount] = localRemainder5; + remainder3Span[chunkCount] = localRemainder3; + remainder7Span[chunkCount] = localRemainder7; + remainder11Span[chunkCount] = localRemainder11; + + processedCount++; + lastProcessed = candidate; + iterationCounter++; + chunkCount++; + + localDivisor += step; + localRemainder10 = AddMod(localRemainder10, step10, (byte)10); + localRemainder8 = AddMod(localRemainder8, step8, (byte)8); + localRemainder5 = AddMod(localRemainder5, step5, (byte)5); + localRemainder3 = AddMod(localRemainder3, step3, (byte)3); + localRemainder7 = AddMod(localRemainder7, step7, (byte)7); + localRemainder11 = AddMod(localRemainder11, step11, (byte)11); + } + + if (chunkCount == 0) + { + break; + } + + divisor = localDivisor; + remainder10 = localRemainder10; + remainder8 = localRemainder8; + remainder5 = localRemainder5; + remainder3 = localRemainder3; + remainder7 = localRemainder7; + remainder11 = localRemainder11; + + // Reuse the candidate buffer in-place so the filtered results stay contiguous without additional arrays. + Span chunkCandidates = candidateSpan[..chunkCount]; + int filteredCount = 0; + + for (int i = 0; i < chunkCount; i++) + { + bool admissible = lastIsSeven + ? (remainder10Span[i] == 3 || remainder10Span[i] == 7 || remainder10Span[i] == 9) + : (remainder10Span[i] == 1 || remainder10Span[i] == 3 || remainder10Span[i] == 9); + + if (!admissible) + { + continue; + } + + if ((remainder8Span[i] != 1 && remainder8Span[i] != 7) + || remainder3Span[i] == 0 + || remainder5Span[i] == 0 + || remainder7Span[i] == 0 + || remainder11Span[i] == 0) + { + continue; + } + + // Store the accepted candidate at the next free slot so the front of the span becomes the filtered set. + chunkCandidates[filteredCount++] = chunkCandidates[i]; + } + + if (filteredCount == 0) + { + continue; + } + + Span filteredSpan = chunkCandidates[..filteredCount]; + Span cycleSpan = cycleBuffer.AsSpan(0, filteredCount); + cycleCache.GetCycleLengths(filteredSpan, cycleSpan); + + int admissibleCount = 0; + Span divisorDataSpan = divisorDataBuffer.AsSpan(0, filteredCount); + + for (int i = 0; i < filteredCount; i++) + { + if (cycleSpan[i] != prime) + { + continue; + } + + ulong candidate = filteredSpan[i]; + // Reuse the filtered span's leading segment to accumulate the admissible divisors for Montgomery checks. + filteredSpan[admissibleCount] = candidate; + divisorDataSpan[admissibleCount] = MontgomeryDivisorData.FromModulus(candidate); + admissibleCount++; + } + + if (admissibleCount == 0) + { + continue; + } + + Span admissibleSpan = filteredSpan[..admissibleCount]; + Span admissibleDataSpan = divisorDataSpan[..admissibleCount]; + + for (int i = 0; i < admissibleCount; i++) + { + if (CheckDivisor(prime, prime, in admissibleDataSpan[i]) != 0) + { + lastProcessed = admissibleSpan[i]; + processedAll = true; + return true; } } - - if (divisorCycle == prime) - { - // A cycle equal to the tested exponent (which is prime in this path) - // guarantees that the candidate divides the corresponding Mersenne - // number because the order of 2 modulo the divisor is exactly p. - processedAll = true; - return true; - } - - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } } - - divisor += step; - remainder10 = AddMod(remainder10, step10, (byte)10); - remainder8 = AddMod(remainder8, step8, (byte)8); - remainder5 = AddMod(remainder5, step5, (byte)5); - remainder3 = AddMod(remainder3, step3, (byte)3); - remainder7 = AddMod(remainder7, step7, (byte)7); - remainder11 = AddMod(remainder11, step11, (byte)11); + } + finally + { + ArrayPool.Shared.Return(candidateBuffer, clearArray: false); + ArrayPool.Shared.Return(cycleBuffer, clearArray: false); + ArrayPool.Shared.Return(divisorDataBuffer, clearArray: false); + ArrayPool.Shared.Return(remainderBuffer, clearArray: false); } processedAll = divisor > limit; return false; } - private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) { ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData); @@ -372,6 +453,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool _hasCachedDivisorData; private ulong _cachedDivisor; private MontgomeryDivisorData _cachedDivisorData; + private ulong[]? _primeDeltas; + private ulong[]? _cycleRemainders; + private ulong[]? _montgomeryDeltas; + private int _bufferCapacity; internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) { @@ -412,6 +497,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv cachedData = _cachedDivisorData; } + ulong modulus = cachedData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits.Clear(); + return; + } + if (divisorCycle == 0UL) { divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); @@ -422,35 +514,167 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } - // Keep these remainder steppers in place so future updates continue reusing the previously computed residues. - // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors. - var exponentStepper = new ExponentRemainderStepper(cachedData); - if (!exponentStepper.IsValidModulus) - { - hits.Clear(); - return; - } + EnsureCapacity(length); - var cycleStepper = new CycleRemainderStepper(divisorCycle); + Span deltaSpan = _primeDeltas!.AsSpan(0, length); + ComputePrimeDeltas(primes, deltaSpan); - ulong remainder = cycleStepper.Initialize(primes[0]); - hits[0] = remainder == 0UL - ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0) - : (byte)0; + Span remainderSpan = _cycleRemainders!.AsSpan(0, length); + ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); + + Span montgomerySpan = _montgomeryDeltas!.AsSpan(0, length); + ComputeMontgomeryDeltas(cachedData, primes, deltaSpan, montgomerySpan); + + ulong nPrime = cachedData.NPrime; + ulong montgomeryOne = cachedData.MontgomeryOne; + ulong currentMontgomery = montgomerySpan[0]; + + hits[0] = remainderSpan[0] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0; for (int i = 1; i < length; i++) { - remainder = cycleStepper.ComputeNext(primes[i]); - if (remainder != 0UL) + currentMontgomery = currentMontgomery.MontgomeryMultiply(montgomerySpan[i], modulus, nPrime); + hits[i] = remainderSpan[i] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0; + } + } + + private void EnsureCapacity(int requiredLength) + { + if (requiredLength <= 0) + { + return; + } + + if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null) + { + return; + } + + int newCapacity = _bufferCapacity; + if (newCapacity == 0) + { + newCapacity = 1; + } + + while (newCapacity < requiredLength) + { + newCapacity <<= 1; + } + + ReturnBuffers(); + + _primeDeltas = ArrayPool.Shared.Rent(newCapacity); + _cycleRemainders = ArrayPool.Shared.Rent(newCapacity); + _montgomeryDeltas = ArrayPool.Shared.Rent(newCapacity); + _bufferCapacity = newCapacity; + } + + private void ReturnBuffers() + { + if (_primeDeltas is not null) + { + ArrayPool.Shared.Return(_primeDeltas, clearArray: false); + _primeDeltas = null; + } + + if (_cycleRemainders is not null) + { + ArrayPool.Shared.Return(_cycleRemainders, clearArray: false); + _cycleRemainders = null; + } + + if (_montgomeryDeltas is not null) + { + ArrayPool.Shared.Return(_montgomeryDeltas, clearArray: false); + _montgomeryDeltas = null; + } + + _bufferCapacity = 0; + } + + private static void ComputePrimeDeltas(ReadOnlySpan primes, Span deltas) + { + if (primes.Length == 0) + { + return; + } + + ulong previous = primes[0]; + deltas[0] = previous; + + for (int i = 1; i < primes.Length; i++) + { + ulong current = primes[i]; + if (current <= previous) { - hits[i] = 0; - continue; + throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing."); } - hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; + deltas[i] = current - previous; + previous = current; } } + private static void ComputeCycleRemainders(ulong cycleLength, ReadOnlySpan primes, ReadOnlySpan deltas, Span remainders) + { + if (cycleLength == 0UL || primes.Length == 0) + { + remainders.Clear(); + return; + } + + ulong remainder = primes[0] % cycleLength; + remainders[0] = remainder; + + for (int i = 1; i < primes.Length; i++) + { + remainder = AddMod(remainder, deltas[i], cycleLength); + remainders[i] = remainder; + } + } + + private static void ComputeMontgomeryDeltas(in MontgomeryDivisorData divisorData, ReadOnlySpan primes, ReadOnlySpan deltas, Span montgomeryValues) + { + if (primes.Length == 0) + { + montgomeryValues.Clear(); + return; + } + + montgomeryValues[0] = primes[0].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true); + + for (int i = 1; i < primes.Length; i++) + { + montgomeryValues[i] = deltas[i].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true); + } + } + + private static ulong AddMod(ulong value, ulong delta, ulong modulus) + { + if (modulus == 0UL) + { + return 0UL; + } + + if (ulong.MaxValue - value < delta) + { + UInt128 extended = (UInt128)value + delta; + return (ulong)(extended % modulus); + } + + ulong sum = value + delta; + if (sum >= modulus) + { + sum -= modulus; + if (sum >= modulus) + { + sum %= modulus; + } + } + + return sum; + } + public void Dispose() { if (_disposed) From 0902b18721a64a9060501f4b6c134b26b792ef00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 01:35:21 +0200 Subject: [PATCH 2/8] Add GPU delta routing option for by-divisor mode --- EvenPerfectBitScanner/Program.cs | 33 ++++++- PerfectNumbers.Core/ByDivisorDeltasDevice.cs | 7 ++ ...MersenneNumberDivisorByDivisorCpuTester.cs | 87 ++++++++++++++++++- 3 files changed, 123 insertions(+), 4 deletions(-) create mode 100644 PerfectNumbers.Core/ByDivisorDeltasDevice.cs diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 7cefd12..e813d6b 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -75,6 +75,7 @@ internal static class Program bool useResidue = false; // M_p test via residue divisors (replaces LL/incremental) bool useDivisor = false; // M_p divisibility by specific divisor bool useByDivisor = false; // Iterative divisor scan across primes UInt128 divisor = UInt128.Zero; + ByDivisorDeltasDevice byDivisorDeltasDevice = ByDivisorDeltasDevice.Cpu; // Device routing bool useGpuCycles = true; bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device @@ -182,6 +183,23 @@ internal static class Program { mersenneOnGpu = !arg.AsSpan(arg.IndexOf('=') + 1).Equals("cpu", StringComparison.OrdinalIgnoreCase); } + else if (arg.StartsWith("--bydivisor-deltas-device=", StringComparison.OrdinalIgnoreCase)) + { + ReadOnlySpan deviceValue = arg.AsSpan(arg.IndexOf('=') + 1); + if (deviceValue.Equals("gpu", StringComparison.OrdinalIgnoreCase)) + { + byDivisorDeltasDevice = ByDivisorDeltasDevice.Gpu; + } + else if (deviceValue.Equals("cpu", StringComparison.OrdinalIgnoreCase)) + { + byDivisorDeltasDevice = ByDivisorDeltasDevice.Cpu; + } + else + { + Console.WriteLine("Invalid --bydivisor-deltas-device value. Expected cpu or gpu."); + return; + } + } else if (arg.StartsWith("--prime-test-limit=", StringComparison.OrdinalIgnoreCase)) { ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); @@ -531,6 +549,10 @@ internal static class Program ? new MersenneNumberDivisorByDivisorGpuTester() : new MersenneNumberDivisorByDivisorCpuTester(); _byDivisorTester.BatchSize = scanBatchSize; + if (_byDivisorTester is MersenneNumberDivisorByDivisorCpuTester cpuByDivisorTester) + { + cpuByDivisorTester.DeltasDevice = byDivisorDeltasDevice; + } } // Load RLE blacklist (optional) @@ -557,6 +579,9 @@ internal static class Program Console.WriteLine("Initialization..."); // Compose a results file name that encodes configuration (before opening file) + string byDivisorDeltasLabel = mersenneOnGpu + ? "gpu" + : (byDivisorDeltasDevice == ByDivisorDeltasDevice.Gpu ? "cpu-gpu" : "cpu"); var builtName = BuildResultsFileName( useBitTransform, threadCount, @@ -577,7 +602,8 @@ internal static class Program NttGpuMath.ReductionMode, mersenneOnGpu ? "gpu" : "cpu", (GpuContextPool.ForceCpu ? "cpu" : "gpu"), - orderOnGpu ? "gpu" : "cpu"); + orderOnGpu ? "gpu" : "cpu", + byDivisorDeltasLabel); if (!string.IsNullOrEmpty(_resultsPrefix)) { @@ -984,6 +1010,7 @@ internal static class Program Console.WriteLine(" --divisor-cycles-continue continue divisor cycles generation"); Console.WriteLine(" --divisor-cycles-limit= cycle search iterations when --mersenne=divisor"); Console.WriteLine(" --prime-test-limit= Time limit for --mersenne=bydivisor prime checks (e.g. 5s, 500ms)"); + Console.WriteLine(" --bydivisor-deltas-device=cpu|gpu Device for --mersenne=bydivisor residue batches (default: cpu)"); Console.WriteLine(" --use-order test primality via q order"); Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU"); // mod-automaton removed @@ -993,7 +1020,7 @@ internal static class Program Console.WriteLine(" --help, -help, --?, -?, /? show this help message"); } - private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool useByDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice) + private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool useByDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice, string byDivisorDeltasDevice) { string inc = bitInc ? "bit" : "add"; string mers = useDivisorFlag @@ -1001,7 +1028,7 @@ internal static class Program : (useLucasFlag ? "lucas" : (useByDivisorFlag - ? "bydivisor" + ? $"bydivisor-{byDivisorDeltasDevice}" : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental"))); string ntt = nttBackend == NttBackend.Staged ? "staged" : "reference"; string red = reduction switch { ModReductionMode.Mont64 => "mont64", ModReductionMode.Barrett128 => "barrett128", ModReductionMode.GpuUInt128 => "uint128", _ => "auto" }; diff --git a/PerfectNumbers.Core/ByDivisorDeltasDevice.cs b/PerfectNumbers.Core/ByDivisorDeltasDevice.cs new file mode 100644 index 0000000..2011b25 --- /dev/null +++ b/PerfectNumbers.Core/ByDivisorDeltasDevice.cs @@ -0,0 +1,7 @@ +namespace PerfectNumbers.Core; + +public enum ByDivisorDeltasDevice +{ + Cpu, + Gpu +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index d64ef23..4ee09a7 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -3,6 +3,7 @@ using System.Buffers; using System.Collections.Concurrent; using System.Numerics; using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Cpu; @@ -15,20 +16,51 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool _isConfigured; private int _batchSize = 1_024; + private ulong _configuredMaxPrime; + private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu; + private MersenneNumberDivisorByDivisorGpuTester? _gpuDeltasTester; + private bool _gpuTesterConfigured; + public int BatchSize { get => _batchSize; - set => _batchSize = Math.Max(1, value); + set + { + int sanitized = Math.Max(1, value); + _batchSize = sanitized; + + lock (_sync) + { + if (_gpuDeltasTester is not null) + { + _gpuDeltasTester.GpuBatchSize = sanitized; + } + } + } + } + + public ByDivisorDeltasDevice DeltasDevice + { + get => _deltasDevice; + set + { + lock (_sync) + { + _deltasDevice = value; + } + } } public void ConfigureFromMaxPrime(ulong maxPrime) { lock (_sync) { + _configuredMaxPrime = maxPrime; _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); _lastStatusDivisor = 0UL; _isConfigured = true; + _gpuTesterConfigured = false; } } @@ -135,6 +167,31 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } + internal IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateGpuDeltasSession() + { + MersenneNumberDivisorByDivisorGpuTester tester; + + lock (_sync) + { + if (!_isConfigured) + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } + + tester = _gpuDeltasTester ??= new MersenneNumberDivisorByDivisorGpuTester(); + + if (!_gpuTesterConfigured) + { + tester.ConfigureFromMaxPrime(_configuredMaxPrime); + _gpuTesterConfigured = true; + } + + tester.GpuBatchSize = _batchSize; + } + + return tester.CreateDivisorSession(); + } + public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession() { lock (_sync) @@ -457,6 +514,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private ulong[]? _cycleRemainders; private ulong[]? _montgomeryDeltas; private int _bufferCapacity; + private IMersenneNumberDivisorByDivisorTester.IDivisorScanSession? _gpuSession; internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) { @@ -514,6 +572,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } + if (_owner._deltasDevice == ByDivisorDeltasDevice.Gpu) + { + EnsureGpuSession(); + _gpuSession!.CheckDivisor(divisor, cachedData, divisorCycle, primes, hits); + return; + } + EnsureCapacity(length); Span deltaSpan = _primeDeltas!.AsSpan(0, length); @@ -538,6 +603,25 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } + private void EnsureGpuSession() + { + if (_gpuSession is null) + { + _gpuSession = _owner.CreateGpuDeltasSession(); + } + } + + private void ReleaseGpuSession() + { + if (_gpuSession is null) + { + return; + } + + _gpuSession.Dispose(); + _gpuSession = null; + } + private void EnsureCapacity(int requiredLength) { if (requiredLength <= 0) @@ -683,6 +767,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } _disposed = true; + ReleaseGpuSession(); _owner.ReturnSession(this); } } From 41e6703926cbeaad4e0a1d51519fc39cb5e0906a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 01:35:26 +0200 Subject: [PATCH 3/8] Add dedicated GPU residue evaluator for by-divisor deltas --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 79 ++++-------- ...ersenneNumberDivisorResidueGpuEvaluator.cs | 114 ++++++++++++++++++ 2 files changed, 140 insertions(+), 53 deletions(-) create mode 100644 PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 4ee09a7..70cc274 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -18,8 +18,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private ulong _configuredMaxPrime; private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu; - private MersenneNumberDivisorByDivisorGpuTester? _gpuDeltasTester; - private bool _gpuTesterConfigured; public int BatchSize @@ -29,14 +27,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv { int sanitized = Math.Max(1, value); _batchSize = sanitized; - - lock (_sync) - { - if (_gpuDeltasTester is not null) - { - _gpuDeltasTester.GpuBatchSize = sanitized; - } - } } } @@ -60,7 +50,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); _lastStatusDivisor = 0UL; _isConfigured = true; - _gpuTesterConfigured = false; } } @@ -167,31 +156,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } - internal IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateGpuDeltasSession() - { - MersenneNumberDivisorByDivisorGpuTester tester; - - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - tester = _gpuDeltasTester ??= new MersenneNumberDivisorByDivisorGpuTester(); - - if (!_gpuTesterConfigured) - { - tester.ConfigureFromMaxPrime(_configuredMaxPrime); - _gpuTesterConfigured = true; - } - - tester.GpuBatchSize = _batchSize; - } - - return tester.CreateDivisorSession(); - } - public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession() { lock (_sync) @@ -514,7 +478,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private ulong[]? _cycleRemainders; private ulong[]? _montgomeryDeltas; private int _bufferCapacity; - private IMersenneNumberDivisorByDivisorTester.IDivisorScanSession? _gpuSession; + private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator; internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) { @@ -572,13 +536,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } - if (_owner._deltasDevice == ByDivisorDeltasDevice.Gpu) - { - EnsureGpuSession(); - _gpuSession!.CheckDivisor(divisor, cachedData, divisorCycle, primes, hits); - return; - } - EnsureCapacity(length); Span deltaSpan = _primeDeltas!.AsSpan(0, length); @@ -587,11 +544,27 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv Span remainderSpan = _cycleRemainders!.AsSpan(0, length); ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); + ulong montgomeryOne = cachedData.MontgomeryOne; + + if (_owner._deltasDevice == ByDivisorDeltasDevice.Gpu) + { + EnsureGpuEvaluator(); + Span residueSpan = _montgomeryDeltas!.AsSpan(0, length); + // Reuse the Montgomery delta buffer to hold the final residues returned by the GPU evaluator. + _gpuResidueEvaluator!.ComputeResidues(primes, cachedData, residueSpan); + + for (int i = 0; i < length; i++) + { + hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0; + } + + return; + } + Span montgomerySpan = _montgomeryDeltas!.AsSpan(0, length); ComputeMontgomeryDeltas(cachedData, primes, deltaSpan, montgomerySpan); ulong nPrime = cachedData.NPrime; - ulong montgomeryOne = cachedData.MontgomeryOne; ulong currentMontgomery = montgomerySpan[0]; hits[0] = remainderSpan[0] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0; @@ -603,23 +576,23 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } - private void EnsureGpuSession() + private void EnsureGpuEvaluator() { - if (_gpuSession is null) + if (_gpuResidueEvaluator is null) { - _gpuSession = _owner.CreateGpuDeltasSession(); + _gpuResidueEvaluator = new MersenneNumberDivisorResidueGpuEvaluator(_owner._batchSize); } } - private void ReleaseGpuSession() + private void ReleaseGpuEvaluator() { - if (_gpuSession is null) + if (_gpuResidueEvaluator is null) { return; } - _gpuSession.Dispose(); - _gpuSession = null; + _gpuResidueEvaluator.Dispose(); + _gpuResidueEvaluator = null; } private void EnsureCapacity(int requiredLength) @@ -767,7 +740,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } _disposed = true; - ReleaseGpuSession(); + ReleaseGpuEvaluator(); _owner.ReturnSession(this); } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs new file mode 100644 index 0000000..8fdd163 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs @@ -0,0 +1,114 @@ +using System; +using System.Buffers; +using System.Runtime.InteropServices; +using PerfectNumbers.Core; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable +{ + private readonly GpuContextPool.GpuContextLease _lease; + private readonly Accelerator _accelerator; + private readonly Action, ArrayView> _kernel; + private MemoryBuffer1D _exponentBuffer; + private MemoryBuffer1D _resultBuffer; + private ulong[] _hostBuffer; + private int _capacity; + private bool _disposed; + + internal MersenneNumberDivisorResidueGpuEvaluator(int batchSize) + { + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(ComputeResiduesKernel); + _capacity = Math.Max(1, batchSize); + _exponentBuffer = _accelerator.Allocate1D(_capacity); + _resultBuffer = _accelerator.Allocate1D(_capacity); + _hostBuffer = ArrayPool.Shared.Rent(_capacity); + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _exponentBuffer.Dispose(); + _resultBuffer.Dispose(); + ArrayPool.Shared.Return(_hostBuffer, clearArray: false); + _lease.Dispose(); + } + + internal void ComputeResidues(ReadOnlySpan exponents, in MontgomeryDivisorData divisorData, Span destination) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(MersenneNumberDivisorResidueGpuEvaluator)); + } + + int length = exponents.Length; + if (length == 0) + { + return; + } + + EnsureCapacity(length); + + Span hostSpan = _hostBuffer.AsSpan(0, length); + exponents.CopyTo(hostSpan); + + ArrayView1D exponentView = _exponentBuffer.View.SubView(0, length); + exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length); + + ArrayView1D resultView = _resultBuffer.View.SubView(0, length); + _kernel(length, divisorData, exponentView, resultView); + resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostSpan), length); + + hostSpan.CopyTo(destination); + } + + private void EnsureCapacity(int required) + { + if (required <= _capacity) + { + return; + } + + int newCapacity = _capacity; + if (newCapacity == 0) + { + newCapacity = 1; + } + + while (newCapacity < required) + { + newCapacity <<= 1; + } + + _exponentBuffer.Dispose(); + _resultBuffer.Dispose(); + ArrayPool.Shared.Return(_hostBuffer, clearArray: false); + + _capacity = newCapacity; + _exponentBuffer = _accelerator.Allocate1D(_capacity); + _resultBuffer = _accelerator.Allocate1D(_capacity); + _hostBuffer = ArrayPool.Shared.Rent(_capacity); + } + + private static void ComputeResiduesKernel(Index1D index, MontgomeryDivisorData divisorData, ArrayView exponents, ArrayView results) + { + int idx = index; + int length = (int)exponents.Length; + if (idx >= length) + { + return; + } + + ulong exponent = exponents[idx]; + results[idx] = exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: true); + } +} From cfedfd12c8db7ec7c1be103c774cc7c1e6fc9ce7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 01:35:32 +0200 Subject: [PATCH 4/8] Restore CPU by-divisor flow when GPU deltas disabled --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 258 +++++++++++++++--- 1 file changed, 215 insertions(+), 43 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 70cc274..81f5df7 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,6 +1,7 @@ using System; using System.Buffers; using System.Collections.Concurrent; +using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; @@ -187,6 +188,162 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv out bool processedAll, out ulong processedCount, TimeSpan? timeLimit) + { + if (_deltasDevice == ByDivisorDeltasDevice.Cpu) + { + return CheckDivisorsCpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); + } + + return CheckDivisorsGpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); + } + + private bool CheckDivisorsCpu( + ulong prime, + ulong allowedMax, + out ulong lastProcessed, + out bool processedAll, + out ulong processedCount, + TimeSpan? timeLimit) + { + lastProcessed = 0UL; + processedCount = 0UL; + processedAll = false; + + PrimeTestTimeLimit limitGuard; + if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard)) + { + return false; + } + + bool enforceLimit = limitGuard.IsActive; + ulong iterationCounter = 0UL; + + if (allowedMax < 3UL) + { + return false; + } + + UInt128 step = (UInt128)prime << 1; + if (step == UInt128.Zero) + { + processedAll = true; + return false; + } + + UInt128 limit = allowedMax; + UInt128 divisor = step + UInt128.One; + if (divisor > limit) + { + processedAll = true; + return false; + } + + Dictionary? factorCache = null; + DivisorCycleCache cycleCache = DivisorCycleCache.Shared; + + ulong cachedModulus = 0UL; + MontgomeryDivisorData cachedDivisorData = default; + bool hasCachedDivisorData = false; + + byte step10 = (byte)(step % 10UL); + byte step8 = (byte)(step % 8UL); + byte step5 = (byte)(step % 5UL); + byte step3 = (byte)(step % 3UL); + byte step7 = (byte)(step % 7UL); + byte step11 = (byte)(step % 11UL); + + byte remainder10 = (byte)(divisor % 10UL); + byte remainder8 = (byte)(divisor % 8UL); + byte remainder5 = (byte)(divisor % 5UL); + byte remainder3 = (byte)(divisor % 3UL); + byte remainder7 = (byte)(divisor % 7UL); + byte remainder11 = (byte)(divisor % 11UL); + + // Keep the divisibility filters aligned with the divisor-cycle generator so the + // CPU path never requests cycles that were skipped during cache creation. + bool lastIsSeven = (prime & 3UL) == 3UL; + + while (divisor <= limit) + { + if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) + { + processedAll = false; + return false; + } + + iterationCounter++; + ulong candidate = (ulong)divisor; + processedCount++; + lastProcessed = candidate; + + bool admissible = lastIsSeven + ? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9) + : (remainder10 == 1 || remainder10 == 3 || remainder10 == 9); + + if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0) + { + if (!hasCachedDivisorData || cachedModulus != candidate) + { + cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate); + cachedModulus = candidate; + hasCachedDivisorData = true; + } + + MontgomeryDivisorData divisorData = cachedDivisorData; + ulong divisorCycle; + + if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) + { + divisorCycle = cycleCache.GetCycleLength(candidate); + } + else + { + factorCache ??= new Dictionary(8); + // EvenPerfectBitScanner reaches this path only for divisors above the precomputed cycle snapshot. + // Those candidates follow q = 2 * prime * k + 1 with prime >= 5, so the downstream order calculator + // can rely on an odd divisor and a prime exponent when deriving the cycle length. + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + { + divisorCycle = cycleCache.GetCycleLength(candidate); + } + } + + if (divisorCycle == prime) + { + // A cycle equal to the tested exponent (which is prime in this path) + // guarantees that the candidate divides the corresponding Mersenne + // number because the order of 2 modulo the divisor is exactly p. + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + divisor += step; + remainder10 = AddMod(remainder10, step10, (byte)10); + remainder8 = AddMod(remainder8, step8, (byte)8); + remainder5 = AddMod(remainder5, step5, (byte)5); + remainder3 = AddMod(remainder3, step3, (byte)3); + remainder7 = AddMod(remainder7, step7, (byte)7); + remainder11 = AddMod(remainder11, step11, (byte)11); + } + + processedAll = divisor > limit; + return false; + } + + + private bool CheckDivisorsGpu( + ulong prime, + ulong allowedMax, + out ulong lastProcessed, + out bool processedAll, + out ulong processedCount, + TimeSpan? timeLimit) { lastProcessed = 0UL; processedCount = 0UL; @@ -536,6 +693,58 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } + if (_owner._deltasDevice == ByDivisorDeltasDevice.Cpu) + { + CheckDivisorCpu(cachedData, divisorCycle, primes, hits); + return; + } + + CheckDivisorGpu(cachedData, divisorCycle, primes, hits, length); + } + + private static void CheckDivisorCpu( + in MontgomeryDivisorData divisorData, + ulong divisorCycle, + ReadOnlySpan primes, + Span hits) + { + // Keep these remainder steppers in place so future updates continue reusing the previously computed residues. + // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors. + var exponentStepper = new ExponentRemainderStepper(divisorData); + if (!exponentStepper.IsValidModulus) + { + hits.Clear(); + return; + } + + var cycleStepper = new CycleRemainderStepper(divisorCycle); + + ulong remainder = cycleStepper.Initialize(primes[0]); + hits[0] = remainder == 0UL + ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0) + : (byte)0; + + int length = primes.Length; + for (int i = 1; i < length; i++) + { + remainder = cycleStepper.ComputeNext(primes[i]); + if (remainder != 0UL) + { + hits[i] = 0; + continue; + } + + hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; + } + } + + private void CheckDivisorGpu( + in MontgomeryDivisorData divisorData, + ulong divisorCycle, + ReadOnlySpan primes, + Span hits, + int length) + { EnsureCapacity(length); Span deltaSpan = _primeDeltas!.AsSpan(0, length); @@ -544,35 +753,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv Span remainderSpan = _cycleRemainders!.AsSpan(0, length); ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); - ulong montgomeryOne = cachedData.MontgomeryOne; + EnsureGpuEvaluator(); + Span residueSpan = _montgomeryDeltas!.AsSpan(0, length); + _gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan); - if (_owner._deltasDevice == ByDivisorDeltasDevice.Gpu) + ulong montgomeryOne = divisorData.MontgomeryOne; + for (int i = 0; i < length; i++) { - EnsureGpuEvaluator(); - Span residueSpan = _montgomeryDeltas!.AsSpan(0, length); - // Reuse the Montgomery delta buffer to hold the final residues returned by the GPU evaluator. - _gpuResidueEvaluator!.ComputeResidues(primes, cachedData, residueSpan); - - for (int i = 0; i < length; i++) - { - hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0; - } - - return; - } - - Span montgomerySpan = _montgomeryDeltas!.AsSpan(0, length); - ComputeMontgomeryDeltas(cachedData, primes, deltaSpan, montgomerySpan); - - ulong nPrime = cachedData.NPrime; - ulong currentMontgomery = montgomerySpan[0]; - - hits[0] = remainderSpan[0] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0; - - for (int i = 1; i < length; i++) - { - currentMontgomery = currentMontgomery.MontgomeryMultiply(montgomerySpan[i], modulus, nPrime); - hits[i] = remainderSpan[i] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0; + hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0; } } @@ -690,22 +878,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } - private static void ComputeMontgomeryDeltas(in MontgomeryDivisorData divisorData, ReadOnlySpan primes, ReadOnlySpan deltas, Span montgomeryValues) - { - if (primes.Length == 0) - { - montgomeryValues.Clear(); - return; - } - - montgomeryValues[0] = primes[0].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true); - - for (int i = 1; i < primes.Length; i++) - { - montgomeryValues[i] = deltas[i].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true); - } - } - private static ulong AddMod(ulong value, ulong delta, ulong modulus) { if (modulus == 0UL) From cd47b122a57d051faf1200f7e60d32a0574c180f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 01:35:36 +0200 Subject: [PATCH 5/8] Batch GPU remainder filtering for by-divisor scan --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 259 +++++++++-------- ...senneNumberDivisorCandidateGpuEvaluator.cs | 265 ++++++++++++++++++ 2 files changed, 402 insertions(+), 122 deletions(-) create mode 100644 PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 81f5df7..65f0a9c 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -363,21 +363,29 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - UInt128 step = (UInt128)prime << 1; - if (step == UInt128.Zero) + UInt128 step128 = (UInt128)prime << 1; + if (step128 == UInt128.Zero) { processedAll = true; return false; } + if (step128 > ulong.MaxValue) + { + return CheckDivisorsCpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); + } + UInt128 limit = allowedMax; - UInt128 divisor = step + UInt128.One; - if (divisor > limit) + UInt128 divisor128 = step128 + UInt128.One; + if (divisor128 > limit) { processedAll = true; return false; } + ulong step = (ulong)step128; + ulong limit64 = allowedMax; + int batchCapacity = Math.Max(1, _batchSize); DivisorCycleCache cycleCache = DivisorCycleCache.Shared; int preferredCycleBatch = cycleCache.PreferredBatchSize; @@ -386,17 +394,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv batchCapacity = preferredCycleBatch; } - ulong[] candidateBuffer = ArrayPool.Shared.Rent(batchCapacity); - ulong[] cycleBuffer = ArrayPool.Shared.Rent(batchCapacity); - MontgomeryDivisorData[] divisorDataBuffer = ArrayPool.Shared.Rent(batchCapacity); - byte[] remainderBuffer = ArrayPool.Shared.Rent(batchCapacity * 6); + var candidateEvaluator = new MersenneNumberDivisorCandidateGpuEvaluator(batchCapacity); - Span remainder10Span = remainderBuffer.AsSpan(0, batchCapacity); - Span remainder8Span = remainderBuffer.AsSpan(batchCapacity, batchCapacity); - Span remainder5Span = remainderBuffer.AsSpan(batchCapacity * 2, batchCapacity); - Span remainder3Span = remainderBuffer.AsSpan(batchCapacity * 3, batchCapacity); - Span remainder7Span = remainderBuffer.AsSpan(batchCapacity * 4, batchCapacity); - Span remainder11Span = remainderBuffer.AsSpan(batchCapacity * 5, batchCapacity); + ulong[] candidateBuffer = ArrayPool.Shared.Rent(batchCapacity); + byte[] maskBuffer = ArrayPool.Shared.Rent(batchCapacity); + + Dictionary? factorCache = null; + + ulong cachedModulus = 0UL; + MontgomeryDivisorData cachedDivisorData = default; + bool hasCachedDivisorData = false; byte step10 = (byte)(step % 10UL); byte step8 = (byte)(step % 8UL); @@ -405,31 +412,49 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv byte step7 = (byte)(step % 7UL); byte step11 = (byte)(step % 11UL); - byte remainder10 = (byte)(divisor % 10UL); - byte remainder8 = (byte)(divisor % 8UL); - byte remainder5 = (byte)(divisor % 5UL); - byte remainder3 = (byte)(divisor % 3UL); - byte remainder7 = (byte)(divisor % 7UL); - byte remainder11 = (byte)(divisor % 11UL); + byte remainder10 = (byte)((ulong)divisor128 % 10UL); + byte remainder8 = (byte)((ulong)divisor128 % 8UL); + byte remainder5 = (byte)((ulong)divisor128 % 5UL); + byte remainder3 = (byte)((ulong)divisor128 % 3UL); + byte remainder7 = (byte)((ulong)divisor128 % 7UL); + byte remainder11 = (byte)((ulong)divisor128 % 11UL); bool lastIsSeven = (prime & 3UL) == 3UL; try { - while (divisor <= limit) + while (divisor128 <= limit) { - Span candidateSpan = candidateBuffer.AsSpan(0, batchCapacity); - int chunkCount = 0; - UInt128 localDivisor = divisor; + int chunkCount = ComputeChunkCount(divisor128, limit, step128, batchCapacity); + if (chunkCount <= 0) + { + break; + } - byte localRemainder10 = remainder10; - byte localRemainder8 = remainder8; - byte localRemainder5 = remainder5; - byte localRemainder3 = remainder3; - byte localRemainder7 = remainder7; - byte localRemainder11 = remainder11; + Span candidateSpan = candidateBuffer.AsSpan(0, chunkCount); + Span maskSpan = maskBuffer.AsSpan(0, chunkCount); - while (chunkCount < batchCapacity && localDivisor <= limit) + candidateEvaluator.EvaluateCandidates( + (ulong)divisor128, + step, + limit64, + remainder10, + remainder8, + remainder5, + remainder3, + remainder7, + remainder11, + step10, + step8, + step5, + step3, + step7, + step11, + lastIsSeven, + candidateSpan, + maskSpan); + + for (int i = 0; i < chunkCount; i++) { if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) { @@ -437,126 +462,116 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - ulong candidate = (ulong)localDivisor; - candidateSpan[chunkCount] = candidate; - remainder10Span[chunkCount] = localRemainder10; - remainder8Span[chunkCount] = localRemainder8; - remainder5Span[chunkCount] = localRemainder5; - remainder3Span[chunkCount] = localRemainder3; - remainder7Span[chunkCount] = localRemainder7; - remainder11Span[chunkCount] = localRemainder11; + iterationCounter++; + ulong candidate = candidateSpan[i]; processedCount++; lastProcessed = candidate; - iterationCounter++; - chunkCount++; - localDivisor += step; - localRemainder10 = AddMod(localRemainder10, step10, (byte)10); - localRemainder8 = AddMod(localRemainder8, step8, (byte)8); - localRemainder5 = AddMod(localRemainder5, step5, (byte)5); - localRemainder3 = AddMod(localRemainder3, step3, (byte)3); - localRemainder7 = AddMod(localRemainder7, step7, (byte)7); - localRemainder11 = AddMod(localRemainder11, step11, (byte)11); - } - - if (chunkCount == 0) - { - break; - } - - divisor = localDivisor; - remainder10 = localRemainder10; - remainder8 = localRemainder8; - remainder5 = localRemainder5; - remainder3 = localRemainder3; - remainder7 = localRemainder7; - remainder11 = localRemainder11; - - // Reuse the candidate buffer in-place so the filtered results stay contiguous without additional arrays. - Span chunkCandidates = candidateSpan[..chunkCount]; - int filteredCount = 0; - - for (int i = 0; i < chunkCount; i++) - { - bool admissible = lastIsSeven - ? (remainder10Span[i] == 3 || remainder10Span[i] == 7 || remainder10Span[i] == 9) - : (remainder10Span[i] == 1 || remainder10Span[i] == 3 || remainder10Span[i] == 9); - - if (!admissible) + if (maskSpan[i] == 0) { continue; } - if ((remainder8Span[i] != 1 && remainder8Span[i] != 7) - || remainder3Span[i] == 0 - || remainder5Span[i] == 0 - || remainder7Span[i] == 0 - || remainder11Span[i] == 0) + if (!hasCachedDivisorData || cachedModulus != candidate) { - continue; + cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate); + cachedModulus = candidate; + hasCachedDivisorData = true; } - // Store the accepted candidate at the next free slot so the front of the span becomes the filtered set. - chunkCandidates[filteredCount++] = chunkCandidates[i]; - } + MontgomeryDivisorData divisorData = cachedDivisorData; + ulong divisorCycle; - if (filteredCount == 0) - { - continue; - } - - Span filteredSpan = chunkCandidates[..filteredCount]; - Span cycleSpan = cycleBuffer.AsSpan(0, filteredCount); - cycleCache.GetCycleLengths(filteredSpan, cycleSpan); - - int admissibleCount = 0; - Span divisorDataSpan = divisorDataBuffer.AsSpan(0, filteredCount); - - for (int i = 0; i < filteredCount; i++) - { - if (cycleSpan[i] != prime) + if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) { - continue; + divisorCycle = cycleCache.GetCycleLength(candidate); + } + else + { + factorCache ??= new Dictionary(8); + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + { + divisorCycle = cycleCache.GetCycleLength(candidate); + } } - ulong candidate = filteredSpan[i]; - // Reuse the filtered span's leading segment to accumulate the admissible divisors for Montgomery checks. - filteredSpan[admissibleCount] = candidate; - divisorDataSpan[admissibleCount] = MontgomeryDivisorData.FromModulus(candidate); - admissibleCount++; - } - - if (admissibleCount == 0) - { - continue; - } - - Span admissibleSpan = filteredSpan[..admissibleCount]; - Span admissibleDataSpan = divisorDataSpan[..admissibleCount]; - - for (int i = 0; i < admissibleCount; i++) - { - if (CheckDivisor(prime, prime, in admissibleDataSpan[i]) != 0) + if (divisorCycle == prime) { - lastProcessed = admissibleSpan[i]; processedAll = true; return true; } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } } + + divisor128 += step128 * (UInt128)chunkCount; + remainder10 = AdvanceRemainder(remainder10, step10, 10, chunkCount); + remainder8 = AdvanceRemainder(remainder8, step8, 8, chunkCount); + remainder5 = AdvanceRemainder(remainder5, step5, 5, chunkCount); + remainder3 = AdvanceRemainder(remainder3, step3, 3, chunkCount); + remainder7 = AdvanceRemainder(remainder7, step7, 7, chunkCount); + remainder11 = AdvanceRemainder(remainder11, step11, 11, chunkCount); } } finally { ArrayPool.Shared.Return(candidateBuffer, clearArray: false); - ArrayPool.Shared.Return(cycleBuffer, clearArray: false); - ArrayPool.Shared.Return(divisorDataBuffer, clearArray: false); - ArrayPool.Shared.Return(remainderBuffer, clearArray: false); + ArrayPool.Shared.Return(maskBuffer, clearArray: false); + candidateEvaluator.Dispose(); } - processedAll = divisor > limit; + processedAll = divisor128 > limit; return false; } + + private static int ComputeChunkCount(UInt128 currentDivisor, UInt128 limit, UInt128 step, int maxCount) + { + if (maxCount <= 0 || step == UInt128.Zero || currentDivisor > limit) + { + return 0; + } + + UInt128 distance = limit - currentDivisor; + UInt128 stepsAvailable = distance / step; + UInt128 total = stepsAvailable + UInt128.One; + if (total > (UInt128)maxCount) + { + total = (UInt128)maxCount; + } + + if (total > (UInt128)int.MaxValue) + { + return int.MaxValue; + } + + return (int)total; + } + + private static byte AdvanceRemainder(byte remainder, byte step, byte modulus, int count) + { + if (modulus == 0 || count <= 0) + { + return remainder; + } + + int modulusValue = modulus; + int remainderValue = remainder % modulusValue; + int stepValue = step % modulusValue; + + long delta = (long)stepValue * count; + int adjusted = (int)((remainderValue + delta) % modulusValue); + if (adjusted < 0) + { + adjusted += modulusValue; + } + + return (byte)adjusted; + } + private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) { ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs new file mode 100644 index 0000000..9a75acb --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs @@ -0,0 +1,265 @@ +using System; +using System.Buffers; +using System.Runtime.InteropServices; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable +{ + private readonly GpuContextPool.GpuContextLease _lease; + private readonly Accelerator _accelerator; + private readonly Action, ArrayView> _kernel; + private MemoryBuffer1D _candidateBuffer; + private MemoryBuffer1D _maskBuffer; + private ulong[] _hostCandidates; + private byte[] _hostMask; + private int _capacity; + private bool _disposed; + + internal MersenneNumberDivisorCandidateGpuEvaluator(int batchSize) + { + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(EvaluateCandidatesKernel); + _capacity = Math.Max(1, batchSize); + _candidateBuffer = _accelerator.Allocate1D(_capacity); + _maskBuffer = _accelerator.Allocate1D(_capacity); + _hostCandidates = ArrayPool.Shared.Rent(_capacity); + _hostMask = ArrayPool.Shared.Rent(_capacity); + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _candidateBuffer.Dispose(); + _maskBuffer.Dispose(); + ArrayPool.Shared.Return(_hostCandidates, clearArray: false); + ArrayPool.Shared.Return(_hostMask, clearArray: false); + _lease.Dispose(); + } + + internal void EvaluateCandidates( + ulong startDivisor, + ulong step, + ulong limit, + byte remainder10, + byte remainder8, + byte remainder5, + byte remainder3, + byte remainder7, + byte remainder11, + byte step10, + byte step8, + byte step5, + byte step3, + byte step7, + byte step11, + bool lastIsSeven, + Span candidates, + Span mask) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(MersenneNumberDivisorCandidateGpuEvaluator)); + } + + int count = candidates.Length; + if (count == 0) + { + return; + } + + EnsureCapacity(count); + + Span hostCandidateSpan = _hostCandidates.AsSpan(0, count); + Span hostMaskSpan = _hostMask.AsSpan(0, count); + + ArrayView1D candidateView = _candidateBuffer.View.SubView(0, count); + ArrayView1D maskView = _maskBuffer.View.SubView(0, count); + + var args = new CandidateFilterArgs( + startDivisor, + step, + limit, + remainder10, + remainder8, + remainder5, + remainder3, + remainder7, + remainder11, + step10, + step8, + step5, + step3, + step7, + step11, + lastIsSeven ? (byte)1 : (byte)0); + + _kernel(count, args, candidateView, maskView); + + candidateView.CopyToCPU(ref MemoryMarshal.GetReference(hostCandidateSpan), count); + maskView.CopyToCPU(ref MemoryMarshal.GetReference(hostMaskSpan), count); + + hostCandidateSpan.CopyTo(candidates); + hostMaskSpan.CopyTo(mask); + } + + private void EnsureCapacity(int required) + { + if (required <= _capacity) + { + return; + } + + int newCapacity = _capacity == 0 ? 1 : _capacity; + while (newCapacity < required) + { + newCapacity <<= 1; + } + + _candidateBuffer.Dispose(); + _maskBuffer.Dispose(); + ArrayPool.Shared.Return(_hostCandidates, clearArray: false); + ArrayPool.Shared.Return(_hostMask, clearArray: false); + + _capacity = newCapacity; + _candidateBuffer = _accelerator.Allocate1D(_capacity); + _maskBuffer = _accelerator.Allocate1D(_capacity); + _hostCandidates = ArrayPool.Shared.Rent(_capacity); + _hostMask = ArrayPool.Shared.Rent(_capacity); + } + + private static void EvaluateCandidatesKernel( + Index1D index, + CandidateFilterArgs args, + ArrayView candidates, + ArrayView mask) + { + int idx = index; + if (idx >= mask.Length) + { + return; + } + + ulong offset = (ulong)idx; + ulong candidate = args.StartDivisor + args.Step * offset; + candidates[idx] = candidate; + + if (candidate > args.Limit || (args.Step != 0UL && candidate < args.StartDivisor)) + { + mask[idx] = 0; + return; + } + + int remainderValue10 = ComputeRemainder(args.Remainder10, args.Step10, idx, 10); + int remainderValue8 = ComputeRemainder(args.Remainder8, args.Step8, idx, 8); + int remainderValue5 = ComputeRemainder(args.Remainder5, args.Step5, idx, 5); + int remainderValue3 = ComputeRemainder(args.Remainder3, args.Step3, idx, 3); + int remainderValue7 = ComputeRemainder(args.Remainder7, args.Step7, idx, 7); + int remainderValue11 = ComputeRemainder(args.Remainder11, args.Step11, idx, 11); + + bool lastIsSeven = args.LastIsSevenFlag != 0; + bool admissible = lastIsSeven + ? (remainderValue10 == 3 || remainderValue10 == 7 || remainderValue10 == 9) + : (remainderValue10 == 1 || remainderValue10 == 3 || remainderValue10 == 9); + + if (!admissible) + { + mask[idx] = 0; + return; + } + + if ((remainderValue8 != 1 && remainderValue8 != 7) + || remainderValue3 == 0 + || remainderValue5 == 0 + || remainderValue7 == 0 + || remainderValue11 == 0) + { + mask[idx] = 0; + return; + } + + mask[idx] = 1; + } + + private static int ComputeRemainder(byte remainder, byte step, int index, int modulus) + { + if (modulus == 0) + { + return 0; + } + + long value = remainder; + value += (long)step * index; + int result = (int)(value % modulus); + if (result < 0) + { + result += modulus; + } + + return result; + } + + private readonly struct CandidateFilterArgs + { + public CandidateFilterArgs( + ulong startDivisor, + ulong step, + ulong limit, + byte remainder10, + byte remainder8, + byte remainder5, + byte remainder3, + byte remainder7, + byte remainder11, + byte step10, + byte step8, + byte step5, + byte step3, + byte step7, + byte step11, + byte lastIsSevenFlag) + { + StartDivisor = startDivisor; + Step = step; + Limit = limit; + Remainder10 = remainder10; + Remainder8 = remainder8; + Remainder5 = remainder5; + Remainder3 = remainder3; + Remainder7 = remainder7; + Remainder11 = remainder11; + Step10 = step10; + Step8 = step8; + Step5 = step5; + Step3 = step3; + Step7 = step7; + Step11 = step11; + LastIsSevenFlag = lastIsSevenFlag; + } + + public ulong StartDivisor { get; } + public ulong Step { get; } + public ulong Limit { get; } + public byte Remainder10 { get; } + public byte Remainder8 { get; } + public byte Remainder5 { get; } + public byte Remainder3 { get; } + public byte Remainder7 { get; } + public byte Remainder11 { get; } + public byte Step10 { get; } + public byte Step8 { get; } + public byte Step5 { get; } + public byte Step3 { get; } + public byte Step7 { get; } + public byte Step11 { get; } + public byte LastIsSevenFlag { get; } + } +} From ea29035c554c7e02ab60cc32804119b9447aaab5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 01:35:40 +0200 Subject: [PATCH 6/8] Use GpuUInt128 in by-divisor GPU scan --- EvenPerfectBitScanner/Program.cs | 32 ++++----- ...MersenneNumberDivisorByDivisorCpuTester.cs | 67 ++++++++++--------- 2 files changed, 51 insertions(+), 48 deletions(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index e813d6b..03e11d9 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -579,9 +579,6 @@ internal static class Program Console.WriteLine("Initialization..."); // Compose a results file name that encodes configuration (before opening file) - string byDivisorDeltasLabel = mersenneOnGpu - ? "gpu" - : (byDivisorDeltasDevice == ByDivisorDeltasDevice.Gpu ? "cpu-gpu" : "cpu"); var builtName = BuildResultsFileName( useBitTransform, threadCount, @@ -591,19 +588,18 @@ internal static class Program useDivisor, useByDivisor, mersenneOnGpu, - useOrder, - useModuloWorkaround, - _useGcdFilter, - NttGpuMath.GpuTransformBackend, - gpuPrimeThreads, - sliceSize, - scanBatchSize, - _orderWarmupLimitOverride ?? 5_000_000UL, - NttGpuMath.ReductionMode, - mersenneOnGpu ? "gpu" : "cpu", - (GpuContextPool.ForceCpu ? "cpu" : "gpu"), - orderOnGpu ? "gpu" : "cpu", - byDivisorDeltasLabel); + useOrder, + useModuloWorkaround, + _useGcdFilter, + NttGpuMath.GpuTransformBackend, + gpuPrimeThreads, + sliceSize, + scanBatchSize, + _orderWarmupLimitOverride ?? 5_000_000UL, + NttGpuMath.ReductionMode, + mersenneOnGpu ? "gpu" : "cpu", + (GpuContextPool.ForceCpu ? "cpu" : "gpu"), + orderOnGpu ? "gpu" : "cpu"); if (!string.IsNullOrEmpty(_resultsPrefix)) { @@ -1020,7 +1016,7 @@ internal static class Program Console.WriteLine(" --help, -help, --?, -?, /? show this help message"); } - private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool useByDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice, string byDivisorDeltasDevice) + private static string BuildResultsFileName(bool bitInc, int threads, int block, GpuKernelType kernelType, bool useLucasFlag, bool useDivisorFlag, bool useByDivisorFlag, bool mersenneOnGpu, bool useOrder, bool useModWorkaround, bool useGcd, NttBackend nttBackend, int gpuPrimeThreads, int llSlice, int gpuScanBatch, ulong warmupLimit, ModReductionMode reduction, string mersenneDevice, string primesDevice, string orderDevice) { string inc = bitInc ? "bit" : "add"; string mers = useDivisorFlag @@ -1028,7 +1024,7 @@ internal static class Program : (useLucasFlag ? "lucas" : (useByDivisorFlag - ? $"bydivisor-{byDivisorDeltasDevice}" + ? "bydivisor" : (kernelType == GpuKernelType.Pow2Mod ? "pow2mod" : "incremental"))); string ntt = nttBackend == NttBackend.Staged ? "staged" : "reference"; string red = reduction switch { ModReductionMode.Mont64 => "mont64", ModReductionMode.Barrett128 => "barrett128", ModReductionMode.GpuUInt128 => "uint128", _ => "auto" }; diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 65f0a9c..41889aa 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -363,27 +363,24 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - UInt128 step128 = (UInt128)prime << 1; - if (step128 == UInt128.Zero) + GpuUInt128 step128 = new GpuUInt128(prime); + step128.ShiftLeft(1); + if (step128.IsZero) { processedAll = true; return false; } - if (step128 > ulong.MaxValue) - { - return CheckDivisorsCpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); - } - - UInt128 limit = allowedMax; - UInt128 divisor128 = step128 + UInt128.One; - if (divisor128 > limit) + GpuUInt128 limit128 = new GpuUInt128(allowedMax); + GpuUInt128 divisor128 = step128; + divisor128.Add(1UL); + if (divisor128.CompareTo(limit128) > 0) { processedAll = true; return false; } - ulong step = (ulong)step128; + ulong step = step128.Low; ulong limit64 = allowedMax; int batchCapacity = Math.Max(1, _batchSize); @@ -412,20 +409,27 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv byte step7 = (byte)(step % 7UL); byte step11 = (byte)(step % 11UL); - byte remainder10 = (byte)((ulong)divisor128 % 10UL); - byte remainder8 = (byte)((ulong)divisor128 % 8UL); - byte remainder5 = (byte)((ulong)divisor128 % 5UL); - byte remainder3 = (byte)((ulong)divisor128 % 3UL); - byte remainder7 = (byte)((ulong)divisor128 % 7UL); - byte remainder11 = (byte)((ulong)divisor128 % 11UL); + ulong startDivisor64 = divisor128.Low; + byte remainder10 = (byte)(startDivisor64 % 10UL); + byte remainder8 = (byte)(startDivisor64 % 8UL); + byte remainder5 = (byte)(startDivisor64 % 5UL); + byte remainder3 = (byte)(startDivisor64 % 3UL); + byte remainder7 = (byte)(startDivisor64 % 7UL); + byte remainder11 = (byte)(startDivisor64 % 11UL); bool lastIsSeven = (prime & 3UL) == 3UL; try { - while (divisor128 <= limit) + while (divisor128.CompareTo(limit128) <= 0) { - int chunkCount = ComputeChunkCount(divisor128, limit, step128, batchCapacity); + if (divisor128.High != 0UL) + { + break; + } + + ulong currentDivisor = divisor128.Low; + int chunkCount = ComputeChunkCount(currentDivisor, limit64, step, batchCapacity); if (chunkCount <= 0) { break; @@ -435,7 +439,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv Span maskSpan = maskBuffer.AsSpan(0, chunkCount); candidateEvaluator.EvaluateCandidates( - (ulong)divisor128, + currentDivisor, step, limit64, remainder10, @@ -508,7 +512,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } - divisor128 += step128 * (UInt128)chunkCount; + GpuUInt128 increment = step128; + increment.Mul((ulong)chunkCount); + divisor128.Add(increment); + remainder10 = AdvanceRemainder(remainder10, step10, 10, chunkCount); remainder8 = AdvanceRemainder(remainder8, step8, 8, chunkCount); remainder5 = AdvanceRemainder(remainder5, step5, 5, chunkCount); @@ -524,26 +531,26 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv candidateEvaluator.Dispose(); } - processedAll = divisor128 > limit; + processedAll = divisor128.CompareTo(limit128) > 0; return false; } - private static int ComputeChunkCount(UInt128 currentDivisor, UInt128 limit, UInt128 step, int maxCount) + private static int ComputeChunkCount(ulong currentDivisor, ulong limit, ulong step, int maxCount) { - if (maxCount <= 0 || step == UInt128.Zero || currentDivisor > limit) + if (maxCount <= 0 || step == 0UL || currentDivisor > limit) { return 0; } - UInt128 distance = limit - currentDivisor; - UInt128 stepsAvailable = distance / step; - UInt128 total = stepsAvailable + UInt128.One; - if (total > (UInt128)maxCount) + ulong distance = limit - currentDivisor; + ulong stepsAvailable = step == 0UL ? 0UL : distance / step; + ulong total = stepsAvailable + 1UL; + if (total > (ulong)maxCount) { - total = (UInt128)maxCount; + total = (ulong)maxCount; } - if (total > (UInt128)int.MaxValue) + if (total > (ulong)int.MaxValue) { return int.MaxValue; } From 483584b5d4b687ba75ac788c9738e40e0590853c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 01:35:46 +0200 Subject: [PATCH 7/8] Add GPU options for by-divisor Montgomery data --- EvenPerfectBitScanner/Program.cs | 20 ++ .../ByDivisorMontgomeryDevice.cs | 7 + ...MersenneNumberDivisorByDivisorCpuTester.cs | 185 ++++++++++-------- ...rsenneNumberDivisorMontgomeryGpuBuilder.cs | 179 +++++++++++++++++ ...ersenneNumberDivisorRemainderGpuStepper.cs | 135 +++++++++++++ 5 files changed, 449 insertions(+), 77 deletions(-) create mode 100644 PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs create mode 100644 PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs create mode 100644 PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 03e11d9..06a1bf1 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -76,6 +76,7 @@ internal static class Program bool useDivisor = false; // M_p divisibility by specific divisor bool useByDivisor = false; // Iterative divisor scan across primes UInt128 divisor = UInt128.Zero; ByDivisorDeltasDevice byDivisorDeltasDevice = ByDivisorDeltasDevice.Cpu; + ByDivisorMontgomeryDevice byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Cpu; // Device routing bool useGpuCycles = true; bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device @@ -200,6 +201,23 @@ internal static class Program return; } } + else if (arg.StartsWith("--bydivisor-montgomery-device=", StringComparison.OrdinalIgnoreCase)) + { + ReadOnlySpan deviceValue = arg.AsSpan(arg.IndexOf('=') + 1); + if (deviceValue.Equals("gpu", StringComparison.OrdinalIgnoreCase)) + { + byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Gpu; + } + else if (deviceValue.Equals("cpu", StringComparison.OrdinalIgnoreCase)) + { + byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Cpu; + } + else + { + Console.WriteLine("Invalid --bydivisor-montgomery-device value. Expected cpu or gpu."); + return; + } + } else if (arg.StartsWith("--prime-test-limit=", StringComparison.OrdinalIgnoreCase)) { ReadOnlySpan value = arg.AsSpan(arg.IndexOf('=') + 1); @@ -552,6 +570,7 @@ internal static class Program if (_byDivisorTester is MersenneNumberDivisorByDivisorCpuTester cpuByDivisorTester) { cpuByDivisorTester.DeltasDevice = byDivisorDeltasDevice; + cpuByDivisorTester.MontgomeryDevice = byDivisorMontgomeryDevice; } } @@ -1007,6 +1026,7 @@ internal static class Program Console.WriteLine(" --divisor-cycles-limit= cycle search iterations when --mersenne=divisor"); Console.WriteLine(" --prime-test-limit= Time limit for --mersenne=bydivisor prime checks (e.g. 5s, 500ms)"); Console.WriteLine(" --bydivisor-deltas-device=cpu|gpu Device for --mersenne=bydivisor residue batches (default: cpu)"); + Console.WriteLine(" --bydivisor-montgomery-device=cpu|gpu Device for Montgomery data when using --mersenne=bydivisor (default: cpu)"); Console.WriteLine(" --use-order test primality via q order"); Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU"); // mod-automaton removed diff --git a/PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs b/PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs new file mode 100644 index 0000000..13d7452 --- /dev/null +++ b/PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs @@ -0,0 +1,7 @@ +namespace PerfectNumbers.Core; + +public enum ByDivisorMontgomeryDevice +{ + Cpu = 0, + Gpu = 1, +} diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 41889aa..53242d9 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -19,7 +19,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private ulong _configuredMaxPrime; private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu; - + private ByDivisorMontgomeryDevice _montgomeryDevice = ByDivisorMontgomeryDevice.Cpu; public int BatchSize { @@ -43,6 +43,18 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } + public ByDivisorMontgomeryDevice MontgomeryDevice + { + get => _montgomeryDevice; + set + { + lock (_sync) + { + _montgomeryDevice = value; + } + } + } + public void ConfigureFromMaxPrime(ulong maxPrime) { lock (_sync) @@ -402,26 +414,51 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv MontgomeryDivisorData cachedDivisorData = default; bool hasCachedDivisorData = false; - byte step10 = (byte)(step % 10UL); - byte step8 = (byte)(step % 8UL); - byte step5 = (byte)(step % 5UL); - byte step3 = (byte)(step % 3UL); - byte step7 = (byte)(step % 7UL); - byte step11 = (byte)(step % 11UL); + const int RemainderTableSize = 6; + byte[] remainderBuffer = ArrayPool.Shared.Rent(RemainderTableSize); + byte[] stepBuffer = ArrayPool.Shared.Rent(RemainderTableSize); + Span remainderSpan = remainderBuffer.AsSpan(0, RemainderTableSize); + Span stepSpan = stepBuffer.AsSpan(0, RemainderTableSize); + Span modulusSpan = stackalloc byte[RemainderTableSize] + { + 10, + 8, + 5, + 3, + 7, + 11, + }; ulong startDivisor64 = divisor128.Low; - byte remainder10 = (byte)(startDivisor64 % 10UL); - byte remainder8 = (byte)(startDivisor64 % 8UL); - byte remainder5 = (byte)(startDivisor64 % 5UL); - byte remainder3 = (byte)(startDivisor64 % 3UL); - byte remainder7 = (byte)(startDivisor64 % 7UL); - byte remainder11 = (byte)(startDivisor64 % 11UL); + remainderSpan[0] = (byte)(startDivisor64 % 10UL); + remainderSpan[1] = (byte)(startDivisor64 % 8UL); + remainderSpan[2] = (byte)(startDivisor64 % 5UL); + remainderSpan[3] = (byte)(startDivisor64 % 3UL); + remainderSpan[4] = (byte)(startDivisor64 % 7UL); + remainderSpan[5] = (byte)(startDivisor64 % 11UL); + + stepSpan[0] = (byte)(step % 10UL); + stepSpan[1] = (byte)(step % 8UL); + stepSpan[2] = (byte)(step % 5UL); + stepSpan[3] = (byte)(step % 3UL); + stepSpan[4] = (byte)(step % 7UL); + stepSpan[5] = (byte)(step % 11UL); bool lastIsSeven = (prime & 3UL) == 3UL; + MersenneNumberDivisorRemainderGpuStepper remainderStepper = new(stepSpan, modulusSpan); + MersenneNumberDivisorMontgomeryGpuBuilder? montgomeryBuilder = _montgomeryDevice == ByDivisorMontgomeryDevice.Gpu + ? new MersenneNumberDivisorMontgomeryGpuBuilder(batchCapacity) + : null; + MontgomeryDivisorData[]? montgomeryBuffer = montgomeryBuilder != null + ? ArrayPool.Shared.Rent(batchCapacity) + : null; + + ulong remainingCandidates = ComputeRemainingCandidateCount(startDivisor64, limit64, step); + try { - while (divisor128.CompareTo(limit128) <= 0) + while (remainingCandidates > 0UL && divisor128.CompareTo(limit128) <= 0) { if (divisor128.High != 0UL) { @@ -429,7 +466,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } ulong currentDivisor = divisor128.Low; - int chunkCount = ComputeChunkCount(currentDivisor, limit64, step, batchCapacity); + int chunkCount = remainingCandidates > (ulong)batchCapacity + ? batchCapacity + : (int)remainingCandidates; if (chunkCount <= 0) { break; @@ -442,22 +481,29 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv currentDivisor, step, limit64, - remainder10, - remainder8, - remainder5, - remainder3, - remainder7, - remainder11, - step10, - step8, - step5, - step3, - step7, - step11, + remainderSpan[0], + remainderSpan[1], + remainderSpan[2], + remainderSpan[3], + remainderSpan[4], + remainderSpan[5], + stepSpan[0], + stepSpan[1], + stepSpan[2], + stepSpan[3], + stepSpan[4], + stepSpan[5], lastIsSeven, candidateSpan, maskSpan); + Span montgomerySpan = default; + if (montgomeryBuilder != null && montgomeryBuffer != null) + { + montgomerySpan = montgomeryBuffer.AsSpan(0, chunkCount); + montgomeryBuilder.Build(candidateSpan, montgomerySpan); + } + for (int i = 0; i < chunkCount; i++) { if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) @@ -477,16 +523,28 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv continue; } - if (!hasCachedDivisorData || cachedModulus != candidate) + MontgomeryDivisorData divisorData; + if (!montgomerySpan.IsEmpty) { - cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate); - cachedModulus = candidate; - hasCachedDivisorData = true; + divisorData = montgomerySpan[i]; + if (divisorData.Modulus != candidate) + { + divisorData = MontgomeryDivisorData.FromModulus(candidate); + } + } + else + { + if (!hasCachedDivisorData || cachedModulus != candidate) + { + cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate); + cachedModulus = candidate; + hasCachedDivisorData = true; + } + + divisorData = cachedDivisorData; } - MontgomeryDivisorData divisorData = cachedDivisorData; ulong divisorCycle; - if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) { divisorCycle = cycleCache.GetCycleLength(candidate); @@ -516,67 +574,40 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv increment.Mul((ulong)chunkCount); divisor128.Add(increment); - remainder10 = AdvanceRemainder(remainder10, step10, 10, chunkCount); - remainder8 = AdvanceRemainder(remainder8, step8, 8, chunkCount); - remainder5 = AdvanceRemainder(remainder5, step5, 5, chunkCount); - remainder3 = AdvanceRemainder(remainder3, step3, 3, chunkCount); - remainder7 = AdvanceRemainder(remainder7, step7, 7, chunkCount); - remainder11 = AdvanceRemainder(remainder11, step11, 11, chunkCount); + remainderStepper.Advance(chunkCount, remainderSpan); + remainingCandidates -= (ulong)chunkCount; } } finally { + montgomeryBuilder?.Dispose(); + if (montgomeryBuffer != null) + { + ArrayPool.Shared.Return(montgomeryBuffer, clearArray: false); + } + + remainderStepper.Dispose(); + ArrayPool.Shared.Return(remainderBuffer, clearArray: false); + ArrayPool.Shared.Return(stepBuffer, clearArray: false); ArrayPool.Shared.Return(candidateBuffer, clearArray: false); ArrayPool.Shared.Return(maskBuffer, clearArray: false); candidateEvaluator.Dispose(); } - processedAll = divisor128.CompareTo(limit128) > 0; + processedAll = remainingCandidates == 0UL && divisor128.CompareTo(limit128) > 0; return false; } - private static int ComputeChunkCount(ulong currentDivisor, ulong limit, ulong step, int maxCount) + private static ulong ComputeRemainingCandidateCount(ulong currentDivisor, ulong limit, ulong step) { - if (maxCount <= 0 || step == 0UL || currentDivisor > limit) + if (step == 0UL || currentDivisor > limit) { - return 0; + return 0UL; } ulong distance = limit - currentDivisor; - ulong stepsAvailable = step == 0UL ? 0UL : distance / step; - ulong total = stepsAvailable + 1UL; - if (total > (ulong)maxCount) - { - total = (ulong)maxCount; - } - - if (total > (ulong)int.MaxValue) - { - return int.MaxValue; - } - - return (int)total; - } - - private static byte AdvanceRemainder(byte remainder, byte step, byte modulus, int count) - { - if (modulus == 0 || count <= 0) - { - return remainder; - } - - int modulusValue = modulus; - int remainderValue = remainder % modulusValue; - int stepValue = step % modulusValue; - - long delta = (long)stepValue * count; - int adjusted = (int)((remainderValue + delta) % modulusValue); - if (adjusted < 0) - { - adjusted += modulusValue; - } - - return (byte)adjusted; + ulong stepsAvailable = distance / step; + return stepsAvailable + 1UL; } private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs new file mode 100644 index 0000000..7c944cc --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs @@ -0,0 +1,179 @@ +using System; +using System.Buffers; +using System.Runtime.InteropServices; +using PerfectNumbers.Core; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable +{ + private readonly GpuContextPool.GpuContextLease _lease; + private readonly Accelerator _accelerator; + private readonly Action, ArrayView> _kernel; + private MemoryBuffer1D _modulusBuffer; + private MemoryBuffer1D _resultBuffer; + private ulong[] _hostModulusBuffer; + private MontgomeryDivisorData[] _hostResultBuffer; + private int _capacity; + private bool _disposed; + + internal MersenneNumberDivisorMontgomeryGpuBuilder(int batchSize) + { + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView>(BuildKernel); + _capacity = Math.Max(1, batchSize); + _modulusBuffer = _accelerator.Allocate1D(_capacity); + _resultBuffer = _accelerator.Allocate1D(_capacity); + _hostModulusBuffer = ArrayPool.Shared.Rent(_capacity); + _hostResultBuffer = ArrayPool.Shared.Rent(_capacity); + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _modulusBuffer.Dispose(); + _resultBuffer.Dispose(); + ArrayPool.Shared.Return(_hostModulusBuffer, clearArray: false); + ArrayPool.Shared.Return(_hostResultBuffer, clearArray: false); + _lease.Dispose(); + } + + internal void Build(ReadOnlySpan moduli, Span destination) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(MersenneNumberDivisorMontgomeryGpuBuilder)); + } + + int length = moduli.Length; + if (length == 0) + { + return; + } + + EnsureCapacity(length); + + Span hostModulusSpan = _hostModulusBuffer.AsSpan(0, length); + moduli.CopyTo(hostModulusSpan); + + ArrayView1D modulusView = _modulusBuffer.View.SubView(0, length); + modulusView.CopyFromCPU(ref MemoryMarshal.GetReference(hostModulusSpan), length); + + ArrayView1D resultView = _resultBuffer.View.SubView(0, length); + _kernel(length, modulusView, resultView); + + Span hostResultSpan = _hostResultBuffer.AsSpan(0, length); + resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostResultSpan), length); + + hostResultSpan.CopyTo(destination); + } + + private void EnsureCapacity(int required) + { + if (required <= _capacity) + { + return; + } + + int newCapacity = _capacity == 0 ? 1 : _capacity; + while (newCapacity < required) + { + newCapacity <<= 1; + } + + _modulusBuffer.Dispose(); + _resultBuffer.Dispose(); + ArrayPool.Shared.Return(_hostModulusBuffer, clearArray: false); + ArrayPool.Shared.Return(_hostResultBuffer, clearArray: false); + + _capacity = newCapacity; + _modulusBuffer = _accelerator.Allocate1D(_capacity); + _resultBuffer = _accelerator.Allocate1D(_capacity); + _hostModulusBuffer = ArrayPool.Shared.Rent(_capacity); + _hostResultBuffer = ArrayPool.Shared.Rent(_capacity); + } + + private static void BuildKernel(Index1D index, ArrayView moduli, ArrayView destination) + { + int idx = index; + if (idx >= moduli.Length) + { + return; + } + + ulong modulus = moduli[idx]; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + destination[idx] = default; + return; + } + + ulong nPrime = ComputeMontgomeryNPrime(modulus); + ulong montgomeryOne = ComputeMontgomeryResidue64(modulus); + ulong montgomeryTwo = montgomeryOne.MulMod64GpuCompatible(2UL, modulus); + ulong montgomeryTwoSquared = MontgomeryMultiplyGpu(montgomeryTwo, montgomeryTwo, modulus, nPrime); + + destination[idx] = new MontgomeryDivisorData( + modulus, + nPrime, + montgomeryOne, + montgomeryTwo, + montgomeryTwoSquared); + } + + private static ulong ComputeMontgomeryResidue64(ulong modulus) + { + ulong remainder = ulong.MaxValue % modulus; + remainder++; + if (remainder >= modulus) + { + remainder -= modulus; + } + + return remainder; + } + + private static ulong MontgomeryMultiplyGpu(ulong a, ulong b, ulong modulus, ulong nPrime) + { + GpuUInt128 product = new(a); + product.Mul64(new GpuUInt128(b)); + ulong tLow = product.Low; + ulong tHigh = product.High; + + ulong m = unchecked(tLow * nPrime); + GpuUInt128 mProduct = new(m); + mProduct.Mul64(new GpuUInt128(modulus)); + ulong mLow = mProduct.Low; + ulong mHigh = mProduct.High; + + ulong carry = unchecked((tLow + mLow) < tLow ? 1UL : 0UL); + ulong result = tHigh + mHigh + carry; + if (result >= modulus) + { + result -= modulus; + } + + return result; + } + + private static ulong ComputeMontgomeryNPrime(ulong modulus) + { + ulong inv = modulus; + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + return unchecked(0UL - inv); + } +} diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs new file mode 100644 index 0000000..e110be9 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs @@ -0,0 +1,135 @@ +using System; +using System.Buffers; +using System.Runtime.InteropServices; +using ILGPU; +using ILGPU.Runtime; + +namespace PerfectNumbers.Core.Gpu; + +internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable +{ + private readonly GpuContextPool.GpuContextLease _lease; + private readonly Accelerator _accelerator; + private readonly Action, ArrayView, ArrayView, ArrayView> _kernel; + private MemoryBuffer1D _remainderBuffer; + private MemoryBuffer1D _stepBuffer; + private MemoryBuffer1D _modulusBuffer; + private MemoryBuffer1D _resultBuffer; + private byte[] _hostRemainders; + private byte[] _hostSteps; + private byte[] _hostModuli; + private readonly int _length; + private int _capacity; + private bool _disposed; + + internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan steps, ReadOnlySpan moduli) + { + if (steps.Length != moduli.Length) + { + throw new ArgumentException("Steps and moduli must have the same length."); + } + + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(AdvanceRemaindersKernel); + + _length = steps.Length; + _capacity = Math.Max(1, _length); + _remainderBuffer = _accelerator.Allocate1D(_capacity); + _stepBuffer = _accelerator.Allocate1D(_capacity); + _modulusBuffer = _accelerator.Allocate1D(_capacity); + _resultBuffer = _accelerator.Allocate1D(_capacity); + + _hostRemainders = ArrayPool.Shared.Rent(_capacity); + _hostSteps = ArrayPool.Shared.Rent(_capacity); + _hostModuli = ArrayPool.Shared.Rent(_capacity); + + steps.CopyTo(_hostSteps.AsSpan(0, _length)); + moduli.CopyTo(_hostModuli.AsSpan(0, _length)); + + ArrayView1D stepView = _stepBuffer.View.SubView(0, _length); + stepView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostSteps.AsSpan(0, _length)), _length); + + ArrayView1D modulusView = _modulusBuffer.View.SubView(0, _length); + modulusView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostModuli.AsSpan(0, _length)), _length); + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _remainderBuffer.Dispose(); + _stepBuffer.Dispose(); + _modulusBuffer.Dispose(); + _resultBuffer.Dispose(); + ArrayPool.Shared.Return(_hostRemainders, clearArray: false); + ArrayPool.Shared.Return(_hostSteps, clearArray: false); + ArrayPool.Shared.Return(_hostModuli, clearArray: false); + _lease.Dispose(); + } + + internal void Advance(int advanceCount, Span remainders) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(MersenneNumberDivisorRemainderGpuStepper)); + } + + if (advanceCount <= 0 || remainders.Length == 0) + { + return; + } + + if (remainders.Length != _length) + { + throw new ArgumentException("Remainder buffer length must match the configured length."); + } + + Span hostRemainderSpan = _hostRemainders.AsSpan(0, _length); + remainders.CopyTo(hostRemainderSpan); + + ArrayView1D remainderView = _remainderBuffer.View.SubView(0, _length); + remainderView.CopyFromCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length); + + ArrayView1D stepView = _stepBuffer.View.SubView(0, _length); + ArrayView1D modulusView = _modulusBuffer.View.SubView(0, _length); + ArrayView1D resultView = _resultBuffer.View.SubView(0, _length); + + _kernel(_length, advanceCount, remainderView, stepView, modulusView, resultView); + + resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length); + hostRemainderSpan.CopyTo(remainders); + } + + private static void AdvanceRemaindersKernel(Index1D index, int advanceCount, ArrayView remainders, ArrayView steps, ArrayView moduli, ArrayView results) + { + int idx = index; + if (idx >= remainders.Length) + { + return; + } + + int modulus = moduli[idx]; + if (modulus == 0) + { + results[idx] = remainders[idx]; + return; + } + + int remainder = remainders[idx] % modulus; + int step = steps[idx] % modulus; + + long delta = (long)step * advanceCount; + int adjusted = (int)((remainder + delta) % modulus); + if (adjusted < 0) + { + adjusted += modulus; + } + + results[idx] = (byte)adjusted; + } +} From 9e06a21718cc707f63f720cc66bbd4040510faa4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Mon, 13 Oct 2025 01:35:50 +0200 Subject: [PATCH 8/8] Streamline GPU remainder batching tables --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 26 ++--------- ...senneNumberDivisorCandidateGpuEvaluator.cs | 45 +++++++++---------- ...ersenneNumberDivisorRemainderGpuStepper.cs | 43 ++++++++++-------- 3 files changed, 49 insertions(+), 65 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 53242d9..0dfcc24 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -419,15 +419,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv byte[] stepBuffer = ArrayPool.Shared.Rent(RemainderTableSize); Span remainderSpan = remainderBuffer.AsSpan(0, RemainderTableSize); Span stepSpan = stepBuffer.AsSpan(0, RemainderTableSize); - Span modulusSpan = stackalloc byte[RemainderTableSize] - { - 10, - 8, - 5, - 3, - 7, - 11, - }; + // The remainder and step tables follow the fixed modulus order: 10, 8, 5, 3, 7, 11. ulong startDivisor64 = divisor128.Low; remainderSpan[0] = (byte)(startDivisor64 % 10UL); @@ -446,7 +438,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv bool lastIsSeven = (prime & 3UL) == 3UL; - MersenneNumberDivisorRemainderGpuStepper remainderStepper = new(stepSpan, modulusSpan); + MersenneNumberDivisorRemainderGpuStepper remainderStepper = new(stepSpan); MersenneNumberDivisorMontgomeryGpuBuilder? montgomeryBuilder = _montgomeryDevice == ByDivisorMontgomeryDevice.Gpu ? new MersenneNumberDivisorMontgomeryGpuBuilder(batchCapacity) : null; @@ -481,18 +473,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv currentDivisor, step, limit64, - remainderSpan[0], - remainderSpan[1], - remainderSpan[2], - remainderSpan[3], - remainderSpan[4], - remainderSpan[5], - stepSpan[0], - stepSpan[1], - stepSpan[2], - stepSpan[3], - stepSpan[4], - stepSpan[5], + remainderSpan, + stepSpan, lastIsSeven, candidateSpan, maskSpan); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs index 9a75acb..e4e0fe1 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs @@ -8,6 +8,8 @@ namespace PerfectNumbers.Core.Gpu; internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable { + private const int RemainderCount = 6; + private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; private readonly Action, ArrayView> _kernel; @@ -49,18 +51,8 @@ internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable ulong startDivisor, ulong step, ulong limit, - byte remainder10, - byte remainder8, - byte remainder5, - byte remainder3, - byte remainder7, - byte remainder11, - byte step10, - byte step8, - byte step5, - byte step3, - byte step7, - byte step11, + ReadOnlySpan remainders, + ReadOnlySpan steps, bool lastIsSeven, Span candidates, Span mask) @@ -70,6 +62,11 @@ internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable throw new ObjectDisposedException(nameof(MersenneNumberDivisorCandidateGpuEvaluator)); } + if (remainders.Length != RemainderCount || steps.Length != RemainderCount) + { + throw new ArgumentException("Remainder and step tables must contain exactly six entries."); + } + int count = candidates.Length; if (count == 0) { @@ -88,18 +85,18 @@ internal sealed class MersenneNumberDivisorCandidateGpuEvaluator : IDisposable startDivisor, step, limit, - remainder10, - remainder8, - remainder5, - remainder3, - remainder7, - remainder11, - step10, - step8, - step5, - step3, - step7, - step11, + remainders[0], + remainders[1], + remainders[2], + remainders[3], + remainders[4], + remainders[5], + steps[0], + steps[1], + steps[2], + steps[3], + steps[4], + steps[5], lastIsSeven ? (byte)1 : (byte)0); _kernel(count, args, candidateView, maskView); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs index e110be9..e293430 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs @@ -8,50 +8,44 @@ namespace PerfectNumbers.Core.Gpu; internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable { + private const int RemainderCount = 6; + private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; - private readonly Action, ArrayView, ArrayView, ArrayView> _kernel; + private readonly Action, ArrayView, ArrayView> _kernel; private MemoryBuffer1D _remainderBuffer; private MemoryBuffer1D _stepBuffer; - private MemoryBuffer1D _modulusBuffer; private MemoryBuffer1D _resultBuffer; private byte[] _hostRemainders; private byte[] _hostSteps; - private byte[] _hostModuli; private readonly int _length; private int _capacity; private bool _disposed; - internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan steps, ReadOnlySpan moduli) + internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan steps) { - if (steps.Length != moduli.Length) + if (steps.Length != RemainderCount) { - throw new ArgumentException("Steps and moduli must have the same length."); + throw new ArgumentException("Step table must contain exactly six entries."); } _lease = GpuContextPool.RentPreferred(preferCpu: false); _accelerator = _lease.Accelerator; - _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView>(AdvanceRemaindersKernel); + _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(AdvanceRemaindersKernel); _length = steps.Length; _capacity = Math.Max(1, _length); _remainderBuffer = _accelerator.Allocate1D(_capacity); _stepBuffer = _accelerator.Allocate1D(_capacity); - _modulusBuffer = _accelerator.Allocate1D(_capacity); _resultBuffer = _accelerator.Allocate1D(_capacity); _hostRemainders = ArrayPool.Shared.Rent(_capacity); _hostSteps = ArrayPool.Shared.Rent(_capacity); - _hostModuli = ArrayPool.Shared.Rent(_capacity); steps.CopyTo(_hostSteps.AsSpan(0, _length)); - moduli.CopyTo(_hostModuli.AsSpan(0, _length)); ArrayView1D stepView = _stepBuffer.View.SubView(0, _length); stepView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostSteps.AsSpan(0, _length)), _length); - - ArrayView1D modulusView = _modulusBuffer.View.SubView(0, _length); - modulusView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostModuli.AsSpan(0, _length)), _length); } public void Dispose() @@ -64,11 +58,9 @@ internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable _disposed = true; _remainderBuffer.Dispose(); _stepBuffer.Dispose(); - _modulusBuffer.Dispose(); _resultBuffer.Dispose(); ArrayPool.Shared.Return(_hostRemainders, clearArray: false); ArrayPool.Shared.Return(_hostSteps, clearArray: false); - ArrayPool.Shared.Return(_hostModuli, clearArray: false); _lease.Dispose(); } @@ -96,16 +88,15 @@ internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable remainderView.CopyFromCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length); ArrayView1D stepView = _stepBuffer.View.SubView(0, _length); - ArrayView1D modulusView = _modulusBuffer.View.SubView(0, _length); ArrayView1D resultView = _resultBuffer.View.SubView(0, _length); - _kernel(_length, advanceCount, remainderView, stepView, modulusView, resultView); + _kernel(_length, advanceCount, remainderView, stepView, resultView); resultView.CopyToCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length); hostRemainderSpan.CopyTo(remainders); } - private static void AdvanceRemaindersKernel(Index1D index, int advanceCount, ArrayView remainders, ArrayView steps, ArrayView moduli, ArrayView results) + private static void AdvanceRemaindersKernel(Index1D index, int advanceCount, ArrayView remainders, ArrayView steps, ArrayView results) { int idx = index; if (idx >= remainders.Length) @@ -113,7 +104,7 @@ internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable return; } - int modulus = moduli[idx]; + int modulus = GetModulus(idx); if (modulus == 0) { results[idx] = remainders[idx]; @@ -132,4 +123,18 @@ internal sealed class MersenneNumberDivisorRemainderGpuStepper : IDisposable results[idx] = (byte)adjusted; } + + private static int GetModulus(int index) + { + return index switch + { + 0 => 10, + 1 => 8, + 2 => 5, + 3 => 3, + 4 => 7, + 5 => 11, + _ => 0, + }; + } }