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; + } +}