From 00341fbb8c8259a829729b78f2ec3295ad0b44de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 22 Oct 2025 03:49:04 +0200 Subject: [PATCH] Simplify GPU by-divisor configuration guards --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 95 +-- ...MersenneNumberDivisorByDivisorGpuTester.cs | 647 ++++++++++-------- PerfectNumbers.Core/ThreadStaticPools.cs | 12 + 3 files changed, 390 insertions(+), 364 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 13c6e93..068e718 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -10,14 +10,11 @@ internal static class DivisorByDivisorKernels { MontgomeryDivisorData divisor = divisors[index]; ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } - + // The GPU by-divisor pipeline only materializes odd moduli greater than one (q = 2kp + 1), + // so the defensive guard for invalid values stays disabled here to keep the kernel branch-free. ulong exponent = exponents[index]; - hits[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent) == 1UL ? (byte)1 : (byte)0; + ulong montgomeryResult = ULongExtensions.Pow2MontgomeryModWindowedGpuConvertToStandard(divisor, exponent); + hits[index] = montgomeryResult == 1UL ? (byte)1 : (byte)0; } public static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) @@ -50,42 +47,25 @@ internal static class DivisorByDivisorKernels int mod = modulus; while (offset < groupSize) { - int addend = 0; - if (localIndex >= offset) - { - addend = shared[localIndex - offset]; - } + bool apply = localIndex >= offset; + int addend = apply ? shared[localIndex - offset] : 0; Group.Barrier(); - if (localIndex >= offset) - { - int sum = shared[localIndex] + addend; - if (sum >= mod) - { - sum %= mod; - } - - shared[localIndex] = sum; - } + int current = shared[localIndex]; + int sum = apply ? current + addend : current; + int wrap = sum >= mod ? 1 : 0; + sum -= wrap * mod; + shared[localIndex] = apply ? sum : current; Group.Barrier(); offset <<= 1; } - if (globalIndex == 0) - { - remainders[0] = baseRemainder; - return; - } - int remainder = baseRemainder + shared[localIndex]; - if (remainder >= mod) - { - remainder %= mod; - } - + int remainderWrap = remainder >= mod ? 1 : 0; + remainder -= remainderWrap * mod; remainders[globalIndex] = (byte)remainder; } @@ -106,49 +86,24 @@ internal static class DivisorByDivisorKernels } byte value10 = remainder10[globalIndex]; - bool accept10; - if (lastIsSevenFlag != 0) - { - accept10 = value10 == 3 || value10 == 7 || value10 == 9; - } - else - { - accept10 = value10 == 1 || value10 == 3 || value10 == 7 || value10 == 9; - } - - if (!accept10) - { - mask[globalIndex] = 0; - return; - } - + bool lastIsSeven = lastIsSevenFlag != 0; + bool acceptSeven = value10 == 3 || value10 == 7 || value10 == 9; + bool acceptNonSeven = value10 == 1 || value10 == 3 || value10 == 7 || value10 == 9; + bool accept10 = lastIsSeven ? acceptSeven : acceptNonSeven; byte value8 = remainder8[globalIndex]; - if (value8 != 1 && value8 != 7) - { - mask[globalIndex] = 0; - return; - } - - if (remainder3[globalIndex] == 0 || remainder5[globalIndex] == 0) - { - mask[globalIndex] = 0; - return; - } - - mask[globalIndex] = 1; + bool accept8 = value8 == 1 || value8 == 7; + bool accept5 = remainder5[globalIndex] != 0; + bool accept3 = remainder3[globalIndex] != 0; + byte accepted = (byte)((accept10 && accept8 && accept3 && accept5) ? 1 : 0); + mask[globalIndex] = accepted; } public static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[index] = 0UL; - return; - } - + // Each divisor flowing through this kernel follows q = 2kp + 1 with k >= 1, so modulus is always odd and exceeds one. ulong exponent = exponents[index]; - - results[index] = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponent); + ulong montgomeryResult = ULongExtensions.Pow2MontgomeryModWindowedGpuKeepMontgomery(divisor, exponent); + results[index] = montgomeryResult; } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index cb57cfd..4fbe697 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -3,6 +3,7 @@ using System.Collections.Concurrent; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Threading; +using PerfectNumbers.Core; using ILGPU; using ILGPU.Runtime; @@ -16,9 +17,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private const int RemainderSlotCount = 4; private int _gpuBatchSize = GpuConstants.ScanBatchSize; - private readonly object _sync = new(); + // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, + // so the synchronization fields from the previous implementation remain commented out here. + // private readonly object _sync = new(); private ulong _divisorLimit; - private bool _isConfigured; + // private bool _isConfigured; private readonly ConcurrentDictionary, ArrayView, ArrayView>> _kernelCache = new(); private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); @@ -51,7 +54,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv public int GpuBatchSize { get => _gpuBatchSize; - set => _gpuBatchSize = Math.Max(1, value); + set => _gpuBatchSize = value < 1 ? 1 : value; } int IMersenneNumberDivisorByDivisorTester.BatchSize @@ -62,11 +65,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv public void ConfigureFromMaxPrime(ulong maxPrime) { - lock (_sync) - { - _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); - _isConfigured = true; - } + // EvenPerfectBitScanner configures the GPU tester once before scanning and never mutates the configuration afterwards, + // so synchronization and runtime configuration guards are unnecessary here. + // lock (_sync) + // { + // _divisorLimit = ComputeDivisorLimitFromMaxPrimeGpu(maxPrime); + // _isConfigured = true; + // } + + _divisorLimit = ComputeDivisorLimitFromMaxPrimeGpu(maxPrime); + // _isConfigured = true; } public bool IsPrime(ulong prime, out bool divisorsExhausted) @@ -74,22 +82,27 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ulong allowedMax; int batchCapacity; - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } + // EvenPerfectBitScanner only calls into this tester after configuring it once, so we can read the cached values without locking. + // lock (_sync) + // { + // if (!_isConfigured) + // { + // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + // } - allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); - batchCapacity = _gpuBatchSize; - } + // allowedMax = ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); + // batchCapacity = _gpuBatchSize; + // } - if (allowedMax < 3UL) - { - divisorsExhausted = true; - return true; - } + allowedMax = ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); + batchCapacity = _gpuBatchSize; + + // Production scans never shrink the divisor window below three, so this guard stays commented out. + // if (allowedMax < 3UL) + // { + // divisorsExhausted = true; + // return true; + // } bool composite; bool coveredRange; @@ -145,26 +158,30 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv public void PrepareCandidates(in ReadOnlySpan primes, Span allowedMaxValues) { - if (allowedMaxValues.Length < primes.Length) - { - throw new ArgumentException("allowedMaxValues span must be at least as long as primes span.", nameof(allowedMaxValues)); - } + // The scanner always supplies matching spans, so the previous validation remains commented out. + // if (allowedMaxValues.Length < primes.Length) + // { + // throw new ArgumentException("allowedMaxValues span must be at least as long as primes span.", nameof(allowedMaxValues)); + // } ulong divisorLimit; - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } + // EvenPerfectBitScanner configures the tester exactly once before preparing candidates. + // lock (_sync) + // { + // if (!_isConfigured) + // { + // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + // } - divisorLimit = _divisorLimit; - } + // divisorLimit = _divisorLimit; + // } + + divisorLimit = _divisorLimit; for (int index = 0; index < primes.Length; index++) { - allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit); + allowedMaxValues[index] = ComputeAllowedMaxDivisorGpu(primes[index], divisorLimit); } } @@ -199,48 +216,22 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv processedCount = 0UL; lastProcessed = 0UL; - Span divisorDataSpan; - Span divisorSpan; - Span exponentSpan; - Span hitsSpan; - ArrayView1D divisorView; - ArrayView1D exponentView; - ArrayView1D hitsView; - cycleCapacity = Math.Max(1, cycleCapacity); - int chunkCapacity = Math.Min(batchCapacity, cycleCapacity); int maxThreadsPerGroup = (int)accelerator.MaxNumThreadsPerGroup; - if (maxThreadsPerGroup > 0) - { - chunkCapacity = Math.Min(chunkCapacity, maxThreadsPerGroup); - } - - if (chunkCapacity <= 0) - { - chunkCapacity = cycleCapacity; - } - - ulong[] filteredDivisors = ArrayPool.Shared.Rent(chunkCapacity); - ulong[] divisorGaps = ArrayPool.Shared.Rent(chunkCapacity); + bool clampThreads = maxThreadsPerGroup > 0; + chunkCapacity = clampThreads && chunkCapacity > maxThreadsPerGroup ? maxThreadsPerGroup : chunkCapacity; + chunkCapacity = chunkCapacity <= 0 ? cycleCapacity : chunkCapacity; UInt128 twoP128 = (UInt128)prime << 1; - if (twoP128 == UInt128.Zero) - { - coveredRange = true; - return false; - } - UInt128 allowedMax128 = allowedMax; UInt128 firstDivisor128 = twoP128 + UInt128.One; - if (firstDivisor128 > allowedMax128) - { - coveredRange = true; - return false; - } - - UInt128 maxK128 = (allowedMax128 - UInt128.One) / twoP128; - if (maxK128 == UInt128.Zero) + bool invalidStride = twoP128 == UInt128.Zero; + bool outOfRange = firstDivisor128 > allowedMax128; + UInt128 numerator = allowedMax128 - UInt128.One; + UInt128 maxK128 = invalidStride || outOfRange ? UInt128.Zero : numerator / twoP128; + bool hasCandidates = !invalidStride && !outOfRange && maxK128 != UInt128.Zero; + if (!hasCandidates) { coveredRange = true; return false; @@ -262,175 +253,227 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv bool lastIsSeven = (prime & 3UL) == 3UL; byte lastIsSevenFlag = lastIsSeven ? (byte)1 : (byte)0; DivisorCycleCache cycleCache = DivisorCycleCache.Shared; - ArrayView1D remainderBaseView = remainderBuffer.View; int remainderStride = (int)divisorDeltaBuffer.Length; - if (remainderStride <= 0) - { - remainderStride = chunkCapacity; - } + remainderStride = remainderStride <= 0 ? chunkCapacity : remainderStride; - try + ArrayPool ulongPool = ThreadStaticPools.UlongPool; + ulong[] filteredDivisors = ulongPool.Rent(chunkCapacity); + ulong[] divisorGaps = ulongPool.Rent(chunkCapacity); + 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(); + Span divisorDataStorage = divisorData.AsSpan(); + Span divisorStorage = divisors.AsSpan(); + Span exponentStorage = exponents.AsSpan(); + + while (currentK <= maxK && !composite) { - while (currentK <= maxK) + int chunkCount = Math.Min(chunkCapacity, batchCapacity); + ulong remainingK = maxK - currentK + 1UL; + chunkCount = (ulong)chunkCount > remainingK ? (int)remainingK : chunkCount; + bool hasChunk = chunkCount > 0; + int effectiveCount = hasChunk ? chunkCount : 0; + if (!hasChunk) { - int chunkCount = Math.Min(chunkCapacity, batchCapacity); - ulong remainingK = maxK - currentK + 1UL; - if ((ulong)chunkCount > remainingK) - { - chunkCount = (int)remainingK; - } - - if (chunkCount <= 0) - { - break; - } - - Span candidateSpan = cycleCandidates.AsSpan(0, chunkCount); - Span gapSpan = divisorGaps.AsSpan(0, chunkCount); - // Reuse the hits array as temporary storage for the candidate mask before it stores final GPU hit flags. - Span maskSpan = hits.AsSpan(0, chunkCount); - - UInt128 localDivisor = currentDivisor128; - - for (int i = 0; i < chunkCount; i++) - { - ulong divisorValue = (ulong)localDivisor; - candidateSpan[i] = divisorValue; - - processedCount++; - lastProcessed = divisorValue; - - localDivisor += twoP128; - } - - gapSpan[0] = 0UL; - for (int i = 1; i < chunkCount; i++) - { - gapSpan[i] = candidateSpan[i] - candidateSpan[i - 1]; - } - - currentDivisor128 = localDivisor; - - ArrayView1D gapView = divisorDeltaBuffer.View.SubView(0, chunkCount); - ArrayView1D deltaView = remainderDeltaBuffer.View.SubView(0, chunkCount); - ArrayView1D remainder10View = remainderBaseView.SubView(0, chunkCount); - ArrayView1D remainder8View = remainderBaseView.SubView(remainderStride, chunkCount); - ArrayView1D remainder5View = remainderBaseView.SubView(remainderStride * 2, chunkCount); - ArrayView1D remainder3View = remainderBaseView.SubView(remainderStride * 3, chunkCount); - ArrayView1D maskView = hitsBuffer.View.SubView(0, chunkCount); - - gapView.CopyFromCPU(ref MemoryMarshal.GetReference(gapSpan), chunkCount); - - ComputeRemaindersOnGpu(chunkCount, remainder10, 10, gapView, deltaView, remainder10View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(chunkCount, remainder8, 8, gapView, deltaView, remainder8View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(chunkCount, remainder5, 5, gapView, deltaView, remainder5View, remainderDeltaKernel, remainderScanKernel); - ComputeRemaindersOnGpu(chunkCount, remainder3, 3, gapView, deltaView, remainder3View, remainderDeltaKernel, remainderScanKernel); - - BuildCandidateMaskOnGpu(chunkCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, maskSpan, candidateMaskKernel); - - byte lastRemainder10 = FetchLastRemainder(remainder10View, chunkCount); - byte lastRemainder8 = FetchLastRemainder(remainder8View, chunkCount); - byte lastRemainder5 = FetchLastRemainder(remainder5View, chunkCount); - byte lastRemainder3 = FetchLastRemainder(remainder3View, chunkCount); - - remainder10 = AddMod(lastRemainder10, step10, 10); - remainder8 = AddMod(lastRemainder8, step8, 8); - remainder5 = AddMod(lastRemainder5, step5, 5); - remainder3 = AddMod(lastRemainder3, step3, 3); - - int filteredCount = 0; - Span filteredDivisorsSpan = filteredDivisors.AsSpan(0, chunkCount); - - for (int i = 0; i < chunkCount; i++) - { - if (maskSpan[i] == 0) - { - continue; - } - - filteredDivisorsSpan[filteredCount++] = candidateSpan[i]; - } - - currentK += (ulong)chunkCount; - - if (filteredCount == 0) - { - processedAll = currentK > maxK; - continue; - } - - Span cycleSpan = cycleLengths.AsSpan(0, filteredCount); - cycleCache.GetCycleLengths(filteredDivisorsSpan[..filteredCount], cycleSpan); - - int admissibleCount = 0; - for (int i = 0; i < filteredCount; i++) - { - if (cycleSpan[i] == prime) - { - filteredDivisorsSpan[admissibleCount++] = filteredDivisorsSpan[i]; - } - } - - if (admissibleCount == 0) - { - processedAll = currentK > maxK; - continue; - } - - divisorDataSpan = divisorData.AsSpan(0, admissibleCount); - divisorSpan = divisors.AsSpan(0, admissibleCount); - exponentSpan = exponents.AsSpan(0, admissibleCount); - // Reassign the hits span to store Montgomery kernel results after the mask data has been consumed. - hitsSpan = hits.AsSpan(0, admissibleCount); - divisorView = divisorsBuffer.View.SubView(0, admissibleCount); - exponentView = exponentBuffer.View.SubView(0, admissibleCount); - hitsView = hitsBuffer.View.SubView(0, admissibleCount); - - for (int i = 0; i < admissibleCount; i++) - { - ulong divisorValue = filteredDivisorsSpan[i]; - divisorSpan[i] = divisorValue; - exponentSpan[i] = prime; - divisorDataSpan[i] = MontgomeryDivisorData.FromModulus(divisorValue); - } - - 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); - - for (int i = 0; i < admissibleCount; i++) - { - if (hitsSpan[i] != 0) - { - composite = true; - lastProcessed = divisorSpan[i]; - break; - } - } - - if (!composite) - { - lastProcessed = divisorSpan[admissibleCount - 1]; - } - - if (composite) - { - break; - } - - processedAll = currentK > maxK; + processedAll = true; + break; } + + 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; + ArrayView1D deltaView = viewCache.RemainderDelta; + ArrayView1D remainder10View = viewCache.Remainder10; + ArrayView1D remainder8View = viewCache.Remainder8; + ArrayView1D remainder5View = viewCache.Remainder5; + ArrayView1D remainder3View = viewCache.Remainder3; + ArrayView1D maskView = viewCache.Hits; + + gapView.CopyFromCPU(ref MemoryMarshal.GetReference(gapSpan), effectiveCount); + + 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); + ComputeRemaindersOnGpu(effectiveCount, remainder3, 3, gapView, deltaView, remainder3View, remainderDeltaKernel, remainderScanKernel); + + BuildCandidateMaskOnGpu(effectiveCount, remainder10View, remainder8View, remainder5View, remainder3View, lastIsSevenFlag, maskView, maskSpan, candidateMaskKernel); + + byte lastRemainder10 = FetchLastRemainder(remainder10View, effectiveCount); + byte lastRemainder8 = FetchLastRemainder(remainder8View, effectiveCount); + byte lastRemainder5 = FetchLastRemainder(remainder5View, effectiveCount); + byte lastRemainder3 = FetchLastRemainder(remainder3View, effectiveCount); + + remainder10 = AddModGpu(lastRemainder10, step10, 10); + remainder8 = AddModGpu(lastRemainder8, step8, 8); + remainder5 = AddModGpu(lastRemainder5, step5, 5); + remainder3 = AddModGpu(lastRemainder3, step3, 3); + + int filteredCount = 0; + Span filteredDivisorsSpan = filteredStorage[..effectiveCount]; + + for (int i = 0; i < effectiveCount; i++) + { + int add = maskSpan[i] == 0 ? 0 : 1; + filteredDivisorsSpan[filteredCount] = candidateSpan[i]; + filteredCount += add; + } + + currentK += (ulong)effectiveCount; + processedAll = currentK > maxK; + + if (filteredCount == 0) + { + continue; + } + + Span cycleSpan = cycleLengthStorage[..filteredCount]; + cycleCache.GetCycleLengths(filteredDivisorsSpan[..filteredCount], cycleSpan); + + int admissibleCount = 0; + for (int i = 0; i < filteredCount; i++) + { + int add = cycleSpan[i] == prime ? 1 : 0; + filteredDivisorsSpan[admissibleCount] = filteredDivisorsSpan[i]; + admissibleCount += add; + } + + if (admissibleCount == 0) + { + continue; + } + + 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); + + for (int i = 0; i < admissibleCount; i++) + { + ulong divisorValue = filteredDivisorsSpan[i]; + divisorSpan[i] = divisorValue; + exponentSpan[i] = prime; + divisorDataSpan[i] = MontgomeryDivisorData.FromModulus(divisorValue); + } + + 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); + bool hitFound = hitIndex >= 0; + composite = hitFound; + int lastIndex = admissibleCount - 1; + lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex]; } - finally - { - ArrayPool.Shared.Return(filteredDivisors, clearArray: false); - ArrayPool.Shared.Return(divisorGaps, clearArray: false); - } + + ulongPool.Return(filteredDivisors, clearArray: false); + ulongPool.Return(divisorGaps, clearArray: false); coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); return composite; } + + private ref struct BufferSliceCache + { + private readonly MemoryBuffer1D _divisorDeltaBuffer; + private readonly MemoryBuffer1D _remainderDeltaBuffer; + private readonly MemoryBuffer1D _remainderBuffer; + private readonly MemoryBuffer1D _hitsBuffer; + private readonly int _remainderStride; + private int _length; + private ArrayView1D _divisorDelta; + private ArrayView1D _remainderDelta; + private ArrayView1D _remainder10; + private ArrayView1D _remainder8; + private ArrayView1D _remainder5; + private ArrayView1D _remainder3; + private ArrayView1D _hits; + + internal BufferSliceCache( + MemoryBuffer1D divisorDeltaBuffer, + MemoryBuffer1D remainderDeltaBuffer, + MemoryBuffer1D remainderBuffer, + MemoryBuffer1D hitsBuffer, + int remainderStride) + { + _divisorDeltaBuffer = divisorDeltaBuffer; + _remainderDeltaBuffer = remainderDeltaBuffer; + _remainderBuffer = remainderBuffer; + _hitsBuffer = hitsBuffer; + _remainderStride = remainderStride; + _length = -1; + _divisorDelta = default; + _remainderDelta = default; + _remainder10 = default; + _remainder8 = default; + _remainder5 = default; + _remainder3 = default; + _hits = default; + } + + internal void Ensure(int length) + { + if (_length == length) + { + return; + } + + _length = length; + ArrayView1D remainderBaseView = _remainderBuffer.View; + _divisorDelta = _divisorDeltaBuffer.View.SubView(0, length); + _remainderDelta = _remainderDeltaBuffer.View.SubView(0, length); + _remainder10 = remainderBaseView.SubView(0, length); + _remainder8 = remainderBaseView.SubView(_remainderStride, length); + _remainder5 = remainderBaseView.SubView(_remainderStride * 2, length); + _remainder3 = remainderBaseView.SubView(_remainderStride * 3, length); + _hits = _hitsBuffer.View.SubView(0, length); + } + + internal ArrayView1D DivisorDelta => _divisorDelta; + + internal ArrayView1D RemainderDelta => _remainderDelta; + + internal ArrayView1D Remainder10 => _remainder10; + + internal ArrayView1D Remainder8 => _remainder8; + + internal ArrayView1D Remainder5 => _remainder5; + + internal ArrayView1D Remainder3 => _remainder3; + + internal ArrayView1D Hits => _hits; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ComputeRemaindersOnGpu( int count, @@ -490,28 +533,24 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return value; } - private static byte AddMod(byte value, byte delta, byte modulus) + private static byte AddModGpu(byte value, byte delta, byte modulus) { int result = value + delta; - if (result >= modulus) - { - result -= modulus; - if (result >= modulus) - { - result %= modulus; - } - } - + int firstWrap = result >= modulus ? 1 : 0; + result -= firstWrap * modulus; + int secondWrap = result >= modulus ? 1 : 0; + result -= secondWrap * modulus; return (byte)result; } - private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) + private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime) { - if (maxPrime <= 1UL) - { - return 0UL; - } - + // The --mersenne=bydivisor flow in EvenPerfectBitScanner only calls ConfigureFromMaxPrime with primes greater than 1, + // so the guard below never trips in production runs. + // if (maxPrime <= 1UL) + // { + // return 0UL; + // } if (maxPrime - 1UL >= 64UL) { return ulong.MaxValue; @@ -520,19 +559,20 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return (1UL << (int)(maxPrime - 1UL)) - 1UL; } - private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) + private static ulong ComputeAllowedMaxDivisorGpu(ulong prime, ulong divisorLimit) { - if (prime <= 1UL) - { - return 0UL; - } - + // Production --mersenne=bydivisor runs only pass prime exponents, so the guard below never executes outside tests. + // if (prime <= 1UL) + // { + // return 0UL; + // } if (prime - 1UL >= 64UL) { return divisorLimit; } - return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); + ulong computedLimit = (1UL << (int)(prime - 1UL)) - 1UL; + return computedLimit < divisorLimit ? computedLimit : divisorLimit; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -565,7 +605,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _kernel = _owner.GetKernelByPrimeExponent(_accelerator); _exponentsBuffer = _accelerator.Allocate1D(_capacity); _resultsBuffer = _accelerator.Allocate1D(_capacity); - _hostBuffer = ArrayPool.Shared.Rent(_capacity); + _hostBuffer = ThreadStaticPools.UlongPool.Rent(_capacity); } public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, in ReadOnlySpan primes, Span hits) @@ -640,15 +680,18 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { get { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } + // The EvenPerfectBitScanner driver configures the GPU tester before exposing it to callers, so the previous synchronization guard remains commented out. + // lock (_sync) + // { + // if (!_isConfigured) + // { + // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + // } - return _divisorLimit; - } + // return _divisorLimit; + // } + + return _divisorLimit; } } @@ -695,12 +738,13 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderDeltaBuffer = accelerator.Allocate1D(actualCapacity); RemainderBuffer = accelerator.Allocate1D(actualCapacity * RemainderSlotCount); - Divisors = ArrayPool.Shared.Rent(actualCapacity); - Exponents = ArrayPool.Shared.Rent(actualCapacity); - Hits = ArrayPool.Shared.Rent(actualCapacity); - DivisorData = ArrayPool.Shared.Rent(actualCapacity); - CycleCandidates = ArrayPool.Shared.Rent(cycleCapacity); - CycleLengths = ArrayPool.Shared.Rent(cycleCapacity); + 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); RemainderDeltaKernel = owner.GetRemainderDeltaKernel(accelerator); RemainderScanKernel = owner.GetRemainderScanKernel(accelerator); CandidateMaskKernel = owner.GetCandidateMaskKernel(accelerator); @@ -749,12 +793,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorDeltaBuffer.Dispose(); RemainderDeltaBuffer.Dispose(); RemainderBuffer.Dispose(); - ArrayPool.Shared.Return(Divisors, clearArray: false); - ArrayPool.Shared.Return(Exponents, clearArray: false); - ArrayPool.Shared.Return(Hits, clearArray: false); - ArrayPool.Shared.Return(DivisorData, clearArray: false); - ArrayPool.Shared.Return(CycleCandidates, clearArray: false); - ArrayPool.Shared.Return(CycleLengths, clearArray: false); + 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); } } @@ -769,34 +816,46 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv public ulong GetAllowedMaxDivisor(ulong prime) { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } + // The configuration is immutable after setup, so the cached limit can be read without locking. + // lock (_sync) + // { + // if (!_isConfigured) + // { + // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + // } - return ComputeAllowedMaxDivisor(prime, _divisorLimit); - } + // return ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); + // } + + return ComputeAllowedMaxDivisorGpu(prime, _divisorLimit); } public DivisorScanSession CreateDivisorSession() { - lock (_sync) + // The tester is configured once at startup and the session pool relies on thread-safe collections, so the previous synchronization guard stays commented out. + // lock (_sync) + // { + // if (!_isConfigured) + // { + // throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + // } + + // if (_sessionPool.TryTake(out DivisorScanSession? session)) + // { + // session.Reset(); + // return session; + // } + + // return new DivisorScanSession(this); + // } + + if (_sessionPool.TryTake(out DivisorScanSession? session)) { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - if (_sessionPool.TryTake(out DivisorScanSession? session)) - { - session.Reset(); - return session; - } - - return new DivisorScanSession(this); + session.Reset(); + return session; } + + return new DivisorScanSession(this); } IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() => CreateDivisorSession(); diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index f8b07cc..2da421e 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -90,6 +90,18 @@ namespace PerfectNumbers.Core } } + [ThreadStatic] + private static ArrayPool? _montgomeryDivisorDataPool; + + public static ArrayPool MontgomeryDivisorDataPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _montgomeryDivisorDataPool ??= ArrayPool.Create(); + } + } + [ThreadStatic] private static ArrayPool? _intPool;