From bfc5712d81cb37642ba51e112d2c18fe75957632 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 05:24:01 +0200 Subject: [PATCH] Add GPU kernel for divisor candidate generation --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 20 +++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 124 ++++++++++++++---- 2 files changed, 120 insertions(+), 24 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 068e718..3f026c0 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/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 candidates, ArrayView 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 gaps, ArrayView deltas, byte modulus) { deltas[index] = (byte)(gaps[index] % modulus); diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 4fbe697..a0d6219 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -28,6 +28,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly ConcurrentDictionary, ArrayView, byte>> _remainderDeltaKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView, byte, byte>> _remainderScanKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView, byte, ArrayView>> _candidateMaskKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView, ulong, ulong>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentBag _sessionPool = new(); @@ -51,6 +52,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, byte, ArrayView>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); + private Action, ArrayView, ulong, ulong> GetCandidateGenerationKernel(Accelerator accelerator) => + _candidateGenerationKernelCache.GetOrAdd( + accelerator, + acc => acc.LoadAutoGroupedStreamKernel, ArrayView, 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 divisorsBuffer, MemoryBuffer1D exponentBuffer, MemoryBuffer1D hitsBuffer, + MemoryBuffer1D candidateBuffer, MemoryBuffer1D divisorDeltaBuffer, MemoryBuffer1D remainderDeltaBuffer, MemoryBuffer1D remainderBuffer, Action, ArrayView, byte> remainderDeltaKernel, Action, ArrayView, byte, byte> remainderScanKernel, Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> candidateMaskKernel, + Action, ArrayView, ulong, ulong> candidateGenerationKernel, in ulong[] divisors, ulong[] exponents, byte[] hits, @@ -258,10 +268,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayPool 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 candidateStorage = cycleCandidates.AsSpan(); - Span gapStorage = divisorGaps.AsSpan(); Span hitStorage = hits.AsSpan(); Span filteredStorage = filteredDivisors.AsSpan(); Span cycleLengthStorage = cycleLengths.AsSpan(); @@ -283,28 +292,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } Span candidateSpan = candidateStorage[..effectiveCount]; - Span gapSpan = gapStorage[..effectiveCount]; Span 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 gapView = viewCache.DivisorDelta; @@ -315,7 +304,68 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D remainder3View = viewCache.Remainder3; ArrayView1D 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 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 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 Hits => _hits; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void GenerateCandidatesOnGpu( + int count, + ulong startValue, + ulong stride, + ArrayView1D candidateView, + ArrayView1D gapView, + Action, ArrayView, 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(actualCapacity); ExponentBuffer = accelerator.Allocate1D(actualCapacity); HitsBuffer = accelerator.Allocate1D(actualCapacity); + CandidateBuffer = accelerator.Allocate1D(actualCapacity); DivisorDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderBuffer = accelerator.Allocate1D(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 HitsBuffer; + internal readonly MemoryBuffer1D CandidateBuffer; + internal readonly MemoryBuffer1D DivisorDeltaBuffer; internal readonly MemoryBuffer1D RemainderDeltaBuffer; @@ -780,6 +854,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly Action, ArrayView, byte, byte> RemainderScanKernel; internal readonly Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> CandidateMaskKernel; + internal readonly Action, ArrayView, 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();