Merge pull request #162 from mlemanczyk/codex/analyze-execution-path-for-gpu-optimization

Batch CPU divisor scanning in reusable batches

PerfectNumbers.Core/ByDivisorDeltasDevice.cs

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

PerfectNumbers.Core/ByDivisorMontgomeryDevice.cs

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

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<ulong>.Shared.Rent(batchCapacity);
+        byte[] maskBuffer = ArrayPool<byte>.Shared.Rent(batchCapacity);
+
+        Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>? factorCache = null;
+
+        ulong cachedModulus = 0UL;
+        MontgomeryDivisorData cachedDivisorData = default;
+        bool hasCachedDivisorData = false;
+
+        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);
+        // 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<MontgomeryDivisorData>.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<ulong> candidateSpan = candidateBuffer.AsSpan(0, chunkCount);
+                Span<byte> maskSpan = maskBuffer.AsSpan(0, chunkCount);
+
+                candidateEvaluator.EvaluateCandidates(
+                    currentDivisor,
+                    step,
+                    limit64,
+                    remainderSpan,
+                    stepSpan,
+                    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())
+                    {
+                        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<ulong, MersenneDivisorCycles.FactorCacheEntry>(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<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 = 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<ulong> primes,
+            Span<byte> 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<ulong> primes,
+            Span<byte> hits,
+            int length)
+        {
+            EnsureCapacity(length);
+
+            Span<ulong> deltaSpan = _primeDeltas!.AsSpan(0, length);
+            ComputePrimeDeltas(primes, deltaSpan);
+
+            Span<ulong> remainderSpan = _cycleRemainders!.AsSpan(0, length);
+            ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan);
+
+            EnsureGpuEvaluator();
+            Span<ulong> 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<ulong>.Shared.Rent(newCapacity);
+            _cycleRemainders = ArrayPool<ulong>.Shared.Rent(newCapacity);
+            _montgomeryDeltas = ArrayPool<ulong>.Shared.Rent(newCapacity);
+            _bufferCapacity = newCapacity;
+        }
+
+        private void ReturnBuffers()
+        {
+            if (_primeDeltas is not null)
+            {
+                ArrayPool<ulong>.Shared.Return(_primeDeltas, clearArray: false);
+                _primeDeltas = null;
+            }
+
+            if (_cycleRemainders is not null)
+            {
+                ArrayPool<ulong>.Shared.Return(_cycleRemainders, clearArray: false);
+                _cycleRemainders = null;
+            }
+
+            if (_montgomeryDeltas is not null)
+            {
+                ArrayPool<ulong>.Shared.Return(_montgomeryDeltas, clearArray: false);
+                _montgomeryDeltas = null;
+            }
+
+            _bufferCapacity = 0;
+        }
+
+        private static void ComputePrimeDeltas(ReadOnlySpan<ulong> primes, Span<ulong> 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<ulong> primes, ReadOnlySpan<ulong> deltas, Span<ulong> 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);
         }
     }

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<Index1D, CandidateFilterArgs, ArrayView<ulong>, ArrayView<byte>> _kernel;
+    private MemoryBuffer1D<ulong, Stride1D.Dense> _candidateBuffer;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _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<Index1D, CandidateFilterArgs, ArrayView<ulong>, ArrayView<byte>>(EvaluateCandidatesKernel);
+        _capacity = Math.Max(1, batchSize);
+        _candidateBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _maskBuffer = _accelerator.Allocate1D<byte>(_capacity);
+        _hostCandidates = ArrayPool<ulong>.Shared.Rent(_capacity);
+        _hostMask = ArrayPool<byte>.Shared.Rent(_capacity);
+    }
+
+    public void Dispose()
+    {
+        if (_disposed)
+        {
+            return;
+        }
+
+        _disposed = true;
+        _candidateBuffer.Dispose();
+        _maskBuffer.Dispose();
+        ArrayPool<ulong>.Shared.Return(_hostCandidates, clearArray: false);
+        ArrayPool<byte>.Shared.Return(_hostMask, clearArray: false);
+        _lease.Dispose();
+    }
+
+    internal void EvaluateCandidates(
+        ulong startDivisor,
+        ulong step,
+        ulong limit,
+        ReadOnlySpan<byte> remainders,
+        ReadOnlySpan<byte> steps,
+        bool lastIsSeven,
+        Span<ulong> candidates,
+        Span<byte> 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<ulong> hostCandidateSpan = _hostCandidates.AsSpan(0, count);
+        Span<byte> hostMaskSpan = _hostMask.AsSpan(0, count);
+
+        ArrayView1D<ulong, Stride1D.Dense> candidateView = _candidateBuffer.View.SubView(0, count);
+        ArrayView1D<byte, Stride1D.Dense> 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<ulong>.Shared.Return(_hostCandidates, clearArray: false);
+        ArrayPool<byte>.Shared.Return(_hostMask, clearArray: false);
+
+        _capacity = newCapacity;
+        _candidateBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _maskBuffer = _accelerator.Allocate1D<byte>(_capacity);
+        _hostCandidates = ArrayPool<ulong>.Shared.Rent(_capacity);
+        _hostMask = ArrayPool<byte>.Shared.Rent(_capacity);
+    }
+
+    private static void EvaluateCandidatesKernel(
+        Index1D index,
+        CandidateFilterArgs args,
+        ArrayView<ulong> candidates,
+        ArrayView<byte> 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; }
+    }
+}

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,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<Index1D, int, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>> _kernel;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _remainderBuffer;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _stepBuffer;
+    private MemoryBuffer1D<byte, Stride1D.Dense> _resultBuffer;
+    private byte[] _hostRemainders;
+    private byte[] _hostSteps;
+    private readonly int _length;
+    private int _capacity;
+    private bool _disposed;
+
+    internal MersenneNumberDivisorRemainderGpuStepper(ReadOnlySpan<byte> 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<Index1D, int, 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);
+        _resultBuffer = _accelerator.Allocate1D<byte>(_capacity);
+
+        _hostRemainders = ArrayPool<byte>.Shared.Rent(_capacity);
+        _hostSteps = ArrayPool<byte>.Shared.Rent(_capacity);
+
+        steps.CopyTo(_hostSteps.AsSpan(0, _length));
+
+        ArrayView1D<byte, Stride1D.Dense> 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<byte>.Shared.Return(_hostRemainders, clearArray: false);
+        ArrayPool<byte>.Shared.Return(_hostSteps, 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> 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<byte> remainders, ArrayView<byte> steps, ArrayView<byte> 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,
+        };
+    }
+}

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<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel;
+    private MemoryBuffer1D<ulong, Stride1D.Dense> _exponentBuffer;
+    private MemoryBuffer1D<ulong, Stride1D.Dense> _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<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputeResiduesKernel);
+        _capacity = Math.Max(1, batchSize);
+        _exponentBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _resultBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
+    }
+
+    public void Dispose()
+    {
+        if (_disposed)
+        {
+            return;
+        }
+
+        _disposed = true;
+        _exponentBuffer.Dispose();
+        _resultBuffer.Dispose();
+        ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
+        _lease.Dispose();
+    }
+
+    internal void ComputeResidues(ReadOnlySpan<ulong> exponents, in MontgomeryDivisorData divisorData, Span<ulong> destination)
+    {
+        if (_disposed)
+        {
+            throw new ObjectDisposedException(nameof(MersenneNumberDivisorResidueGpuEvaluator));
+        }
+
+        int length = exponents.Length;
+        if (length == 0)
+        {
+            return;
+        }
+
+        EnsureCapacity(length);
+
+        Span<ulong> hostSpan = _hostBuffer.AsSpan(0, length);
+        exponents.CopyTo(hostSpan);
+
+        ArrayView1D<ulong, Stride1D.Dense> exponentView = _exponentBuffer.View.SubView(0, length);
+        exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), length);
+
+        ArrayView1D<ulong, Stride1D.Dense> 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<ulong>.Shared.Return(_hostBuffer, clearArray: false);
+
+        _capacity = newCapacity;
+        _exponentBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _resultBuffer = _accelerator.Allocate1D<ulong>(_capacity);
+        _hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
+    }
+
+    private static void ComputeResiduesKernel(Index1D index, MontgomeryDivisorData divisorData, ArrayView<ulong> exponents, ArrayView<ulong> results)
+    {
+        int idx = index;
+        int length = (int)exponents.Length;
+        if (idx >= length)
+        {
+            return;
+        }
+
+        ulong exponent = exponents[idx];
+        results[idx] = exponent.Pow2MontgomeryModWindowedGpu(divisorData, keepMontgomery: true);
+    }
+}