diff --git a/AGENTS.md b/AGENTS.md index 8b1606a..0cb4f58 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -23,11 +23,31 @@ To make any changes to the files in the repository, you should always utilize a dotnet script apply_patch.cs -- [target-directory] ``` -It supports the standard `.diff` patch files as input. You should be able to run it both in the local and remote `Codex` environments. You should expect it to work both in `Windows` and `Linux` environments. +It supports the standard `git diff` patch files as input. You should be able to run it both in the local and remote `Codex` environments. You should expect it to work both in `Windows` and `Linux` environments. Make sure that you include `diff --git` lines, headers and file prefixes in the patch files to avoid applying issues. + +`apply_patch.cs` supports several parameters, depending on your needs: +* `--check` - allows checking if the patch will work without making changes to the files +* `--ignore-whitespace` +* `--ignore-eol` +* `--ignore-line-endings` +* `--target` - allows specifying target root directory for the patch + +* Build patches in standard `git diff` format. + * It supports `---`, `+++`, and `@@`. + * It doesn't support `***` or `*** Begin/End Patch` wrappers. + * It requires `diff --git` at the beginning. + * Patch header needs the actual line numbers. +* Check if the patch will work by using `--check` parameter, without making changes to the files, which would force patch changes. +* Build patches in small, line-precise hunks; confirm target lines with `Get-Content` or `nl -ba` before writing the diff to avoid context mismatches. +* Regenerate the diff whenever the file shifts—never reuse an old hunk after other edits. +* Keep patch filenames unique and delete applied files to prevent accidental re-use. +* When a hunk fails, inspect the around-lines immediately and adjust rather than retrying the same diff. +* Prefer single-purpose patches (one logical change per file) so rollbacks or fixes stay focused. +* Use `nl -ba` to show line numbers, when needed. If you want to make any changes, create a patch file and run the script to apply it. Always run it from the root directory and reference any modified files with relative paths respectively in the patch. -Do **not** attempt to apply the changes manually, or with `Python`, or with `PowerShell` scripts, or `cat` or any other way. If you identify any issues with `apply_patch.cs` script, you should propose fixes to resolve them and continue using the script after updates. If you spot any missing, but required features, propose enhancements and continue using the script after updates, too. +Do **not** attempt to apply the changes manually, or with `Python`, or with `PowerShell` scripts, or `cat` or any other way, unless you apply changes to `apply_patch.cs` script itself. If you identify any issues with `apply_patch.cs` script, you should propose fixes to resolve them and continue using the script after updates. If you spot any missing, but required features, propose enhancements and continue using the script after updates, too. ### Test execution time policy 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); } + } } - - -