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 01/10] 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(); From 4d8ea9dccd594ff075c9c1ee3c0e58a1bec2103a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 06:05:22 +0200 Subject: [PATCH 02/10] Use GpuUInt128 in GPU divisor candidate generation --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 45 +++++++++-- ...MersenneNumberDivisorByDivisorGpuTester.cs | 76 ++++--------------- 2 files changed, 55 insertions(+), 66 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 3f026c0..abeb8da 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -17,7 +17,12 @@ 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) + public static void GenerateCandidatesAndGapsKernel( + Index1D index, + ArrayView candidates, + ArrayView gaps, + GpuUInt128 startValue, + GpuUInt128 stride) { int globalIndex = index; int candidateLength = (int)candidates.Length; @@ -27,14 +32,44 @@ internal static class DivisorByDivisorKernels } ulong offset = (ulong)globalIndex; - ulong candidate = startValue + offset * stride; - candidates[globalIndex] = candidate; + GpuUInt128 candidateValue; + if (offset == 0UL) + { + candidateValue = startValue; + } + else + { + candidateValue = stride; + candidateValue.Mul(offset); + candidateValue.Add(startValue); + } + + candidates[globalIndex] = candidateValue.Low; int gapLength = (int)gaps.Length; - if (globalIndex < gapLength) + if (globalIndex >= gapLength) { - gaps[globalIndex] = offset == 0UL ? 0UL : stride; + return; } + + if (offset == 0UL) + { + gaps[globalIndex] = 0UL; + return; + } + + GpuUInt128 previousValue = startValue; + ulong previousOffset = offset - 1UL; + if (previousOffset != 0UL) + { + GpuUInt128 scaledStride = stride; + scaledStride.Mul(previousOffset); + previousValue.Add(scaledStride); + } + + GpuUInt128 gapValue = candidateValue; + gapValue.Sub(previousValue); + gaps[globalIndex] = gapValue.Low; } public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index a0d6219..5a988bf 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -28,7 +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, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentBag _sessionPool = new(); @@ -52,10 +52,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, byte, ArrayView>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); - private Action, ArrayView, ulong, ulong> GetCandidateGenerationKernel(Accelerator accelerator) => + private Action, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => _candidateGenerationKernelCache.GetOrAdd( accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ulong, ulong>(DivisorByDivisorKernels.GenerateCandidatesAndGapsKernel)); + acc => acc.LoadAutoGroupedStreamKernel, ArrayView, GpuUInt128, GpuUInt128>(DivisorByDivisorKernels.GenerateCandidatesAndGapsKernel)); public int GpuBatchSize { @@ -208,7 +208,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Action, ArrayView, byte> remainderDeltaKernel, Action, ArrayView, byte, byte> remainderScanKernel, Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> candidateMaskKernel, - Action, ArrayView, ulong, ulong> candidateGenerationKernel, + Action, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, in ulong[] divisors, ulong[] exponents, byte[] hits, @@ -268,7 +268,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayPool ulongPool = ThreadStaticPools.UlongPool; ulong[] filteredDivisors = ulongPool.Rent(chunkCapacity); - ulong[]? fallbackDivisorGaps = null; var viewCache = new BufferSliceCache(divisorDeltaBuffer, remainderDeltaBuffer, remainderBuffer, hitsBuffer, remainderStride); Span candidateStorage = cycleCandidates.AsSpan(); Span hitStorage = hits.AsSpan(); @@ -306,58 +305,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv UInt128 localDivisor = currentDivisor128; UInt128 nextDivisor; - bool strideFits = twoP128 <= (UInt128)ulong.MaxValue; - bool baseFits = localDivisor <= (UInt128)ulong.MaxValue; - bool useGpuGeneration = strideFits && baseFits; - if (useGpuGeneration) + ArrayView1D candidateView = candidateBuffer.View.SubView(0, effectiveCount); + GpuUInt128 gpuStartValue = new GpuUInt128(localDivisor); + GpuUInt128 gpuStride = new GpuUInt128(twoP128); + GenerateCandidatesOnGpu(effectiveCount, gpuStartValue, gpuStride, candidateView, gapView, candidateGenerationKernel); + if (effectiveCount > 0) { - 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; + candidateView.CopyToCPU(ref MemoryMarshal.GetReference(candidateSpan), 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; - } + nextDivisor = localDivisor + twoP128 * (UInt128)effectiveCount; processedCount += (ulong)effectiveCount; if (effectiveCount > 0) @@ -447,10 +405,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } ulongPool.Return(filteredDivisors, clearArray: false); - if (fallbackDivisorGaps is not null) - { - ulongPool.Return(fallbackDivisorGaps, clearArray: false); - } coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); return composite; @@ -530,11 +484,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void GenerateCandidatesOnGpu( int count, - ulong startValue, - ulong stride, + GpuUInt128 startValue, + GpuUInt128 stride, ArrayView1D candidateView, ArrayView1D gapView, - Action, ArrayView, ulong, ulong> kernel) + Action, ArrayView, GpuUInt128, GpuUInt128> kernel) { if (count <= 0) { @@ -854,7 +808,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 Action, ArrayView, GpuUInt128, GpuUInt128> CandidateGenerationKernel; internal readonly int Capacity; From 5b973fb513d03c6628c536e3b65b9243498a2a5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 06:43:00 +0200 Subject: [PATCH 03/10] Document and streamline GPU candidate kernel guards --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 57 +++++++------------ 1 file changed, 21 insertions(+), 36 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index abeb8da..35e4170 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -25,47 +25,32 @@ internal static class DivisorByDivisorKernels GpuUInt128 stride) { int globalIndex = index; - int candidateLength = (int)candidates.Length; - if (globalIndex >= candidateLength) - { - return; - } + // EvenPerfectBitScanner always launches this kernel with an extent matching the candidate view length, + // so the bounds guard stays commented out to keep the hot path branch-free. + // int candidateLength = (int)candidates.Length; + // if (globalIndex >= candidateLength) + // { + // return; + // } ulong offset = (ulong)globalIndex; - GpuUInt128 candidateValue; - if (offset == 0UL) - { - candidateValue = startValue; - } - else - { - candidateValue = stride; - candidateValue.Mul(offset); - candidateValue.Add(startValue); - } - + GpuUInt128 candidateValue = stride; + candidateValue.Mul(offset); + candidateValue.Add(startValue); candidates[globalIndex] = candidateValue.Low; - int gapLength = (int)gaps.Length; - if (globalIndex >= gapLength) - { - return; - } + // The scanner supplies gap buffers that mirror the candidate view length, + // so the additional bounds guard remains commented out for the same reason as above. + // int gapLength = (int)gaps.Length; + // if (globalIndex >= gapLength) + // { + // return; + // } - if (offset == 0UL) - { - gaps[globalIndex] = 0UL; - return; - } - - GpuUInt128 previousValue = startValue; - ulong previousOffset = offset - 1UL; - if (previousOffset != 0UL) - { - GpuUInt128 scaledStride = stride; - scaledStride.Mul(previousOffset); - previousValue.Add(scaledStride); - } + ulong previousOffset = offset > 0UL ? offset - 1UL : 0UL; + GpuUInt128 previousValue = stride; + previousValue.Mul(previousOffset); + previousValue.Add(startValue); GpuUInt128 gapValue = candidateValue; gapValue.Sub(previousValue); From bae39ad95e338f25118025daa3757d76d5498a3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 06:43:05 +0200 Subject: [PATCH 04/10] Capture GPU next start for divisor batches --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 17 +++++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 69 +++++++++++-------- 2 files changed, 59 insertions(+), 27 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 35e4170..04f67f7 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -21,6 +21,7 @@ internal static class DivisorByDivisorKernels Index1D index, ArrayView candidates, ArrayView gaps, + ArrayView nextStarts, GpuUInt128 startValue, GpuUInt128 stride) { @@ -55,6 +56,22 @@ internal static class DivisorByDivisorKernels GpuUInt128 gapValue = candidateValue; gapValue.Sub(previousValue); gaps[globalIndex] = gapValue.Low; + + if (globalIndex == 0) + { + // EvenPerfectBitScanner provides a single-element nextStarts view per launch, so the guard stays commented out. + // int nextStartLength = (int)nextStarts.Length; + // if (nextStartLength <= 0) + // { + // return; + // } + + ulong candidateCount = (ulong)candidates.Length; + GpuUInt128 nextValue = stride; + nextValue.Mul(candidateCount); + nextValue.Add(startValue); + nextStarts[0] = nextValue; + } } public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 5a988bf..70e1292 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -28,7 +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, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); private readonly ConcurrentBag _sessionPool = new(); @@ -52,10 +52,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, byte, ArrayView>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); - private Action, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => + private Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => _candidateGenerationKernelCache.GetOrAdd( accelerator, - acc => acc.LoadAutoGroupedStreamKernel, ArrayView, GpuUInt128, GpuUInt128>(DivisorByDivisorKernels.GenerateCandidatesAndGapsKernel)); + acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, GpuUInt128, GpuUInt128>(DivisorByDivisorKernels.GenerateCandidatesAndGapsKernel)); public int GpuBatchSize { @@ -132,6 +132,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.ExponentBuffer, resources.HitsBuffer, resources.CandidateBuffer, + resources.NextStartBuffer, resources.DivisorDeltaBuffer, resources.RemainderDeltaBuffer, resources.RemainderBuffer, @@ -202,13 +203,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D exponentBuffer, MemoryBuffer1D hitsBuffer, MemoryBuffer1D candidateBuffer, + MemoryBuffer1D nextStartBuffer, MemoryBuffer1D divisorDeltaBuffer, MemoryBuffer1D remainderDeltaBuffer, MemoryBuffer1D remainderBuffer, Action, ArrayView, byte> remainderDeltaKernel, Action, ArrayView, byte, byte> remainderScanKernel, Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> candidateMaskKernel, - Action, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, + Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, in ulong[] divisors, ulong[] exponents, byte[] hits, @@ -304,27 +306,27 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D maskView = viewCache.Hits; UInt128 localDivisor = currentDivisor128; - UInt128 nextDivisor; ArrayView1D candidateView = candidateBuffer.View.SubView(0, effectiveCount); + ArrayView1D nextStartView = nextStartBuffer.View.SubView(0, 1); GpuUInt128 gpuStartValue = new GpuUInt128(localDivisor); GpuUInt128 gpuStride = new GpuUInt128(twoP128); - GenerateCandidatesOnGpu(effectiveCount, gpuStartValue, gpuStride, candidateView, gapView, candidateGenerationKernel); + GpuUInt128 nextDivisorGpu = GenerateCandidatesOnGpu( + effectiveCount, + gpuStartValue, + gpuStride, + candidateView, + gapView, + nextStartView, + candidateGenerationKernel); if (effectiveCount > 0) { candidateView.CopyToCPU(ref MemoryMarshal.GetReference(candidateSpan), effectiveCount); - } - - nextDivisor = localDivisor + twoP128 * (UInt128)effectiveCount; - - processedCount += (ulong)effectiveCount; - if (effectiveCount > 0) - { + processedCount += (ulong)effectiveCount; lastProcessed = candidateSpan[effectiveCount - 1]; + currentDivisor128 = (UInt128)nextDivisorGpu; } - currentDivisor128 = nextDivisor; - ComputeRemaindersOnGpu(effectiveCount, remainder10, 10, gapView, deltaView, remainder10View, remainderDeltaKernel, remainderScanKernel); ComputeRemaindersOnGpu(effectiveCount, remainder8, 8, gapView, deltaView, remainder8View, remainderDeltaKernel, remainderScanKernel); ComputeRemaindersOnGpu(effectiveCount, remainder5, 5, gapView, deltaView, remainder5View, remainderDeltaKernel, remainderScanKernel); @@ -482,20 +484,27 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void GenerateCandidatesOnGpu( + private static GpuUInt128 GenerateCandidatesOnGpu( int count, GpuUInt128 startValue, GpuUInt128 stride, ArrayView1D candidateView, ArrayView1D gapView, - Action, ArrayView, GpuUInt128, GpuUInt128> kernel) + ArrayView1D nextStartView, + Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> kernel) { - if (count <= 0) - { - return; - } + // EvenPerfectBitScanner always launches this helper with a positive count, + // so the guard stays commented out to keep the hot path branch-free. + // if (count <= 0) + // { + // return GpuUInt128.Zero; + // } - kernel(count, candidateView, gapView, startValue, stride); + kernel(count, candidateView, gapView, nextStartView, startValue, stride); + + GpuUInt128 nextStart = default; + nextStartView.SubView(0, 1).CopyToCPU(ref nextStart, 1); + return nextStart; } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -509,10 +518,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Action, ArrayView, byte> deltaKernel, Action, ArrayView, byte, byte> scanKernel) { - if (count <= 0) - { - return; - } + // EvenPerfectBitScanner always supplies a positive count here, + // so the guard stays commented out to keep the hot path branch-free. + // if (count <= 0) + // { + // return; + // } deltaKernel(count, gapView, deltaView, modulus); SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count); @@ -760,6 +771,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ExponentBuffer = accelerator.Allocate1D(actualCapacity); HitsBuffer = accelerator.Allocate1D(actualCapacity); CandidateBuffer = accelerator.Allocate1D(actualCapacity); + NextStartBuffer = accelerator.Allocate1D(1); DivisorDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderBuffer = accelerator.Allocate1D(actualCapacity * RemainderSlotCount); @@ -786,6 +798,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly MemoryBuffer1D CandidateBuffer; + internal readonly MemoryBuffer1D NextStartBuffer; + internal readonly MemoryBuffer1D DivisorDeltaBuffer; internal readonly MemoryBuffer1D RemainderDeltaBuffer; @@ -808,7 +822,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, GpuUInt128, GpuUInt128> CandidateGenerationKernel; + internal readonly Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> CandidateGenerationKernel; internal readonly int Capacity; @@ -820,6 +834,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ExponentBuffer.Dispose(); HitsBuffer.Dispose(); CandidateBuffer.Dispose(); + NextStartBuffer.Dispose(); DivisorDeltaBuffer.Dispose(); RemainderDeltaBuffer.Dispose(); RemainderBuffer.Dispose(); From add62df697017f0eaf5fd8661f28dff5c1fb4483 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 18:47:55 +0200 Subject: [PATCH 05/10] Compact GPU divisor masks on device --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 57 +++++++++++++++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 69 +++++++++++++++---- 2 files changed, 114 insertions(+), 12 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 04f67f7..7e6258f 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -155,6 +155,63 @@ internal static class DivisorByDivisorKernels mask[globalIndex] = accepted; } + public static void CompactCandidateMaskKernel( + Index1D index, + ArrayView mask, + ArrayView sourceCandidates, + ArrayView compactedCandidates, + ArrayView compactedCount) + { + int globalIndex = index; + int length = (int)mask.Length; + if (globalIndex >= length) + { + return; + } + + int localIndex = Group.IdxX; + int groupSize = Group.Dimension.X; + var shared = SharedMemory.GetDynamic(); + + byte maskValue = mask[globalIndex]; + int accepted = maskValue != 0 ? 1 : 0; + shared[localIndex] = accepted; + + ulong candidateValue = sourceCandidates[globalIndex]; + + Group.Barrier(); + + int offset = 1; + while (offset < groupSize) + { + int addend = localIndex >= offset ? shared[localIndex - offset] : 0; + + Group.Barrier(); + + int current = shared[localIndex]; + int sum = current + addend; + if (localIndex >= offset) + { + shared[localIndex] = sum; + } + + Group.Barrier(); + + offset <<= 1; + } + + if (accepted != 0) + { + int targetIndex = shared[localIndex] - 1; + compactedCandidates[targetIndex] = candidateValue; + } + + if (globalIndex == length - 1) + { + compactedCount[0] = shared[localIndex]; + } + } + public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 70e1292..a5b7955 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, ArrayView, ArrayView>> _candidateCompactionKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); @@ -52,6 +53,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, byte, ArrayView>(DivisorByDivisorKernels.EvaluateCandidateMaskKernel)); + private Action, ArrayView, ArrayView, ArrayView> GetCandidateCompactionKernel(Accelerator accelerator) => + _candidateCompactionKernelCache.GetOrAdd( + accelerator, + acc => acc.LoadStreamKernel, ArrayView, ArrayView, ArrayView>(DivisorByDivisorKernels.CompactCandidateMaskKernel)); + private Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => _candidateGenerationKernelCache.GetOrAdd( accelerator, @@ -136,9 +142,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.DivisorDeltaBuffer, resources.RemainderDeltaBuffer, resources.RemainderBuffer, + resources.FilterCountBuffer, resources.RemainderDeltaKernel, resources.RemainderScanKernel, resources.CandidateMaskKernel, + resources.CandidateCompactionKernel, resources.CandidateGenerationKernel, resources.Divisors, resources.Exponents, @@ -207,9 +215,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D divisorDeltaBuffer, MemoryBuffer1D remainderDeltaBuffer, MemoryBuffer1D remainderBuffer, + MemoryBuffer1D filterCountBuffer, Action, ArrayView, byte> remainderDeltaKernel, Action, ArrayView, byte, byte> remainderScanKernel, Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> candidateMaskKernel, + Action, ArrayView, ArrayView, ArrayView> candidateCompactionKernel, Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, in ulong[] divisors, ulong[] exponents, @@ -293,7 +303,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } Span candidateSpan = candidateStorage[..effectiveCount]; - Span maskSpan = hitStorage[..effectiveCount]; viewCache.Ensure(effectiveCount); @@ -304,6 +313,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D remainder5View = viewCache.Remainder5; ArrayView1D remainder3View = viewCache.Remainder3; ArrayView1D maskView = viewCache.Hits; + ArrayView1D filteredCountView = filterCountBuffer.View; UInt128 localDivisor = currentDivisor128; @@ -332,7 +342,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ComputeRemaindersOnGpu(effectiveCount, remainder5, 5, gapView, deltaView, remainder5View, remainderDeltaKernel, remainderScanKernel); ComputeRemaindersOnGpu(effectiveCount, remainder3, 3, gapView, deltaView, remainder3View, remainderDeltaKernel, remainderScanKernel); - BuildCandidateMaskOnGpu(effectiveCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, maskSpan, candidateMaskKernel); + BuildCandidateMaskOnGpu(effectiveCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, candidateMaskKernel); byte lastRemainder10 = FetchLastRemainder(remainder10View, effectiveCount); byte lastRemainder8 = FetchLastRemainder(remainder8View, effectiveCount); @@ -344,14 +354,21 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv remainder5 = AddModGpu(lastRemainder5, step5, 5); remainder3 = AddModGpu(lastRemainder3, step3, 3); - int filteredCount = 0; - Span filteredDivisorsSpan = filteredStorage[..effectiveCount]; + int filteredCount = CompactCandidatesOnGpu( + effectiveCount, + maskView, + candidateView, + candidateView, + filteredCountView, + candidateCompactionKernel); - for (int i = 0; i < effectiveCount; i++) + if (filteredCount > 0) { - int add = maskSpan[i] == 0 ? 0 : 1; - filteredDivisorsSpan[filteredCount] = candidateSpan[i]; - filteredCount += add; + Span filteredSlice = filteredStorage[..filteredCount]; + ArrayView1D filteredView = candidateView.SubView(0, filteredCount); + // The compaction kernel rewrites the accepted divisors into the front of the candidate buffer, + // so copying the compacted prefix is enough to build the CPU span for cycle checks. + filteredView.CopyToCPU(ref MemoryMarshal.GetReference(filteredSlice), filteredCount); } currentK += (ulong)effectiveCount; @@ -362,8 +379,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv continue; } + Span filteredDivisorsSpan = filteredStorage[..filteredCount]; Span cycleSpan = cycleLengthStorage[..filteredCount]; - cycleCache.GetCycleLengths(filteredDivisorsSpan[..filteredCount], cycleSpan); + cycleCache.GetCycleLengths(filteredDivisorsSpan, cycleSpan); int admissibleCount = 0; for (int i = 0; i < filteredCount; i++) @@ -543,17 +561,39 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D remainder3View, byte lastIsSevenFlag, ArrayView1D maskView, - Span destination, Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> maskKernel) { if (count <= 0) { - destination.Clear(); return; } maskKernel(count, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView); - maskView.CopyToCPU(ref MemoryMarshal.GetReference(destination), count); + } + + private static int CompactCandidatesOnGpu( + int count, + ArrayView1D maskView, + ArrayView1D sourceView, + ArrayView1D destinationView, + ArrayView1D countView, + Action, ArrayView, ArrayView, ArrayView> compactionKernel) + { + if (count <= 0) + { + return 0; + } + + SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count); + Index1D gridDim = new Index1D(1); + Index1D groupDim = new Index1D(count); + KernelConfig kernelConfig = new KernelConfig(gridDim, groupDim, sharedMemoryConfig); + + compactionKernel(kernelConfig, new Index1D(count), maskView, sourceView, destinationView, countView); + + int filteredCount = 0; + countView.SubView(0, 1).CopyToCPU(ref filteredCount, 1); + return filteredCount; } private static byte FetchLastRemainder(ArrayView1D remainderView, int count) @@ -775,6 +815,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderBuffer = accelerator.Allocate1D(actualCapacity * RemainderSlotCount); + FilterCountBuffer = accelerator.Allocate1D(1); ArrayPool ulongPool = ThreadStaticPools.UlongPool; Divisors = ulongPool.Rent(actualCapacity); Exponents = ulongPool.Rent(actualCapacity); @@ -785,6 +826,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv RemainderDeltaKernel = owner.GetRemainderDeltaKernel(accelerator); RemainderScanKernel = owner.GetRemainderScanKernel(accelerator); CandidateMaskKernel = owner.GetCandidateMaskKernel(accelerator); + CandidateCompactionKernel = owner.GetCandidateCompactionKernel(accelerator); CandidateGenerationKernel = owner.GetCandidateGenerationKernel(accelerator); Capacity = actualCapacity; CycleCapacity = cycleCapacity; @@ -805,6 +847,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly MemoryBuffer1D RemainderDeltaBuffer; internal readonly MemoryBuffer1D RemainderBuffer; + internal readonly MemoryBuffer1D FilterCountBuffer; internal readonly ulong[] Divisors; @@ -822,6 +865,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, ArrayView, ArrayView> CandidateCompactionKernel; internal readonly Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> CandidateGenerationKernel; internal readonly int Capacity; @@ -838,6 +882,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorDeltaBuffer.Dispose(); RemainderDeltaBuffer.Dispose(); RemainderBuffer.Dispose(); + FilterCountBuffer.Dispose(); ArrayPool ulongPool = ThreadStaticPools.UlongPool; ArrayPool bytePool = ThreadStaticPools.BytePool; ArrayPool divisorPool = ThreadStaticPools.MontgomeryDivisorDataPool; From 586eb8d30396a6cff6054c3a17ba3e2122e78c7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 18:48:03 +0200 Subject: [PATCH 06/10] Document hot-path guards in GPU mask compaction helpers --- ...MersenneNumberDivisorByDivisorGpuTester.cs | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index a5b7955..92e4b15 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -344,10 +344,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv BuildCandidateMaskOnGpu(effectiveCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, candidateMaskKernel); - byte lastRemainder10 = FetchLastRemainder(remainder10View, effectiveCount); - byte lastRemainder8 = FetchLastRemainder(remainder8View, effectiveCount); - byte lastRemainder5 = FetchLastRemainder(remainder5View, effectiveCount); - byte lastRemainder3 = FetchLastRemainder(remainder3View, effectiveCount); + byte lastRemainder10 = effectiveCount > 0 ? FetchLastRemainder(remainder10View, effectiveCount) : (byte)0; + byte lastRemainder8 = effectiveCount > 0 ? FetchLastRemainder(remainder8View, effectiveCount) : (byte)0; + byte lastRemainder5 = effectiveCount > 0 ? FetchLastRemainder(remainder5View, effectiveCount) : (byte)0; + byte lastRemainder3 = effectiveCount > 0 ? FetchLastRemainder(remainder3View, effectiveCount) : (byte)0; remainder10 = AddModGpu(lastRemainder10, step10, 10); remainder8 = AddModGpu(lastRemainder8, step8, 8); @@ -563,10 +563,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D maskView, Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> maskKernel) { - if (count <= 0) - { - return; - } + // EvenPerfectBitScanner always supplies a positive count here, + // so the guard stays commented out to keep the hot path branch-free. + // if (count <= 0) + // { + // return; + // } maskKernel(count, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView); } @@ -579,10 +581,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D countView, Action, ArrayView, ArrayView, ArrayView> compactionKernel) { - if (count <= 0) - { - return 0; - } + // EvenPerfectBitScanner always compacts with a positive count, + // so the guard stays commented out to keep the hot path branch-free. + // if (count <= 0) + // { + // return 0; + // } SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count); Index1D gridDim = new Index1D(1); @@ -598,11 +602,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static byte FetchLastRemainder(ArrayView1D remainderView, int count) { - if (count <= 0) - { - return 0; - } - + // Callers guard against zero counts before invoking this helper to keep the hot path branch-free. byte value = 0; remainderView.SubView(count - 1, 1).CopyToCPU(ref value, 1); return value; From ea0770a84a80d0ac0ed7cfacd0a240d4efe15352 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 21:41:51 +0200 Subject: [PATCH 07/10] Clarify stream kernel bounds assumptions --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 34 +++++++++++-------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 7e6258f..e9f8a10 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -87,11 +87,13 @@ internal static class DivisorByDivisorKernels byte modulus) { int globalIndex = index; - int length = (int)remainders.Length; - if (globalIndex >= length) - { - return; - } + // EvenPerfectBitScanner launches this stream kernel with an extent matching the remainder span length, + // so the defensive guard stays commented out to keep the kernel branch-free. + // int length = (int)remainders.Length; + // if (globalIndex >= length) + // { + // return; + // } int localIndex = Group.IdxX; int groupSize = Group.Dimension.X; @@ -136,11 +138,13 @@ internal static class DivisorByDivisorKernels ArrayView mask) { int globalIndex = index; - int length = (int)mask.Length; - if (globalIndex >= length) - { - return; - } + // EvenPerfectBitScanner launches this auto-grouped kernel with an extent matching the mask length, + // so the defensive guard stays commented out to keep the kernel branch-free. + // int length = (int)mask.Length; + // if (globalIndex >= length) + // { + // return; + // } byte value10 = remainder10[globalIndex]; bool lastIsSeven = lastIsSevenFlag != 0; @@ -164,10 +168,12 @@ internal static class DivisorByDivisorKernels { int globalIndex = index; int length = (int)mask.Length; - if (globalIndex >= length) - { - return; - } + // EvenPerfectBitScanner launches this stream kernel with an extent matching the mask length, + // so the defensive guard stays commented out to keep the kernel branch-free. + // if (globalIndex >= length) + // { + // return; + // } int localIndex = Group.IdxX; int groupSize = Group.Dimension.X; From 4d9158e8f390a3d975b4d7569f647e96ec970250 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 22:01:39 +0200 Subject: [PATCH 08/10] Make GPU mask compaction branchless with shared staging --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 54 ++++++++++++------- ...MersenneNumberDivisorByDivisorGpuTester.cs | 3 +- 2 files changed, 38 insertions(+), 19 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index e9f8a10..e7546cc 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -1,4 +1,5 @@ using ILGPU; +using ILGPU.Algorithms; using ILGPU.Runtime; using PerfectNumbers.Core; @@ -179,43 +180,60 @@ internal static class DivisorByDivisorKernels int groupSize = Group.Dimension.X; var shared = SharedMemory.GetDynamic(); + int prefixBase = 0; + int mappingBase = length; + int fallbackBase = length * 2; + // The dynamic shared-memory slice stores prefix sums, compacted index mappings, and fallback slots + // in separate regions so the scatter stage can avoid divergent branches. + byte maskValue = mask[globalIndex]; int accepted = maskValue != 0 ? 1 : 0; - shared[localIndex] = accepted; - - ulong candidateValue = sourceCandidates[globalIndex]; + shared[prefixBase + localIndex] = accepted; + shared[fallbackBase + localIndex] = globalIndex; Group.Barrier(); int offset = 1; while (offset < groupSize) { - int addend = localIndex >= offset ? shared[localIndex - offset] : 0; + int addendIndex = prefixBase + localIndex - offset; + int addend = localIndex >= offset ? shared[addendIndex] : 0; Group.Barrier(); - int current = shared[localIndex]; + int currentIndex = prefixBase + localIndex; + int current = shared[currentIndex]; int sum = current + addend; - if (localIndex >= offset) - { - shared[localIndex] = sum; - } + int shouldUpdate = localIndex >= offset ? 1 : 0; + shared[currentIndex] = shouldUpdate != 0 ? sum : current; Group.Barrier(); offset <<= 1; } - if (accepted != 0) - { - int targetIndex = shared[localIndex] - 1; - compactedCandidates[targetIndex] = candidateValue; - } + int inclusivePrefix = shared[prefixBase + localIndex]; + int exclusiveIndex = inclusivePrefix - accepted; + int targetBase = accepted != 0 ? mappingBase : fallbackBase; + int targetOffset = accepted != 0 ? exclusiveIndex : localIndex; + shared[targetBase + targetOffset] = globalIndex; - if (globalIndex == length - 1) - { - compactedCount[0] = shared[localIndex]; - } + Group.Barrier(); + + int totalAccepted = shared[prefixBase + length - 1]; + int writerMask = localIndex < totalAccepted ? 1 : 0; + int lookupBase = writerMask != 0 ? mappingBase : fallbackBase; + int sourceIndex = shared[lookupBase + localIndex]; + + ulong candidateValue = sourceCandidates[sourceIndex]; + int destinationIndex = localIndex; + ulong existingValue = compactedCandidates[destinationIndex]; + compactedCandidates[destinationIndex] = writerMask != 0 ? candidateValue : existingValue; + + int lastLaneMask = globalIndex == length - 1 ? 1 : 0; + int countValue = shared[prefixBase + localIndex]; + int storedCount = lastLaneMask != 0 ? countValue : 0; + Atomic.Max(ref compactedCount[0], storedCount); } public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 92e4b15..6dd5c47 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -588,7 +588,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // return 0; // } - SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count); + countView.MemSetToZero(); + SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count * 3); Index1D gridDim = new Index1D(1); Index1D groupDim = new Index1D(count); KernelConfig kernelConfig = new KernelConfig(gridDim, groupDim, sharedMemoryConfig); From faab60c2cd78aae5825bd3fbbafb2ddf2451bb66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 23:51:19 +0200 Subject: [PATCH 09/10] Reduce host hit detection overhead --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 19 +++++++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 53 +++++++++++++++---- 2 files changed, 61 insertions(+), 11 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index e7546cc..b15a091 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -236,6 +236,25 @@ internal static class DivisorByDivisorKernels Atomic.Max(ref compactedCount[0], storedCount); } + public static void PublishFirstHitKernel( + Index1D index, + ArrayView hits, + ArrayView firstHit) + { + int globalIndex = index; + // EvenPerfectBitScanner launches this auto-grouped kernel with an extent matching the hit view length, + // so the defensive guard stays commented out to keep the hot path branch-free. + // int length = (int)hits.Length; + // if (globalIndex >= length) + // { + // return; + // } + + byte hitValue = hits[globalIndex]; + int candidate = hitValue != 0 ? globalIndex : int.MaxValue; + Atomic.Min(ref firstHit[0], candidate); + } + public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 6dd5c47..342a935 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -29,6 +29,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv 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, ArrayView, ArrayView>> _candidateCompactionKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView>> _hitDetectionKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); @@ -58,6 +59,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator, acc => acc.LoadStreamKernel, ArrayView, ArrayView, ArrayView>(DivisorByDivisorKernels.CompactCandidateMaskKernel)); + private Action, ArrayView> GetHitDetectionKernel(Accelerator accelerator) => + _hitDetectionKernelCache.GetOrAdd( + accelerator, + acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.PublishFirstHitKernel)); + private Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => _candidateGenerationKernelCache.GetOrAdd( accelerator, @@ -137,6 +143,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.DivisorsBuffer, resources.ExponentBuffer, resources.HitsBuffer, + resources.HitIndexBuffer, resources.CandidateBuffer, resources.NextStartBuffer, resources.DivisorDeltaBuffer, @@ -148,9 +155,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.CandidateMaskKernel, resources.CandidateCompactionKernel, resources.CandidateGenerationKernel, + resources.HitDetectionKernel, resources.Divisors, resources.Exponents, - resources.Hits, resources.DivisorData, resources.CycleCandidates, resources.CycleLengths, @@ -210,6 +217,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D divisorsBuffer, MemoryBuffer1D exponentBuffer, MemoryBuffer1D hitsBuffer, + MemoryBuffer1D hitIndexBuffer, MemoryBuffer1D candidateBuffer, MemoryBuffer1D nextStartBuffer, MemoryBuffer1D divisorDeltaBuffer, @@ -221,9 +229,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> candidateMaskKernel, Action, ArrayView, ArrayView, ArrayView> candidateCompactionKernel, Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> candidateGenerationKernel, + Action, ArrayView> hitDetectionKernel, in ulong[] divisors, ulong[] exponents, - byte[] hits, MontgomeryDivisorData[] divisorData, ulong[] cycleCandidates, ulong[] cycleLengths, @@ -282,7 +290,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong[] filteredDivisors = ulongPool.Rent(chunkCapacity); var viewCache = new BufferSliceCache(divisorDeltaBuffer, remainderDeltaBuffer, remainderBuffer, hitsBuffer, remainderStride); Span candidateStorage = cycleCandidates.AsSpan(); - Span hitStorage = hits.AsSpan(); Span filteredStorage = filteredDivisors.AsSpan(); Span cycleLengthStorage = cycleLengths.AsSpan(); Span divisorDataStorage = divisorData.AsSpan(); @@ -399,10 +406,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Span divisorDataSpan = divisorDataStorage[..admissibleCount]; Span divisorSpan = divisorStorage[..admissibleCount]; Span exponentSpan = exponentStorage[..admissibleCount]; - Span hitsSpan = hitStorage[..admissibleCount]; ArrayView1D divisorView = divisorsBuffer.View.SubView(0, admissibleCount); ArrayView1D exponentView = exponentBuffer.View.SubView(0, admissibleCount); ArrayView1D hitsView = hitsBuffer.View.SubView(0, admissibleCount); + ArrayView1D hitIndexView = hitIndexBuffer.View.SubView(0, 1); for (int i = 0; i < admissibleCount; i++) { @@ -415,9 +422,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), admissibleCount); exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), admissibleCount); kernel(admissibleCount, divisorView, exponentView, hitsView); - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), admissibleCount); - int hitIndex = System.MemoryExtensions.IndexOf(hitsSpan, (byte)1); + int hitIndex = FindFirstHitOnGpu(admissibleCount, hitsView, hitIndexView, hitDetectionKernel); bool hitFound = hitIndex >= 0; composite = hitFound; int lastIndex = admissibleCount - 1; @@ -601,6 +607,30 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return filteredCount; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int FindFirstHitOnGpu( + int count, + ArrayView1D hitsView, + ArrayView1D hitIndexView, + Action, ArrayView> hitKernel) + { + // EvenPerfectBitScanner only calls this helper with positive counts, + // so the guard stays commented out to keep the hot path branch-free. + // if (count <= 0) + // { + // return -1; + // } + + int sentinel = int.MaxValue; + hitIndexView.SubView(0, 1).CopyFromCPU(ref sentinel, 1); + + hitKernel(count, hitsView, hitIndexView); + + int firstHit = sentinel; + hitIndexView.SubView(0, 1).CopyToCPU(ref firstHit, 1); + return firstHit >= count ? -1 : firstHit; + } + private static byte FetchLastRemainder(ArrayView1D remainderView, int count) { // Callers guard against zero counts before invoking this helper to keep the hot path branch-free. @@ -811,6 +841,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorsBuffer = accelerator.Allocate1D(actualCapacity); ExponentBuffer = accelerator.Allocate1D(actualCapacity); HitsBuffer = accelerator.Allocate1D(actualCapacity); + HitIndexBuffer = accelerator.Allocate1D(1); CandidateBuffer = accelerator.Allocate1D(actualCapacity); NextStartBuffer = accelerator.Allocate1D(1); DivisorDeltaBuffer = accelerator.Allocate1D(actualCapacity); @@ -820,7 +851,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayPool ulongPool = ThreadStaticPools.UlongPool; Divisors = ulongPool.Rent(actualCapacity); Exponents = ulongPool.Rent(actualCapacity); - Hits = ThreadStaticPools.BytePool.Rent(actualCapacity); DivisorData = ThreadStaticPools.MontgomeryDivisorDataPool.Rent(actualCapacity); CycleCandidates = ulongPool.Rent(cycleCapacity); CycleLengths = ulongPool.Rent(cycleCapacity); @@ -829,6 +859,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv CandidateMaskKernel = owner.GetCandidateMaskKernel(accelerator); CandidateCompactionKernel = owner.GetCandidateCompactionKernel(accelerator); CandidateGenerationKernel = owner.GetCandidateGenerationKernel(accelerator); + HitDetectionKernel = owner.GetHitDetectionKernel(accelerator); Capacity = actualCapacity; CycleCapacity = cycleCapacity; } @@ -839,6 +870,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly MemoryBuffer1D HitsBuffer; + internal readonly MemoryBuffer1D HitIndexBuffer; + internal readonly MemoryBuffer1D CandidateBuffer; internal readonly MemoryBuffer1D NextStartBuffer; @@ -854,8 +887,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly ulong[] Exponents; - internal readonly byte[] Hits; - internal readonly MontgomeryDivisorData[] DivisorData; internal readonly ulong[] CycleCandidates; @@ -868,6 +899,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly Action, ArrayView, ArrayView, ArrayView, byte, ArrayView> CandidateMaskKernel; internal readonly Action, ArrayView, ArrayView, ArrayView> CandidateCompactionKernel; internal readonly Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> CandidateGenerationKernel; + internal readonly Action, ArrayView> HitDetectionKernel; internal readonly int Capacity; @@ -878,6 +910,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorsBuffer.Dispose(); ExponentBuffer.Dispose(); HitsBuffer.Dispose(); + HitIndexBuffer.Dispose(); CandidateBuffer.Dispose(); NextStartBuffer.Dispose(); DivisorDeltaBuffer.Dispose(); @@ -885,11 +918,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv RemainderBuffer.Dispose(); FilterCountBuffer.Dispose(); ArrayPool ulongPool = ThreadStaticPools.UlongPool; - ArrayPool bytePool = ThreadStaticPools.BytePool; ArrayPool divisorPool = ThreadStaticPools.MontgomeryDivisorDataPool; ulongPool.Return(Divisors, clearArray: false); ulongPool.Return(Exponents, clearArray: false); - bytePool.Return(Hits, clearArray: false); divisorPool.Return(DivisorData, clearArray: false); ulongPool.Return(CycleCandidates, clearArray: false); ulongPool.Return(CycleLengths, clearArray: false); From ea0158589246e03480ef4dfe44b4cddd1ea6f35d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 23 Oct 2025 00:34:17 +0200 Subject: [PATCH 10/10] Convert divisor session results to hits on the GPU --- .../MersenneNumberDivisorGpuTesterTests.cs | 24 +------- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 18 ++++++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 56 ++++++++++++------- 3 files changed, 56 insertions(+), 42 deletions(-) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index 29988a2..8f67d06 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -197,17 +197,7 @@ public class MersenneNumberDivisorGpuTesterTests session.CheckDivisor(7UL, divisorData, cycle, exponents, hits); - var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession) - .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!; - ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!; - Span residues = hostBuffer.AsSpan(0, exponents.Length); - byte[] computedHits = new byte[exponents.Length]; - for (int i = 0; i < residues.Length; i++) - { - computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0; - } - - computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); + hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); } [Fact] @@ -229,17 +219,7 @@ public class MersenneNumberDivisorGpuTesterTests session.CheckDivisor(11UL, divisorData, cycle, exponents, hits); - var hostBufferField = typeof(MersenneNumberDivisorByDivisorGpuTester.DivisorScanSession) - .GetField("_hostBuffer", BindingFlags.NonPublic | BindingFlags.Instance)!; - ulong[] hostBuffer = (ulong[])hostBufferField.GetValue(session)!; - Span residues = hostBuffer.AsSpan(0, exponents.Length); - byte[] computedHits = new byte[exponents.Length]; - for (int i = 0; i < residues.Length; i++) - { - computedHits[i] = residues[i] == divisorData.MontgomeryOne ? (byte)1 : (byte)0; - } - - computedHits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); + hits.Should().Equal(new byte[] { 1, 0, 1, 0, 1 }); } [Fact] diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index b15a091..1540c09 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -255,6 +255,24 @@ internal static class DivisorByDivisorKernels Atomic.Min(ref firstHit[0], candidate); } + public static void ConvertMontgomeryResultsToHitsKernel( + Index1D index, + ArrayView results, + ArrayView hits) + { + int globalIndex = index; + // EvenPerfectBitScanner launches this kernel with an extent matching the result and hit spans, + // so the defensive bounds guard stays commented out to keep the conversion branch-free. + // int length = (int)results.Length; + // if (globalIndex >= length) + // { + // return; + // } + + ulong montgomeryResult = results[globalIndex]; + hits[globalIndex] = montgomeryResult == 1UL ? (byte)1 : (byte)0; + } + public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 342a935..0af37c0 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -30,6 +30,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView, byte, ArrayView>> _candidateMaskKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView, ArrayView, ArrayView>> _candidateCompactionKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView>> _hitDetectionKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView>> _hitConversionKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary, ArrayView, ArrayView, GpuUInt128, GpuUInt128>> _candidateGenerationKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _acceleratorPool = new(); @@ -64,6 +65,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.PublishFirstHitKernel)); + private Action, ArrayView> GetResultConversionKernel(Accelerator accelerator) => + _hitConversionKernelCache.GetOrAdd( + accelerator, + acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(DivisorByDivisorKernels.ConvertMontgomeryResultsToHitsKernel)); + private Action, ArrayView, ArrayView, GpuUInt128, GpuUInt128> GetCandidateGenerationKernel(Accelerator accelerator) => _candidateGenerationKernelCache.GetOrAdd( accelerator, @@ -686,10 +692,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; - private Action, ArrayView> _kernel; - private MemoryBuffer1D _exponentsBuffer; - private MemoryBuffer1D _resultsBuffer; - private ulong[] _hostBuffer; + private Action, ArrayView> _kernel = null!; + private Action, ArrayView> _hitConversionKernel = null!; + private MemoryBuffer1D _exponentsBuffer = null!; + private MemoryBuffer1D _resultsBuffer = null!; + private MemoryBuffer1D _hitBuffer = null!; + private ulong[] _hostBuffer = null!; private int _capacity; private bool _disposed; @@ -709,8 +717,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private void EnsureExecutionResourcesLocked() { _kernel = _owner.GetKernelByPrimeExponent(_accelerator); + _hitConversionKernel = _owner.GetResultConversionKernel(_accelerator); _exponentsBuffer = _accelerator.Allocate1D(_capacity); _resultsBuffer = _accelerator.Allocate1D(_capacity); + _hitBuffer = _accelerator.Allocate1D(_capacity); _hostBuffer = ThreadStaticPools.UlongPool.Rent(_capacity); } @@ -722,26 +732,33 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } int length = primes.Length; - if (length == 0) - { - return; - } + // EvenPerfectBitScanner always supplies at least one exponent per divisor check, so the previous zero-length guard stays commented out. + // if (length == 0) + // { + // return; + // } - if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) - { - hits.Clear(); - return; - } + // The GPU divisor sessions only materialize odd moduli greater than one (q = 2kp + 1), + // so the defensive modulus guard stays commented out to keep the hot path branch-free. + // if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) + // { + // hits.Clear(); + // return; + // } Monitor.Enter(_lease.ExecutionLock); - if (_kernel == null) + if (_kernel == null || _hitConversionKernel == null) { EnsureExecutionResourcesLocked(); } + var kernel = _kernel!; + var hitConversionKernel = _hitConversionKernel!; + ArrayView1D exponentsView = _exponentsBuffer.View; ArrayView1D resultsView = _resultsBuffer.View; + ArrayView1D hitView = _hitBuffer.View; Span hostSpan = _hostBuffer.AsSpan(0, _capacity); int offset = 0; @@ -755,14 +772,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv exponentSlice.CopyFromCPU(ref MemoryMarshal.GetReference(hostSpan), batchSize); ArrayView1D resultSlice = resultsView.SubView(0, batchSize); - _kernel(batchSize, divisorData, exponentSlice, resultSlice); - resultSlice.CopyToCPU(ref MemoryMarshal.GetReference(hostSpan), batchSize); + kernel(batchSize, divisorData, exponentSlice, resultSlice); + + ArrayView1D hitSliceView = hitView.SubView(0, batchSize); + hitConversionKernel(batchSize, resultSlice, hitSliceView); Span hitSlice = hits.Slice(offset, batchSize); - for (int i = 0; i < batchSize; i++) - { - hitSlice[i] = hostSpan[i] == 1UL ? (byte)1 : (byte)0; - } + hitSliceView.CopyToCPU(ref MemoryMarshal.GetReference(hitSlice), batchSize); offset += batchSize; }