Add GPU options for by-divisor Montgomery data

PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs

@@ -0,0 +1,7 @@
+namespace PerfectNumbers.Core;
+
+public enum ByDivisorMontgomeryDevice
+{
+    Cpu = 0,
+    Gpu = 1,
+}

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<byte>.Shared.Rent(RemainderTableSize);
+        byte[] stepBuffer = ArrayPool<byte>.Shared.Rent(RemainderTableSize);
+        Span<byte> remainderSpan = remainderBuffer.AsSpan(0, RemainderTableSize);
+        Span<byte> stepSpan = stepBuffer.AsSpan(0, RemainderTableSize);
+        Span<byte> 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<MontgomeryDivisorData>.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<MontgomeryDivisorData> 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<MontgomeryDivisorData>.Shared.Return(montgomeryBuffer, clearArray: false);
+            }
+
+            remainderStepper.Dispose();
+            ArrayPool<byte>.Shared.Return(remainderBuffer, clearArray: false);
+            ArrayPool<byte>.Shared.Return(stepBuffer, clearArray: false);
             ArrayPool<ulong>.Shared.Return(candidateBuffer, clearArray: false);
             ArrayPool<byte>.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)

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<Index1D, ArrayView<ulong>, ArrayView<MontgomeryDivisorData>> _kernel;
+    private MemoryBuffer1D<ulong, Stride1D.Dense> _modulusBuffer;
+    private MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> _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<Index1D, ArrayView<ulong>, ArrayView<MontgomeryDivisorData>>(BuildKernel);
+        _capacity = Math.Max(1, batchSize);
+        _modulusBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _resultBuffer = _accelerator.Allocate1D<MontgomeryDivisorData>(_capacity);
+        _hostModulusBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
+        _hostResultBuffer = ArrayPool<MontgomeryDivisorData>.Shared.Rent(_capacity);
+    }
+
+    public void Dispose()
+    {
+        if (_disposed)
+        {
+            return;
+        }
+
+        _disposed = true;
+        _modulusBuffer.Dispose();
+        _resultBuffer.Dispose();
+        ArrayPool<ulong>.Shared.Return(_hostModulusBuffer, clearArray: false);
+        ArrayPool<MontgomeryDivisorData>.Shared.Return(_hostResultBuffer, clearArray: false);
+        _lease.Dispose();
+    }
+
+    internal void Build(ReadOnlySpan<ulong> moduli, Span<MontgomeryDivisorData> destination)
+    {
+        if (_disposed)
+        {
+            throw new ObjectDisposedException(nameof(MersenneNumberDivisorMontgomeryGpuBuilder));
+        }
+
+        int length = moduli.Length;
+        if (length == 0)
+        {
+            return;
+        }
+
+        EnsureCapacity(length);
+
+        Span<ulong> hostModulusSpan = _hostModulusBuffer.AsSpan(0, length);
+        moduli.CopyTo(hostModulusSpan);
+
+        ArrayView1D<ulong, Stride1D.Dense> modulusView = _modulusBuffer.View.SubView(0, length);
+        modulusView.CopyFromCPU(ref MemoryMarshal.GetReference(hostModulusSpan), length);
+
+        ArrayView1D<MontgomeryDivisorData, Stride1D.Dense> resultView = _resultBuffer.View.SubView(0, length);
+        _kernel(length, modulusView, resultView);
+
+        Span<MontgomeryDivisorData> 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<ulong>.Shared.Return(_hostModulusBuffer, clearArray: false);
+        ArrayPool<MontgomeryDivisorData>.Shared.Return(_hostResultBuffer, clearArray: false);
+
+        _capacity = newCapacity;
+        _modulusBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _resultBuffer = _accelerator.Allocate1D<MontgomeryDivisorData>(_capacity);
+        _hostModulusBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
+        _hostResultBuffer = ArrayPool<MontgomeryDivisorData>.Shared.Rent(_capacity);
+    }
+
+    private static void BuildKernel(Index1D index, ArrayView<ulong> moduli, ArrayView<MontgomeryDivisorData> 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);
+    }
+}

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<Index1D, int, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>> _kernel;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _remainderBuffer;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _stepBuffer;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _modulusBuffer;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _resultBuffer;
+    private byte[] _hostRemainders;
+    private byte[] _hostSteps;
+    private byte[] _hostModuli;
+    private readonly int _length;
+    private int _capacity;
+    private bool _disposed;
+
+    internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan<byte> steps, ReadOnlySpan<byte> 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<Index1D, int, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>>(AdvanceRemaindersKernel);
+
+        _length = steps.Length;
+        _capacity = Math.Max(1, _length);
+        _remainderBuffer = _accelerator.Allocate1D<byte>(_capacity);
+        _stepBuffer = _accelerator.Allocate1D<byte>(_capacity);
+        _modulusBuffer = _accelerator.Allocate1D<byte>(_capacity);
+        _resultBuffer = _accelerator.Allocate1D<byte>(_capacity);
+
+        _hostRemainders = ArrayPool<byte>.Shared.Rent(_capacity);
+        _hostSteps = ArrayPool<byte>.Shared.Rent(_capacity);
+        _hostModuli = ArrayPool<byte>.Shared.Rent(_capacity);
+
+        steps.CopyTo(_hostSteps.AsSpan(0, _length));
+        moduli.CopyTo(_hostModuli.AsSpan(0, _length));
+
+        ArrayView1D<byte, Stride1D.Dense> stepView = _stepBuffer.View.SubView(0, _length);
+        stepView.CopyFromCPU(ref MemoryMarshal.GetReference(_hostSteps.AsSpan(0, _length)), _length);
+
+        ArrayView1D<byte, Stride1D.Dense> 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<byte>.Shared.Return(_hostRemainders, clearArray: false);
+        ArrayPool<byte>.Shared.Return(_hostSteps, clearArray: false);
+        ArrayPool<byte>.Shared.Return(_hostModuli, clearArray: false);
+        _lease.Dispose();
+    }
+
+    internal void Advance(int advanceCount, Span<byte> 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<byte> hostRemainderSpan = _hostRemainders.AsSpan(0, _length);
+        remainders.CopyTo(hostRemainderSpan);
+
+        ArrayView1D<byte, Stride1D.Dense> remainderView = _remainderBuffer.View.SubView(0, _length);
+        remainderView.CopyFromCPU(ref MemoryMarshal.GetReference(hostRemainderSpan), _length);
+
+        ArrayView1D<byte, Stride1D.Dense> stepView = _stepBuffer.View.SubView(0, _length);
+        ArrayView1D<byte, Stride1D.Dense> modulusView = _modulusBuffer.View.SubView(0, _length);
+        ArrayView1D<byte, Stride1D.Dense> 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<byte> remainders, ArrayView<byte> steps, ArrayView<byte> moduli, ArrayView<byte> 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;
+    }
+}