Add GPU kernel for divisor candidate generation

PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs

@@ -17,6 +17,26 @@ internal static class DivisorByDivisorKernels
         hits[index] = montgomeryResult == 1UL ? (byte)1 : (byte)0;
     }
 
+    public static void GenerateCandidatesAndGapsKernel(Index1D index, ArrayView<ulong> candidates, ArrayView<ulong> gaps, ulong startValue, ulong stride)
+    {
+        int globalIndex = index;
+        int candidateLength = (int)candidates.Length;
+        if (globalIndex >= candidateLength)
+        {
+            return;
+        }
+
+        ulong offset = (ulong)globalIndex;
+        ulong candidate = startValue + offset * stride;
+        candidates[globalIndex] = candidate;
+
+        int gapLength = (int)gaps.Length;
+        if (globalIndex < gapLength)
+        {
+            gaps[globalIndex] = offset == 0UL ? 0UL : stride;
+        }
+    }
+
     public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView<ulong> gaps, ArrayView<byte> deltas, byte modulus)
     {
         deltas[index] = (byte)(gaps[index] % modulus);

PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs

@@ -28,6 +28,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
     private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<ulong>, ArrayView<byte>, byte>> _remainderDeltaKernelCache = new(AcceleratorReferenceComparer.Instance);
     private readonly ConcurrentDictionary<Accelerator, Action<KernelConfig, Index1D, ArrayView<byte>, ArrayView<byte>, byte, byte>> _remainderScanKernelCache = new(AcceleratorReferenceComparer.Instance);
     private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, byte, ArrayView<byte>>> _candidateMaskKernelCache = new(AcceleratorReferenceComparer.Instance);
+    private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<ulong>, ArrayView<ulong>, ulong, ulong>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance);
     private readonly ConcurrentDictionary<Accelerator, ConcurrentBag<BatchResources>> _resourcePools = new(AcceleratorReferenceComparer.Instance);
     private readonly ConcurrentBag<GpuContextPool.GpuContextLease> _acceleratorPool = new();
     private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
@@ -51,6 +52,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             accelerator,
             acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, byte, ArrayView<byte>>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel));
 
+    private Action<Index1D, ArrayView<ulong>, ArrayView<ulong>, ulong, ulong> GetCandidateGenerationKernel(Accelerator accelerator) =>
+        _candidateGenerationKernelCache.GetOrAdd(
+            accelerator,
+            acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<ulong>, ulong, ulong>(DivisorByDivisorKernels.GenerateCandidatesAndGapsKernel));
+
     public int GpuBatchSize
     {
         get => _gpuBatchSize;
@@ -125,12 +131,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             resources.DivisorsBuffer,
             resources.ExponentBuffer,
             resources.HitsBuffer,
+            resources.CandidateBuffer,
             resources.DivisorDeltaBuffer,
             resources.RemainderDeltaBuffer,
             resources.RemainderBuffer,
             resources.RemainderDeltaKernel,
             resources.RemainderScanKernel,
             resources.CandidateMaskKernel,
+            resources.CandidateGenerationKernel,
             resources.Divisors,
             resources.Exponents,
             resources.Hits,
@@ -193,12 +201,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         MemoryBuffer1D<MontgomeryDivisorData, Stride1D.Dense> divisorsBuffer,
         MemoryBuffer1D<ulong, Stride1D.Dense> exponentBuffer,
         MemoryBuffer1D<byte, Stride1D.Dense> hitsBuffer,
+        MemoryBuffer1D<ulong, Stride1D.Dense> candidateBuffer,
         MemoryBuffer1D<ulong, Stride1D.Dense> divisorDeltaBuffer,
         MemoryBuffer1D<byte, Stride1D.Dense> remainderDeltaBuffer,
         MemoryBuffer1D<byte, Stride1D.Dense> remainderBuffer,
         Action<Index1D, ArrayView<ulong>, ArrayView<byte>, byte> remainderDeltaKernel,
         Action<KernelConfig, Index1D, ArrayView<byte>, ArrayView<byte>, byte, byte> remainderScanKernel,
         Action<Index1D, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, byte, ArrayView<byte>> candidateMaskKernel,
+        Action<Index1D, ArrayView<ulong>, ArrayView<ulong>, ulong, ulong> candidateGenerationKernel,
         in ulong[] divisors,
         ulong[] exponents,
         byte[] hits,
@@ -258,10 +268,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         ArrayPool<ulong> ulongPool = ThreadStaticPools.UlongPool;
         ulong[] filteredDivisors = ulongPool.Rent(chunkCapacity);
-        ulong[] divisorGaps = ulongPool.Rent(chunkCapacity);
+        ulong[]? fallbackDivisorGaps = null;
         var viewCache = new BufferSliceCache(divisorDeltaBuffer, remainderDeltaBuffer, remainderBuffer, hitsBuffer, remainderStride);
         Span<ulong> candidateStorage = cycleCandidates.AsSpan();
-        Span<ulong> gapStorage = divisorGaps.AsSpan();
         Span<byte> hitStorage = hits.AsSpan();
         Span<ulong> filteredStorage = filteredDivisors.AsSpan();
         Span<ulong> cycleLengthStorage = cycleLengths.AsSpan();
@@ -283,28 +292,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             }
 
             Span<ulong> candidateSpan = candidateStorage[..effectiveCount];
-            Span<ulong> gapSpan = gapStorage[..effectiveCount];
             Span<byte> maskSpan = hitStorage[..effectiveCount];
 
-            UInt128 localDivisor = currentDivisor128;
-
-            for (int i = 0; i < effectiveCount; i++)
-            {
-                ulong divisorValue = (ulong)localDivisor;
-                candidateSpan[i] = divisorValue;
-                processedCount++;
-                lastProcessed = divisorValue;
-                localDivisor += twoP128;
-            }
-
-            gapSpan[0] = 0UL;
-            for (int i = 1; i < effectiveCount; i++)
-            {
-                gapSpan[i] = candidateSpan[i] - candidateSpan[i - 1];
-            }
-
-            currentDivisor128 = localDivisor;
-
             viewCache.Ensure(effectiveCount);
 
             ArrayView1D<ulong, Stride1D.Dense> gapView = viewCache.DivisorDelta;
@@ -315,7 +304,68 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             ArrayView1D<byte, Stride1D.Dense> remainder3View = viewCache.Remainder3;
             ArrayView1D<byte, Stride1D.Dense> maskView = viewCache.Hits;
 
-            gapView.CopyFromCPU(ref MemoryMarshal.GetReference(gapSpan), effectiveCount);
+            UInt128 localDivisor = currentDivisor128;
+            UInt128 nextDivisor;
+            bool strideFits = twoP128 <= (UInt128)ulong.MaxValue;
+            bool baseFits = localDivisor <= (UInt128)ulong.MaxValue;
+            bool useGpuGeneration = strideFits && baseFits;
+
+            if (useGpuGeneration)
+            {
+                ulong startValue = (ulong)localDivisor;
+                ulong stride = (ulong)twoP128;
+                ArrayView1D<ulong, Stride1D.Dense> candidateView = candidateBuffer.View.SubView(0, effectiveCount);
+                GenerateCandidatesOnGpu(effectiveCount, startValue, stride, candidateView, gapView, candidateGenerationKernel);
+                if (effectiveCount > 0)
+                {
+                    candidateView.CopyToCPU(ref MemoryMarshal.GetReference(candidateSpan), effectiveCount);
+                }
+
+                nextDivisor = localDivisor + twoP128 * (UInt128)(ulong)effectiveCount;
+            }
+            else
+            {
+                if (fallbackDivisorGaps is null || fallbackDivisorGaps.Length < chunkCapacity)
+                {
+                    if (fallbackDivisorGaps is not null)
+                    {
+                        ulongPool.Return(fallbackDivisorGaps, clearArray: false);
+                    }
+
+                    fallbackDivisorGaps = ulongPool.Rent(chunkCapacity);
+                }
+
+                Span<ulong> gapSpan = fallbackDivisorGaps.AsSpan(0, effectiveCount);
+                UInt128 cursor = localDivisor;
+
+                for (int i = 0; i < effectiveCount; i++)
+                {
+                    ulong divisorValue = (ulong)cursor;
+                    candidateSpan[i] = divisorValue;
+                    cursor += twoP128;
+                }
+
+                if (effectiveCount > 0)
+                {
+                    gapSpan[0] = 0UL;
+                    for (int i = 1; i < effectiveCount; i++)
+                    {
+                        gapSpan[i] = candidateSpan[i] - candidateSpan[i - 1];
+                    }
+
+                    gapView.CopyFromCPU(ref MemoryMarshal.GetReference(gapSpan), effectiveCount);
+                }
+
+                nextDivisor = cursor;
+            }
+
+            processedCount += (ulong)effectiveCount;
+            if (effectiveCount > 0)
+            {
+                lastProcessed = candidateSpan[effectiveCount - 1];
+            }
+
+            currentDivisor128 = nextDivisor;
 
             ComputeRemaindersOnGpu(effectiveCount, remainder10, 10, gapView, deltaView, remainder10View, remainderDeltaKernel, remainderScanKernel);
             ComputeRemaindersOnGpu(effectiveCount, remainder8, 8, gapView, deltaView, remainder8View, remainderDeltaKernel, remainderScanKernel);
@@ -397,7 +447,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         }
 
         ulongPool.Return(filteredDivisors, clearArray: false);
-        ulongPool.Return(divisorGaps, clearArray: false);
+        if (fallbackDivisorGaps is not null)
+        {
+            ulongPool.Return(fallbackDivisorGaps, clearArray: false);
+        }
 
         coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128);
         return composite;
@@ -474,6 +527,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
         internal ArrayView1D<byte, Stride1D.Dense> Hits => _hits;
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void GenerateCandidatesOnGpu(
+        int count,
+        ulong startValue,
+        ulong stride,
+        ArrayView1D<ulong, Stride1D.Dense> candidateView,
+        ArrayView1D<ulong, Stride1D.Dense> gapView,
+        Action<Index1D, ArrayView<ulong>, ArrayView<ulong>, ulong, ulong> kernel)
+    {
+        if (count <= 0)
+        {
+            return;
+        }
+
+        kernel(count, candidateView, gapView, startValue, stride);
+    }
+
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static void ComputeRemaindersOnGpu(
         int count,
@@ -735,6 +805,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             DivisorsBuffer = accelerator.Allocate1D<MontgomeryDivisorData>(actualCapacity);
             ExponentBuffer = accelerator.Allocate1D<ulong>(actualCapacity);
             HitsBuffer = accelerator.Allocate1D<byte>(actualCapacity);
+            CandidateBuffer = accelerator.Allocate1D<ulong>(actualCapacity);
             DivisorDeltaBuffer = accelerator.Allocate1D<ulong>(actualCapacity);
             RemainderDeltaBuffer = accelerator.Allocate1D<byte>(actualCapacity);
             RemainderBuffer = accelerator.Allocate1D<byte>(actualCapacity * RemainderSlotCount);
@@ -748,6 +819,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             RemainderDeltaKernel = owner.GetRemainderDeltaKernel(accelerator);
             RemainderScanKernel = owner.GetRemainderScanKernel(accelerator);
             CandidateMaskKernel = owner.GetCandidateMaskKernel(accelerator);
+            CandidateGenerationKernel = owner.GetCandidateGenerationKernel(accelerator);
             Capacity = actualCapacity;
             CycleCapacity = cycleCapacity;
         }
@@ -758,6 +830,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         internal readonly MemoryBuffer1D<byte, Stride1D.Dense> HitsBuffer;
 
+        internal readonly MemoryBuffer1D<ulong, Stride1D.Dense> CandidateBuffer;
+
         internal readonly MemoryBuffer1D<ulong, Stride1D.Dense> DivisorDeltaBuffer;
 
         internal readonly MemoryBuffer1D<byte, Stride1D.Dense> RemainderDeltaBuffer;
@@ -780,6 +854,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
 
         internal readonly Action<KernelConfig, Index1D, ArrayView<byte>, ArrayView<byte>, byte, byte> RemainderScanKernel;
         internal readonly Action<Index1D, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, ArrayView<byte>, byte, ArrayView<byte>> CandidateMaskKernel;
+        internal readonly Action<Index1D, ArrayView<ulong>, ArrayView<ulong>, ulong, ulong> CandidateGenerationKernel;
 
         internal readonly int Capacity;
 
@@ -790,6 +865,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv
             DivisorsBuffer.Dispose();
             ExponentBuffer.Dispose();
             HitsBuffer.Dispose();
+            CandidateBuffer.Dispose();
             DivisorDeltaBuffer.Dispose();
             RemainderDeltaBuffer.Dispose();
             RemainderBuffer.Dispose();