diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 10690e6..5a35e58 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -22,6 +22,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly ConcurrentDictionary, ArrayView, ArrayView>> _kernelCache = new(); private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); + private readonly ConcurrentDictionary, ArrayView, byte>> _remainderDeltaKernelCache = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView, byte, byte>> _remainderScanKernelCache = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _sessionPool = new(); @@ -31,6 +33,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeMontgomeryExponentKernel)); + private Action, ArrayView, byte> GetRemainderDeltaKernel(Accelerator accelerator) => + _remainderDeltaKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, byte>(ComputeRemainderDeltasKernel)); + + private Action, ArrayView, byte, byte> GetRemainderScanKernel(Accelerator accelerator) => + _remainderScanKernelCache.GetOrAdd( + accelerator, + acc => acc.LoadStreamKernel, ArrayView, byte, byte>(AccumulateRemaindersKernel)); + public int GpuBatchSize { get => _gpuBatchSize; @@ -94,6 +104,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.DivisorsBuffer, resources.ExponentBuffer, resources.HitsBuffer, + resources.DivisorDeltaBuffer, + resources.RemainderDeltaBuffer, + resources.RemainderBuffer, + resources.RemainderDeltaKernel, + resources.RemainderScanKernel, resources.Divisors, resources.Exponents, resources.Hits, @@ -152,6 +167,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D divisorsBuffer, MemoryBuffer1D exponentBuffer, MemoryBuffer1D hitsBuffer, + MemoryBuffer1D divisorDeltaBuffer, + MemoryBuffer1D remainderDeltaBuffer, + MemoryBuffer1D remainderBuffer, + Action, ArrayView, byte> remainderDeltaKernel, + Action, ArrayView, byte, byte> remainderScanKernel, in ulong[] divisors, ulong[] exponents, byte[] hits, @@ -163,7 +183,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv int batchCapacity = (int)divisorsBuffer.Length; bool composite = false; bool processedAll = false; - ulong currentDivisor = 3UL; processedCount = 0UL; lastProcessed = 0UL; @@ -177,102 +196,209 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorCycleCache cycleCache = DivisorCycleCache.Shared; int cycleBatchCapacity = Math.Min(batchCapacity, cycleCache.PreferredBatchSize); + int maxKernelGroupSize = (int)accelerator.MaxNumThreadsPerGroup; + if (cycleBatchCapacity > maxKernelGroupSize) + { + cycleBatchCapacity = maxKernelGroupSize; + } + if (cycleBatchCapacity <= 0) { cycleBatchCapacity = batchCapacity; } - ulong[] cycleDivisors = ArrayPool.Shared.Rent(cycleBatchCapacity); - ulong[] cycleLengths = ArrayPool.Shared.Rent(cycleBatchCapacity); + ulong[] candidateDivisors = ArrayPool.Shared.Rent(cycleBatchCapacity); + ulong[] filteredDivisors = ArrayPool.Shared.Rent(cycleBatchCapacity); + ulong[] filteredCycles = ArrayPool.Shared.Rent(cycleBatchCapacity); + ulong[] divisorGaps = ArrayPool.Shared.Rent(cycleBatchCapacity); + byte[] remainder10Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); + byte[] remainder8Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); + byte[] remainder5Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); + byte[] remainder3Buffer = ArrayPool.Shared.Rent(cycleBatchCapacity); + + 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) + { + coveredRange = true; + return false; + } + + ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; + ulong currentK = 1UL; + + byte step10 = (byte)(twoP128 % 10UL); + byte step8 = (byte)(twoP128 % 8UL); + byte step5 = (byte)(twoP128 % 5UL); + byte step3 = (byte)(twoP128 % 3UL); + + UInt128 currentDivisor128 = firstDivisor128; + byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL)); + byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL)); + byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL)); + byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL)); + bool lastIsSeven = (prime & 3UL) == 3UL; try { - while (currentDivisor <= allowedMax) + while (currentK <= maxK) { - int batchSize = 0; - bool reachedEndInBatch = false; - - while (batchSize < batchCapacity && currentDivisor <= allowedMax) + int chunkCount = Math.Min(cycleBatchCapacity, batchCapacity); + ulong remainingK = maxK - currentK + 1UL; + if ((ulong)chunkCount > remainingK) { - int chunkCount = Math.Min(cycleBatchCapacity, batchCapacity - batchSize); - ulong remainingWidth = allowedMax - currentDivisor; - ulong maxOddCount = (remainingWidth >> 1) + 1UL; - if (maxOddCount < (ulong)chunkCount) + chunkCount = (int)remainingK; + } + + if (chunkCount <= 0) + { + break; + } + + Span candidateSpan = candidateDivisors.AsSpan(0, chunkCount); + Span gapSpan = divisorGaps.AsSpan(0, chunkCount); + Span rem10Span = remainder10Buffer.AsSpan(0, chunkCount); + Span rem8Span = remainder8Buffer.AsSpan(0, chunkCount); + Span rem5Span = remainder5Buffer.AsSpan(0, chunkCount); + Span rem3Span = remainder3Buffer.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 remainderView = remainderBuffer.View.SubView(0, chunkCount); + + gapView.CopyFromCPU(ref MemoryMarshal.GetReference(gapSpan), chunkCount); + + ComputeRemaindersOnGpu(accelerator, chunkCount, remainder10, 10, gapView, deltaView, remainderView, remainderDeltaKernel, remainderScanKernel, rem10Span); + ComputeRemaindersOnGpu(accelerator, chunkCount, remainder8, 8, gapView, deltaView, remainderView, remainderDeltaKernel, remainderScanKernel, rem8Span); + ComputeRemaindersOnGpu(accelerator, chunkCount, remainder5, 5, gapView, deltaView, remainderView, remainderDeltaKernel, remainderScanKernel, rem5Span); + ComputeRemaindersOnGpu(accelerator, chunkCount, remainder3, 3, gapView, deltaView, remainderView, remainderDeltaKernel, remainderScanKernel, rem3Span); + + remainder10 = AddMod(rem10Span[chunkCount - 1], step10, 10); + remainder8 = AddMod(rem8Span[chunkCount - 1], step8, 8); + remainder5 = AddMod(rem5Span[chunkCount - 1], step5, 5); + remainder3 = AddMod(rem3Span[chunkCount - 1], step3, 3); + + int filteredCount = 0; + Span filteredDivisorsSpan = filteredDivisors.AsSpan(0, chunkCount); + Span filteredCyclesSpan = filteredCycles.AsSpan(0, chunkCount); + + for (int i = 0; i < chunkCount; i++) + { + byte r10 = rem10Span[i]; + if (lastIsSeven) { - chunkCount = (int)maxOddCount; - } - - if (chunkCount <= 0) - { - reachedEndInBatch = true; - break; - } - - Span chunkDivisors = cycleDivisors.AsSpan(0, chunkCount); - Span chunkCycles = cycleLengths.AsSpan(0, chunkCount); - - ulong fillDivisor = currentDivisor; - for (int i = 0; i < chunkCount; i++) - { - chunkDivisors[i] = fillDivisor; - fillDivisor += 2UL; - } - - cycleCache.GetCycleLengths(chunkDivisors, chunkCycles); - - for (int i = 0; i < chunkCount && batchSize < batchCapacity; i++) - { - ulong divisorValue = chunkDivisors[i]; - ulong cycleLength = chunkCycles[i]; - - processedCount++; - bool includeDivisor = true; - ulong exponent = prime; - - if (cycleLength > 0UL) + if (r10 != 3 && r10 != 7 && r10 != 9) { - ulong remainder = prime; - if (remainder >= cycleLength) - { - remainder %= cycleLength; - } - - if (remainder != 0UL) - { - includeDivisor = false; - } - else - { - exponent = 0UL; - } - } - - if (includeDivisor) - { - divisors[batchSize] = divisorValue; - exponents[batchSize] = exponent; - batchSize++; - } - - lastProcessed = divisorValue; - currentDivisor = divisorValue + 2UL; - - if (currentDivisor > allowedMax) - { - reachedEndInBatch = true; - break; + continue; } } + else + { + if (r10 != 1 && r10 != 3 && r10 != 7 && r10 != 9) + { + continue; + } + } + + byte r8 = rem8Span[i]; + if (r8 != 1 && r8 != 7) + { + continue; + } + + if (rem3Span[i] == 0 || rem5Span[i] == 0) + { + continue; + } + + filteredDivisorsSpan[filteredCount++] = candidateSpan[i]; + } + + currentK += (ulong)chunkCount; + + if (filteredCount == 0) + { + processedAll = currentK > maxK; + continue; + } + + cycleCache.GetCycleLengths(filteredDivisorsSpan[..filteredCount], filteredCyclesSpan[..filteredCount]); + + int batchSize = 0; + for (int i = 0; i < filteredCount; i++) + { + ulong divisorValue = filteredDivisorsSpan[i]; + ulong cycleLength = filteredCyclesSpan[i]; + + bool includeDivisor = true; + ulong exponent = prime; + + if (cycleLength > 0UL) + { + ulong remainder = prime; + if (remainder >= cycleLength) + { + remainder %= cycleLength; + } + + if (remainder != 0UL) + { + includeDivisor = false; + } + else + { + exponent = 0UL; + } + } + + if (!includeDivisor) + { + continue; + } + + divisors[batchSize] = divisorValue; + exponents[batchSize] = exponent; + batchSize++; } if (batchSize == 0) { - if (reachedEndInBatch) - { - processedAll = true; - break; - } - + processedAll = currentK > maxK; continue; } @@ -314,22 +440,68 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv break; } - if (reachedEndInBatch) - { - processedAll = true; - break; - } + processedAll = currentK > maxK; } } finally { - ArrayPool.Shared.Return(cycleDivisors, clearArray: false); - ArrayPool.Shared.Return(cycleLengths, clearArray: false); + ArrayPool.Shared.Return(candidateDivisors, clearArray: false); + ArrayPool.Shared.Return(filteredDivisors, clearArray: false); + ArrayPool.Shared.Return(filteredCycles, clearArray: false); + ArrayPool.Shared.Return(divisorGaps, clearArray: false); + ArrayPool.Shared.Return(remainder10Buffer, clearArray: false); + ArrayPool.Shared.Return(remainder8Buffer, clearArray: false); + ArrayPool.Shared.Return(remainder5Buffer, clearArray: false); + ArrayPool.Shared.Return(remainder3Buffer, clearArray: false); } - coveredRange = composite || processedAll || currentDivisor > allowedMax; + coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); return composite; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void ComputeRemaindersOnGpu( + Accelerator accelerator, + int count, + byte baseRemainder, + byte modulus, + ArrayView1D gapView, + ArrayView1D deltaView, + ArrayView1D remainderView, + Action, ArrayView, byte> deltaKernel, + Action, ArrayView, byte, byte> scanKernel, + Span destination) + { + if (count <= 0) + { + return; + } + + deltaKernel(count, gapView, deltaView, modulus); + SharedMemoryConfig sharedMemoryConfig = SharedMemoryConfig.RequestDynamic(count); + Index1D gridDim = new Index1D(1); + Index1D groupDim = new Index1D(count); + KernelConfig kernelConfig = new KernelConfig(gridDim, groupDim, sharedMemoryConfig); + + scanKernel(kernelConfig, new Index1D(count), deltaView, remainderView, baseRemainder, modulus); + accelerator.Synchronize(); + remainderView.CopyToCPU(ref MemoryMarshal.GetReference(destination), count); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod(byte value, byte delta, byte modulus) + { + int result = value + delta; + if (result >= modulus) + { + result -= modulus; + if (result >= modulus) + { + result %= modulus; + } + } + + return (byte)result; + } private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) { @@ -581,6 +753,75 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } + private static void ComputeRemainderDeltasKernel(Index1D index, ArrayView gaps, ArrayView deltas, byte modulus) + { + deltas[index] = (byte)(gaps[index] % modulus); + } + + private static void AccumulateRemaindersKernel( + Index1D index, + ArrayView deltas, + ArrayView remainders, + byte baseRemainder, + byte modulus) + { + int globalIndex = index; + int length = (int)remainders.Length; + if (globalIndex >= length) + { + return; + } + + int localIndex = Group.IdxX; + int groupSize = Group.Dimension.X; + var shared = SharedMemory.GetDynamic(); + + shared[localIndex] = deltas[globalIndex]; + Group.Barrier(); + + int offset = 1; + int mod = modulus; + while (offset < groupSize) + { + int addend = 0; + if (localIndex >= offset) + { + addend = shared[localIndex - offset]; + } + + Group.Barrier(); + + if (localIndex >= offset) + { + int sum = shared[localIndex] + addend; + if (sum >= mod) + { + sum %= mod; + } + + shared[localIndex] = sum; + } + + Group.Barrier(); + + offset <<= 1; + } + + if (globalIndex == 0) + { + remainders[0] = baseRemainder; + return; + } + + int remainder = baseRemainder + shared[localIndex]; + if (remainder >= mod) + { + remainder %= mod; + } + + remainders[globalIndex] = (byte)remainder; + } + private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus; @@ -619,22 +860,27 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return resources; } - return new BatchResources(accelerator, capacity); + return new BatchResources(this, accelerator, capacity); } private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) => _resourcePools.GetOrAdd(accelerator, static _ => []).Add(resources); private sealed class BatchResources : IDisposable { - internal BatchResources(Accelerator accelerator, int capacity) + internal BatchResources(MersenneNumberDivisorByDivisorGpuTester owner, Accelerator accelerator, int capacity) { DivisorsBuffer = accelerator.Allocate1D(capacity); ExponentBuffer = accelerator.Allocate1D(capacity); HitsBuffer = accelerator.Allocate1D(capacity); + DivisorDeltaBuffer = accelerator.Allocate1D(capacity); + RemainderDeltaBuffer = accelerator.Allocate1D(capacity); + RemainderBuffer = accelerator.Allocate1D(capacity); Divisors = ArrayPool.Shared.Rent(capacity); Exponents = ArrayPool.Shared.Rent(capacity); Hits = ArrayPool.Shared.Rent(capacity); DivisorData = ArrayPool.Shared.Rent(capacity); + RemainderDeltaKernel = owner.GetRemainderDeltaKernel(accelerator); + RemainderScanKernel = owner.GetRemainderScanKernel(accelerator); Capacity = capacity; } @@ -644,6 +890,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly MemoryBuffer1D HitsBuffer; + internal readonly MemoryBuffer1D DivisorDeltaBuffer; + + internal readonly MemoryBuffer1D RemainderDeltaBuffer; + + internal readonly MemoryBuffer1D RemainderBuffer; + internal readonly ulong[] Divisors; internal readonly ulong[] Exponents; @@ -652,6 +904,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal readonly MontgomeryDivisorData[] DivisorData; + internal readonly Action, ArrayView, byte> RemainderDeltaKernel; + + internal readonly Action, ArrayView, byte, byte> RemainderScanKernel; + internal readonly int Capacity; public void Dispose() @@ -659,6 +915,9 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv DivisorsBuffer.Dispose(); ExponentBuffer.Dispose(); HitsBuffer.Dispose(); + DivisorDeltaBuffer.Dispose(); + RemainderDeltaBuffer.Dispose(); + RemainderBuffer.Dispose(); ArrayPool.Shared.Return(Divisors, clearArray: false); ArrayPool.Shared.Return(Exponents, clearArray: false); ArrayPool.Shared.Return(Hits, clearArray: false); diff --git a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs index 806ac5b..a332008 100644 --- a/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs +++ b/PerfectNumbers.Core/MersenneNumberDivisorByDivisorTester.cs @@ -2,734 +2,175 @@ using System; using System.Buffers; using System.Collections.Generic; using System.Globalization; -using System.Numerics; -using System.Runtime.CompilerServices; -using System.Threading; -using System.Threading.Tasks; using Open.Numeric.Primes; namespace PerfectNumbers.Core; public static class MersenneNumberDivisorByDivisorTester { - private const int ByDivisorStateActive = 0; - private const int ByDivisorStateComposite = 1; - private const int ByDivisorStateCompleted = 2; - private const int ByDivisorStateCompletedDetailed = 3; - private const int DivisorAllocationBlockSize = 64; - - public static void Run( - List candidates, - IMersenneNumberDivisorByDivisorTester tester, - Dictionary? previousResults, - ulong startPrime, - Action markComposite, - Action clearComposite, - Action printResult, - int threadCount) + public static void Run( + List candidates, + IMersenneNumberDivisorByDivisorTester tester, + Dictionary? previousResults, + ulong startPrime, + Action markComposite, + Action clearComposite, + Action printResult, + int threadCount) + { + if (candidates.Count == 0) { - if (candidates.Count == 0) - { - Console.WriteLine("No candidates were provided for --mersenne=bydivisor."); - return; - } - - bool applyStartPrime = startPrime > 0UL; - int skippedByPreviousResults = 0; - List states = new(candidates.Count); - ulong maxPrime = 0UL; - - using IEnumerator primeEnumerator = Prime.Numbers.GetEnumerator(); - // TODO: Swap this per-run enumerator creation for the shared PrimeIterator so by-divisor scans reuse the - // benchmarked zero-allocation prime sequence instead of instantiating a new iterator for every invocation. - bool hasPrime = primeEnumerator.MoveNext(); - ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL; - - for (int index = 0; index < candidates.Count; index++) - { - ulong candidate = candidates[index]; - - if (applyStartPrime && candidate < startPrime) - { - continue; - } - - if (previousResults is not null && previousResults.ContainsKey(candidate)) - { - skippedByPreviousResults++; - continue; - } - - if (candidate <= 1UL) - { - markComposite(); - printResult(candidate, false, false, false); - continue; - } - - while (hasPrime && currentPrime < candidate) - { - hasPrime = primeEnumerator.MoveNext(); - if (!hasPrime) - { - break; - } - - currentPrime = primeEnumerator.Current; - } - - if (!hasPrime || currentPrime != candidate) - { - markComposite(); - printResult(candidate, false, false, false); - continue; - } - - if (candidate > maxPrime) - { - maxPrime = candidate; - } - - states.Add(new ByDivisorPrimeState - { - Prime = candidate, - }); - } - - if (skippedByPreviousResults > 0) - { - Console.WriteLine($"Skipped {skippedByPreviousResults.ToString(CultureInfo.InvariantCulture)} candidates excluded by previous results."); - } - - if (states.Count == 0) - { - if (applyStartPrime) - { - Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); - } - else - { - Console.WriteLine("No prime candidates remain for --mersenne=bydivisor after filtering."); - } - - return; - } - - if (maxPrime <= 1UL) - { - Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor."); - return; - } - - tester.ConfigureFromMaxPrime(maxPrime); - - int stateCount = states.Count; - ulong[] primeBatch = ArrayPool.Shared.Rent(stateCount); - ulong[] allowedMaxBatch = ArrayPool.Shared.Rent(stateCount); - - try - { - Span primeBatchSpan = primeBatch.AsSpan(0, stateCount); - Span allowedMaxSpan = allowedMaxBatch.AsSpan(0, stateCount); - - for (int batchIndex = 0; batchIndex < stateCount; batchIndex++) - { - primeBatchSpan[batchIndex] = states[batchIndex].Prime; - } - - tester.PrepareCandidates(primeBatchSpan, allowedMaxSpan); - - List filteredStates = new(stateCount); - for (int stateIndex = 0; stateIndex < stateCount; stateIndex++) - { - ulong prime = primeBatchSpan[stateIndex]; - ulong allowedMax = allowedMaxSpan[stateIndex]; - if (allowedMax < 3UL) - { - clearComposite(); - printResult(prime, true, true, true); - continue; - } - - filteredStates.Add(new ByDivisorPrimeState - { - Prime = prime, - AllowedMax = allowedMax, - Completed = false, - Composite = false, - DetailedCheck = false, - }); - } - - states = filteredStates; - } - finally - { - ArrayPool.Shared.Return(allowedMaxBatch, clearArray: true); - ArrayPool.Shared.Return(primeBatch, clearArray: true); - } - - if (states.Count == 0) - { - if (applyStartPrime) - { - Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); - } - - return; - } - - stateCount = states.Count; - ulong[] primeValues = new ulong[stateCount]; - ulong[] allowedMaxValues = new ulong[stateCount]; - int[] stateFlags = new int[stateCount]; - - for (int stateIndex = 0; stateIndex < stateCount; stateIndex++) - { - ByDivisorPrimeState currentState = states[stateIndex]; - primeValues[stateIndex] = currentState.Prime; - allowedMaxValues[stateIndex] = currentState.AllowedMax; - stateFlags[stateIndex] = ByDivisorStateActive; - } - - states.Clear(); - - ulong divisorLimit = tester.DivisorLimit; - ulong nextDivisor = 3UL; - long finalDivisorBits = unchecked((long)3UL); - int divisorsExhaustedFlag = 0; - int finalizerState = 0; - int finalizationCompleted = 0; - int remainingStates = stateCount; - int activeStartIndex = 0; - long[] activeStateMask = new long[(stateCount + 63) >> 6]; - - for (int maskStateIndex = 0; maskStateIndex < stateCount; maskStateIndex++) - { - int wordIndex = maskStateIndex >> 6; - int bitIndex = maskStateIndex & 63; - activeStateMask[wordIndex] |= 1L << bitIndex; - } - - Task[] workers = new Task[Math.Max(1, threadCount)]; - - for (int workerIndex = 0; workerIndex < workers.Length; workerIndex++) - { - int capturedStateCount = stateCount; - // TODO: Swap Task.Factory.StartNew for the shared low-overhead work queue once the scanner's - // scheduler is exposed so divisor scans avoid the extra allocations and context switches that - // benchmarks showed with Task-based fan-out. - workers[workerIndex] = Task.Factory.StartNew(() => - { - using var session = tester.CreateDivisorSession(); - byte[] hitsBuffer = ArrayPool.Shared.Rent(capturedStateCount); - ulong[] primeBuffer = ArrayPool.Shared.Rent(capturedStateCount); - int[] indexBuffer = ArrayPool.Shared.Rent(capturedStateCount); - PendingResult[] completionsBuffer = ArrayPool.Shared.Rent(capturedStateCount); - PendingResult[] compositesBuffer = ArrayPool.Shared.Rent(capturedStateCount); - int completionsCount = 0; - int compositesCount = 0; - ulong localDivisorCursor = 0UL; - int localDivisorsRemaining = 0; - bool exhausted = false; - ulong divisor; - int activeCount; - Span hitsSpan = default; - int hitIndex = 0; - int index; - Dictionary? factorCache = null; - ulong divisorCycle = 0UL; - - try - { - while (true) - { - if (Volatile.Read(ref remainingStates) == 0) - { - if (Volatile.Read(ref finalizationCompleted) != 0 || Volatile.Read(ref divisorsExhaustedFlag) == 0) - { - break; - } - } - - divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted, ref localDivisorCursor, ref localDivisorsRemaining); - if (divisor != 0UL) - { - - activeCount = BuildPrimeBuffer( - divisor, - primeValues, - allowedMaxValues, - stateFlags, - primeBuffer, - indexBuffer, - completionsBuffer, - ref completionsCount, - ref remainingStates, - activeStateMask, - ref activeStartIndex, - markComposite, - clearComposite, - printResult); - - if (completionsCount > 0) - { - FlushPendingResults(completionsBuffer, ref completionsCount, markComposite, clearComposite, printResult); - } - - if (activeCount == 0) - { - continue; - } - divisorCycle = 0UL; - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); - } - else - { - factorCache ??= new Dictionary(16); - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(divisor, primeBuffer[0], factorCache, out divisorCycle) || divisorCycle == 0UL) - { - divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); - } - } - - hitsSpan = hitsBuffer.AsSpan(0, activeCount); - hitsSpan.Clear(); - session.CheckDivisor(divisor, divisorCycle, primeBuffer.AsSpan(0, activeCount), hitsSpan); - - for (hitIndex = 0; hitIndex < activeCount; hitIndex++) - { - if (hitsSpan[hitIndex] == 0) - { - continue; - } - - index = indexBuffer[hitIndex]; - if (Interlocked.CompareExchange(ref stateFlags[index], ByDivisorStateComposite, ByDivisorStateActive) == ByDivisorStateActive) - { - ClearActiveMask(activeStateMask, index); - Interlocked.Decrement(ref remainingStates); - compositesBuffer[compositesCount++] = new PendingResult(primeValues[index], detailedCheck: true, passedAllTests: false); - } - - if (compositesCount == compositesBuffer.Length) - { - FlushPendingResults(compositesBuffer, ref compositesCount, markComposite, clearComposite, printResult); - } - } - - if (compositesCount > 0) - { - FlushPendingResults(compositesBuffer, ref compositesCount, markComposite, clearComposite, printResult); - } - - continue; - } - if (!exhausted) - { - if (Volatile.Read(ref remainingStates) == 0) - { - break; - } - - Thread.Yield(); - continue; - } - - if (Interlocked.CompareExchange(ref finalizerState, 1, 0) == 0) - { - FinalizeRemainingStates( - primeValues, - allowedMaxValues, - stateFlags, - ref remainingStates, - ref finalDivisorBits, - completionsBuffer, - ref completionsCount, - activeStateMask, - clearComposite, - printResult, - markComposite); - - if (completionsCount > 0) - { - FlushPendingResults(completionsBuffer, ref completionsCount, markComposite, clearComposite, printResult); - } - - Volatile.Write(ref finalizationCompleted, 1); - } - else - { - while (Volatile.Read(ref finalizationCompleted) == 0 && Volatile.Read(ref remainingStates) > 0) - { - Thread.Yield(); - } - } - - if (Volatile.Read(ref remainingStates) == 0) - { - break; - } - } - } - finally - { - if (completionsCount > 0) - { - FlushPendingResults(completionsBuffer, ref completionsCount, markComposite, clearComposite, printResult); - } - - if (compositesCount > 0) - { - FlushPendingResults(compositesBuffer, ref compositesCount, markComposite, clearComposite, printResult); - } - - ArrayPool.Shared.Return(compositesBuffer, clearArray: true); - ArrayPool.Shared.Return(completionsBuffer, clearArray: true); - ArrayPool.Shared.Return(indexBuffer, clearArray: true); - ArrayPool.Shared.Return(primeBuffer, clearArray: true); - ArrayPool.Shared.Return(hitsBuffer, clearArray: true); - } - }, TaskCreationOptions.LongRunning); - } - - Task.WaitAll(workers); + Console.WriteLine("No candidates were provided for --mersenne=bydivisor."); + return; } - private static ulong AcquireNextDivisor( - ref ulong nextDivisor, - ulong divisorLimit, - ref int divisorsExhaustedFlag, - ref long finalDivisorBits, - out bool exhausted, - ref ulong localDivisorCursor, - ref int localDivisorsRemaining) + _ = threadCount; + + bool applyStartPrime = startPrime > 0UL; + int skippedByPreviousResults = 0; + List primesToTest = new(candidates.Count); + ulong maxPrime = 0UL; + + using IEnumerator primeEnumerator = Prime.Numbers.GetEnumerator(); + bool hasPrime = primeEnumerator.MoveNext(); + ulong currentPrime = hasPrime ? primeEnumerator.Current : 0UL; + + for (int index = 0; index < candidates.Count; index++) { - ref long nextDivisorBits = ref Unsafe.As(ref nextDivisor); - ulong blockStride = unchecked(DivisorAllocationBlockSize * 2UL); + ulong candidate = candidates[index]; - // TODO: Replace this odd-only linear allocator with the divisor-cycle aware stride planner once the cache exposes - // skip tables so we can jump directly to viable divisors instead of allocating every odd candidate. - if (TryAcquireLocalDivisor(ref localDivisorCursor, ref localDivisorsRemaining, out ulong localDivisor)) + if (applyStartPrime && candidate < startPrime) + { + continue; + } + + if (previousResults is not null && previousResults.ContainsKey(candidate)) + { + skippedByPreviousResults++; + continue; + } + + if (candidate <= 1UL) + { + markComposite(); + printResult(candidate, false, false, false); + continue; + } + + while (hasPrime && currentPrime < candidate) + { + hasPrime = primeEnumerator.MoveNext(); + if (!hasPrime) { - exhausted = false; - return localDivisor; + break; } - while (true) - { - if (Volatile.Read(ref divisorsExhaustedFlag) != 0) - { - exhausted = true; - return 0UL; - } + currentPrime = primeEnumerator.Current; + } - long currentBits = Volatile.Read(ref nextDivisorBits); - ulong currentValue = unchecked((ulong)currentBits); - if (currentValue > divisorLimit) - { - if (Interlocked.CompareExchange(ref divisorsExhaustedFlag, 1, 0) == 0) - { - Volatile.Write(ref finalDivisorBits, currentBits); - } + if (!hasPrime || currentPrime != candidate) + { + markComposite(); + printResult(candidate, false, false, false); + continue; + } - exhausted = true; - return 0UL; - } + if (candidate > maxPrime) + { + maxPrime = candidate; + } - ulong maximumNext = divisorLimit >= ulong.MaxValue - 1UL ? ulong.MaxValue : divisorLimit + 2UL; - ulong requestedNext = currentValue > ulong.MaxValue - blockStride ? ulong.MaxValue : currentValue + blockStride; - if (requestedNext > maximumNext) - { - requestedNext = maximumNext; - } - - long nextBits = unchecked((long)requestedNext); - if (Interlocked.CompareExchange(ref nextDivisorBits, nextBits, currentBits) != currentBits) - { - continue; - } - - ulong available = requestedNext - currentValue; - if (available == 0UL) - { - continue; - } - - int count = (int)(available >> 1); - if (count <= 0) - { - continue; - } - - localDivisorCursor = currentValue; - localDivisorsRemaining = count; - - if (TryAcquireLocalDivisor(ref localDivisorCursor, ref localDivisorsRemaining, out ulong candidate)) - { - exhausted = false; - return candidate; - } - - localDivisorCursor = 0UL; - localDivisorsRemaining = 0; - } + primesToTest.Add(candidate); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool TryAcquireLocalDivisor(ref ulong cursor, ref int remaining, out ulong divisor) + if (skippedByPreviousResults > 0) { - while (remaining > 0) - { - ulong candidate = cursor; - cursor += 2UL; - remaining--; - - // TODO: Integrate divisor-cycle skip tables here so the local allocator leaps across cached composite - // spans instead of probing every odd candidate sequentially. - if (IsAllowedDivisorCandidate(candidate)) - { - divisor = candidate; - return true; - } - } - - divisor = 0UL; - return false; + Console.WriteLine($"Skipped {skippedByPreviousResults.ToString(CultureInfo.InvariantCulture)} candidates excluded by previous results."); } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool IsAllowedDivisorCandidate(ulong divisor) + if (primesToTest.Count == 0) { - if ((divisor & 1UL) == 0UL) - { - return false; - } + if (applyStartPrime) + { + Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); + } + else + { + Console.WriteLine("No prime candidates remain for --mersenne=bydivisor after filtering."); + } - if ((divisor % 3UL) == 0UL || (divisor % 5UL) == 0UL || (divisor % 7UL) == 0UL || (divisor % 11UL) == 0UL) - { - // TODO: Replace these `%` filters with Mod3/Mod5/Mod7/Mod11 helpers so candidate sieving reuses the - // benchmarked lookup-based residues instead of executing slow modulo instructions for every divisor. - return false; - } - - return true; + return; } - private static int BuildPrimeBuffer( - ulong divisor, - in ulong[] primeValues, - in ulong[] allowedMaxValues, - int[] stateFlags, - ulong[] primeBuffer, - int[] indexBuffer, - PendingResult[] completionsBuffer, - ref int completionsCount, - ref int remainingStates, - long[] activeStateMask, - ref int activeStartIndex, - Action markComposite, - Action clearComposite, - Action printResult) + if (maxPrime <= 1UL) { - int length = primeValues.Length; - int startIndex; - while (true) - { - startIndex = Volatile.Read(ref activeStartIndex); - if (startIndex >= length || allowedMaxValues[startIndex] >= divisor) - { - break; - } - - if (Interlocked.CompareExchange(ref activeStartIndex, startIndex + 1, startIndex) != startIndex) - { - continue; - } - - int state = Volatile.Read(ref stateFlags[startIndex]); - if (state == ByDivisorStateActive) - { - ClearActiveMask(activeStateMask, startIndex); - Interlocked.Decrement(ref remainingStates); - if (completionsCount == completionsBuffer.Length) - { - FlushPendingResults(completionsBuffer, ref completionsCount, markComposite, clearComposite, printResult); - } - completionsBuffer[completionsCount++] = new PendingResult(primeValues[startIndex], detailedCheck: true, passedAllTests: true); - } - else - { - ClearActiveMask(activeStateMask, startIndex); - } - } - - int index = startIndex; - int activeCount = 0; - while (index < length) - { - int wordIndex = index >> 6; - ulong word = unchecked((ulong)activeStateMask[wordIndex]); - if (word == 0UL) - { - word = unchecked((ulong)Volatile.Read(ref activeStateMask[wordIndex])); - if (word == 0UL) - { - index = (wordIndex + 1) << 6; - continue; - } - } - - int bitOffset = index & 63; - if (bitOffset != 0) - { - word &= ulong.MaxValue << bitOffset; - if (word == 0UL) - { - index = (wordIndex + 1) << 6; - continue; - } - } - - while (word != 0UL) - { - int candidateIndex = (wordIndex << 6) + BitOperations.TrailingZeroCount(word); - if (candidateIndex >= length) - { - return activeCount; - } - - if (stateFlags[candidateIndex] == ByDivisorStateActive) - { - primeBuffer[activeCount] = primeValues[candidateIndex]; - indexBuffer[activeCount] = candidateIndex; - activeCount++; - } - else - { - ClearActiveMask(activeStateMask, candidateIndex); - } - - word &= word - 1UL; - } - - index = (wordIndex + 1) << 6; - } - - return activeCount; + Console.WriteLine("The filter specified by --filter-p must contain at least one prime exponent greater than 1 for --mersenne=bydivisor."); + return; } - private static void FinalizeRemainingStates( - ulong[] primeValues, - ulong[] allowedMaxValues, - int[] stateFlags, - ref int remainingStates, - ref long finalDivisorBits, - PendingResult[] completionsBuffer, - ref int completionsCount, - long[] activeStateMask, - Action clearComposite, - Action printResult, - Action markComposite) + tester.ConfigureFromMaxPrime(maxPrime); + + int primeCount = primesToTest.Count; + ulong[] primeBatch = ArrayPool.Shared.Rent(primeCount); + ulong[] allowedMaxBatch = ArrayPool.Shared.Rent(primeCount); + List filteredPrimes = new(primeCount); + + try { - for (int finalizeIndex = 0; finalizeIndex < primeValues.Length; finalizeIndex++) + Span primeSpan = primeBatch.AsSpan(0, primeCount); + Span allowedMaxSpan = allowedMaxBatch.AsSpan(0, primeCount); + + for (int i = 0; i < primeCount; i++) + { + primeSpan[i] = primesToTest[i]; + } + + tester.PrepareCandidates(primeSpan, allowedMaxSpan); + + for (int i = 0; i < primeCount; i++) + { + ulong prime = primeSpan[i]; + ulong allowedMax = allowedMaxSpan[i]; + + if (allowedMax < 3UL) { - if (stateFlags[finalizeIndex] != ByDivisorStateActive) - { - continue; - } - - ulong finalDivisor = unchecked((ulong)Volatile.Read(ref finalDivisorBits)); - bool detailed = finalDivisor > allowedMaxValues[finalizeIndex]; - - if (Interlocked.CompareExchange(ref stateFlags[finalizeIndex], detailed ? ByDivisorStateCompletedDetailed : ByDivisorStateCompleted, ByDivisorStateActive) != ByDivisorStateActive) - { - continue; - } - - ClearActiveMask(activeStateMask, finalizeIndex); - Interlocked.Decrement(ref remainingStates); - if (completionsCount == completionsBuffer.Length) - { - FlushPendingResults(completionsBuffer, ref completionsCount, markComposite, clearComposite, printResult); - } - - completionsBuffer[completionsCount++] = new PendingResult(primeValues[finalizeIndex], detailedCheck: detailed, passedAllTests: true); - } - } - - private static void FlushPendingResults( - PendingResult[] buffer, - ref int count, - Action markComposite, - Action clearComposite, - Action printResult) - { - for (int flushIndex = 0; flushIndex < count; flushIndex++) - { - PendingResult result = buffer[flushIndex]; - if (result.PassedAllTests) - { - clearComposite(); - } - else - { - markComposite(); - } - - printResult(result.Prime, true, result.DetailedCheck, result.PassedAllTests); + clearComposite(); + printResult(prime, true, true, true); + continue; } - count = 0; + filteredPrimes.Add(prime); + } } - - private static void ClearActiveMask(long[] activeStateMask, int index) + finally { - int wordIndex = index >> 6; - long bit = 1L << (index & 63); - - while (true) - { - long current = Volatile.Read(ref activeStateMask[wordIndex]); - if ((current & bit) == 0) - { - return; - } - - if (Interlocked.CompareExchange(ref activeStateMask[wordIndex], current & ~bit, current) == current) - { - return; - } - } + ArrayPool.Shared.Return(allowedMaxBatch, clearArray: true); + ArrayPool.Shared.Return(primeBatch, clearArray: true); } - private struct ByDivisorPrimeState + if (filteredPrimes.Count == 0) { - internal ulong Prime; - internal ulong AllowedMax; - internal bool Completed; - internal bool Composite; - internal bool DetailedCheck; + if (applyStartPrime) + { + Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); + } + + return; } - private readonly struct PendingResult + foreach (ulong prime in filteredPrimes) { - internal PendingResult(ulong prime, bool detailedCheck, bool passedAllTests) - { - Prime = prime; - DetailedCheck = detailedCheck; - PassedAllTests = passedAllTests; - } + bool isPrime = tester.IsPrime(prime, out bool divisorsExhausted); - internal ulong Prime { get; } + if (!isPrime) + { + markComposite(); + printResult(prime, true, true, false); + continue; + } - internal bool DetailedCheck { get; } - - internal bool PassedAllTests { get; } + clearComposite(); + printResult(prime, true, divisorsExhausted, true); } + } } - - -