diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 7cefd12..06a1bf1 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -75,6 +75,8 @@ 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; + ByDivisorMontgomeryDevice byDivisorMontgomeryDevice = ByDivisorMontgomeryDevice.Cpu; // Device routing bool useGpuCycles = true; bool mersenneOnGpu = true; // controls Lucas/incremental/pow2mod device @@ -182,6 +184,40 @@ 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("--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); @@ -531,6 +567,11 @@ internal static class Program ? new MersenneNumberDivisorByDivisorGpuTester() : new MersenneNumberDivisorByDivisorCpuTester(); _byDivisorTester.BatchSize = scanBatchSize; + if (_byDivisorTester is MersenneNumberDivisorByDivisorCpuTester cpuByDivisorTester) + { + cpuByDivisorTester.DeltasDevice = byDivisorDeltasDevice; + cpuByDivisorTester.MontgomeryDevice = byDivisorMontgomeryDevice; + } } // Load RLE blacklist (optional) @@ -566,18 +607,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"); + 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)) { @@ -984,6 +1025,8 @@ 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(" --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/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/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 d269982..0dfcc24 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,8 +1,10 @@ using System; +using System.Buffers; using System.Collections.Concurrent; using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Cpu; @@ -15,17 +17,49 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool _isConfigured; private int _batchSize = 1_024; + private ulong _configuredMaxPrime; + private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu; + private ByDivisorMontgomeryDevice _montgomeryDevice = ByDivisorMontgomeryDevice.Cpu; public int BatchSize { get => _batchSize; - set => _batchSize = Math.Max(1, value); + set + { + int sanitized = Math.Max(1, value); + _batchSize = sanitized; + } + } + + public ByDivisorDeltasDevice DeltasDevice + { + get => _deltasDevice; + set + { + lock (_sync) + { + _deltasDevice = value; + } + } + } + + public ByDivisorMontgomeryDevice MontgomeryDevice + { + get => _montgomeryDevice; + set + { + lock (_sync) + { + _montgomeryDevice = value; + } + } } public void ConfigureFromMaxPrime(ulong maxPrime) { lock (_sync) { + _configuredMaxPrime = maxPrime; _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); _lastStatusDivisor = 0UL; _isConfigured = true; @@ -166,6 +200,22 @@ 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; @@ -298,6 +348,250 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } + + private bool CheckDivisorsGpu( + 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; + } + + GpuUInt128 step128 = new GpuUInt128(prime); + step128.ShiftLeft(1); + if (step128.IsZero) + { + processedAll = true; + return false; + } + + GpuUInt128 limit128 = new GpuUInt128(allowedMax); + GpuUInt128 divisor128 = step128; + divisor128.Add(1UL); + if (divisor128.CompareTo(limit128) > 0) + { + processedAll = true; + return false; + } + + ulong step = step128.Low; + ulong limit64 = allowedMax; + + int batchCapacity = Math.Max(1, _batchSize); + DivisorCycleCache cycleCache = DivisorCycleCache.Shared; + int preferredCycleBatch = cycleCache.PreferredBatchSize; + if (preferredCycleBatch > 0 && preferredCycleBatch < batchCapacity) + { + batchCapacity = preferredCycleBatch; + } + + var candidateEvaluator = new MersenneNumberDivisorCandidateGpuEvaluator(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; + + 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); + // 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); + 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); + 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 (remainingCandidates > 0UL && divisor128.CompareTo(limit128) <= 0) + { + if (divisor128.High != 0UL) + { + break; + } + + ulong currentDivisor = divisor128.Low; + int chunkCount = remainingCandidates > (ulong)batchCapacity + ? batchCapacity + : (int)remainingCandidates; + if (chunkCount <= 0) + { + break; + } + + Span candidateSpan = candidateBuffer.AsSpan(0, chunkCount); + Span maskSpan = maskBuffer.AsSpan(0, chunkCount); + + candidateEvaluator.EvaluateCandidates( + currentDivisor, + step, + limit64, + remainderSpan, + stepSpan, + 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()) + { + processedAll = false; + return false; + } + + iterationCounter++; + + ulong candidate = candidateSpan[i]; + processedCount++; + lastProcessed = candidate; + + if (maskSpan[i] == 0) + { + continue; + } + + MontgomeryDivisorData divisorData; + if (!montgomerySpan.IsEmpty) + { + 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; + } + + ulong divisorCycle; + if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) + { + divisorCycle = cycleCache.GetCycleLength(candidate); + } + else + { + factorCache ??= new Dictionary(8); + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + { + divisorCycle = cycleCache.GetCycleLength(candidate); + } + } + + if (divisorCycle == prime) + { + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + GpuUInt128 increment = step128; + increment.Mul((ulong)chunkCount); + divisor128.Add(increment); + + 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 = remainingCandidates == 0UL && divisor128.CompareTo(limit128) > 0; + return false; + } + + private static ulong ComputeRemainingCandidateCount(ulong currentDivisor, ulong limit, ulong step) + { + if (step == 0UL || currentDivisor > limit) + { + return 0UL; + } + + ulong distance = limit - currentDivisor; + ulong stepsAvailable = distance / step; + return stepsAvailable + 1UL; + } + private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) { ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData); @@ -372,6 +666,11 @@ 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; + private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator; internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) { @@ -412,6 +711,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,9 +728,24 @@ 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(cachedData); + var exponentStepper = new ExponentRemainderStepper(divisorData); if (!exponentStepper.IsValidModulus) { hits.Clear(); @@ -438,6 +759,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv ? (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]); @@ -451,6 +773,172 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } + private void CheckDivisorGpu( + in MontgomeryDivisorData divisorData, + ulong divisorCycle, + ReadOnlySpan primes, + Span hits, + int length) + { + EnsureCapacity(length); + + Span deltaSpan = _primeDeltas!.AsSpan(0, length); + ComputePrimeDeltas(primes, deltaSpan); + + Span remainderSpan = _cycleRemainders!.AsSpan(0, length); + ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); + + EnsureGpuEvaluator(); + Span residueSpan = _montgomeryDeltas!.AsSpan(0, length); + _gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan); + + ulong montgomeryOne = divisorData.MontgomeryOne; + for (int i = 0; i < length; i++) + { + hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0; + } + } + + private void EnsureGpuEvaluator() + { + if (_gpuResidueEvaluator is null) + { + _gpuResidueEvaluator = new MersenneNumberDivisorResidueGpuEvaluator(_owner._batchSize); + } + } + + private void ReleaseGpuEvaluator() + { + if (_gpuResidueEvaluator is null) + { + return; + } + + _gpuResidueEvaluator.Dispose(); + _gpuResidueEvaluator = null; + } + + 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) + { + throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing."); + } + + 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 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) @@ -459,6 +947,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } _disposed = true; + ReleaseGpuEvaluator(); _owner.ReturnSession(this); } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs new file mode 100644 index 0000000..e4e0fe1 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorCandidateGpuEvaluator.cs @@ -0,0 +1,262 @@ +using System; +using System.Buffers; +using System.Runtime.InteropServices; +using ILGPU; +using ILGPU.Runtime; + +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; + 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, + ReadOnlySpan remainders, + ReadOnlySpan steps, + bool lastIsSeven, + Span candidates, + Span mask) + { + if (_disposed) + { + 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) + { + 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, + 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); + + 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; } + } +} 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..e293430 --- /dev/null +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorRemainderGpuStepper.cs @@ -0,0 +1,140 @@ +using System; +using System.Buffers; +using System.Runtime.InteropServices; +using ILGPU; +using ILGPU.Runtime; + +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> _kernel; + private MemoryBuffer1D _remainderBuffer; + private MemoryBuffer1D _stepBuffer; + private MemoryBuffer1D _resultBuffer; + private byte[] _hostRemainders; + private byte[] _hostSteps; + private readonly int _length; + private int _capacity; + private bool _disposed; + + internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan steps) + { + if (steps.Length != RemainderCount) + { + throw new ArgumentException("Step table must contain exactly six entries."); + } + + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = _accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(AdvanceRemaindersKernel); + + _length = steps.Length; + _capacity = Math.Max(1, _length); + _remainderBuffer = _accelerator.Allocate1D(_capacity); + _stepBuffer = _accelerator.Allocate1D(_capacity); + _resultBuffer = _accelerator.Allocate1D(_capacity); + + _hostRemainders = ArrayPool.Shared.Rent(_capacity); + _hostSteps = ArrayPool.Shared.Rent(_capacity); + + steps.CopyTo(_hostSteps.AsSpan(0, _length)); + + ArrayView1D stepView = _stepBuffer.View.SubView(0, _length); + stepView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostSteps.AsSpan(0, _length)), _length); + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _remainderBuffer.Dispose(); + _stepBuffer.Dispose(); + _resultBuffer.Dispose(); + ArrayPool.Shared.Return(_hostRemainders, clearArray: false); + ArrayPool.Shared.Return(_hostSteps, 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 resultView = _resultBuffer.View.SubView(0, _length); + + _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 results) + { + int idx = index; + if (idx >= remainders.Length) + { + return; + } + + int modulus = GetModulus(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; + } + + private static int GetModulus(int index) + { + return index switch + { + 0 => 10, + 1 => 8, + 2 => 5, + 3 => 3, + 4 => 7, + 5 => 11, + _ => 0, + }; + } +} 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); + } +}