diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index d269982..d64ef23 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,6 +1,6 @@ using System; +using System.Buffers; using System.Collections.Concurrent; -using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; @@ -200,12 +200,25 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return false; } - Dictionary? factorCache = null; + int batchCapacity = Math.Max(1, _batchSize); DivisorCycleCache cycleCache = DivisorCycleCache.Shared; + int preferredCycleBatch = cycleCache.PreferredBatchSize; + if (preferredCycleBatch > 0 && preferredCycleBatch < batchCapacity) + { + batchCapacity = preferredCycleBatch; + } - ulong cachedModulus = 0UL; - MontgomeryDivisorData cachedDivisorData = default; - bool hasCachedDivisorData = false; + ulong[] candidateBuffer = ArrayPool.Shared.Rent(batchCapacity); + ulong[] cycleBuffer = ArrayPool.Shared.Rent(batchCapacity); + MontgomeryDivisorData[] divisorDataBuffer = ArrayPool.Shared.Rent(batchCapacity); + byte[] remainderBuffer = ArrayPool.Shared.Rent(batchCapacity * 6); + + Span remainder10Span = remainderBuffer.AsSpan(0, batchCapacity); + Span remainder8Span = remainderBuffer.AsSpan(batchCapacity, batchCapacity); + Span remainder5Span = remainderBuffer.AsSpan(batchCapacity * 2, batchCapacity); + Span remainder3Span = remainderBuffer.AsSpan(batchCapacity * 3, batchCapacity); + Span remainder7Span = remainderBuffer.AsSpan(batchCapacity * 4, batchCapacity); + Span remainder11Span = remainderBuffer.AsSpan(batchCapacity * 5, batchCapacity); byte step10 = (byte)(step % 10UL); byte step8 = (byte)(step % 8UL); @@ -221,83 +234,151 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv byte remainder7 = (byte)(divisor % 7UL); byte remainder11 = (byte)(divisor % 11UL); - // Keep the divisibility filters aligned with the divisor-cycle generator so the - // CPU path never requests cycles that were skipped during cache creation. bool lastIsSeven = (prime & 3UL) == 3UL; - while (divisor <= limit) + try { - if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) + while (divisor <= limit) { - processedAll = false; - return false; - } + Span candidateSpan = candidateBuffer.AsSpan(0, batchCapacity); + int chunkCount = 0; + UInt128 localDivisor = divisor; - iterationCounter++; - ulong candidate = (ulong)divisor; - processedCount++; - lastProcessed = candidate; + byte localRemainder10 = remainder10; + byte localRemainder8 = remainder8; + byte localRemainder5 = remainder5; + byte localRemainder3 = remainder3; + byte localRemainder7 = remainder7; + byte localRemainder11 = remainder11; - bool admissible = lastIsSeven - ? (remainder10 == 3 || remainder10 == 7 || remainder10 == 9) - : (remainder10 == 1 || remainder10 == 3 || remainder10 == 9); - - if (admissible && (remainder8 == 1 || remainder8 == 7) && remainder3 != 0 && remainder5 != 0 && remainder7 != 0 && remainder11 != 0) - { - if (!hasCachedDivisorData || cachedModulus != candidate) + while (chunkCount < batchCapacity && localDivisor <= limit) { - cachedDivisorData = MontgomeryDivisorData.FromModulus(candidate); - cachedModulus = candidate; - hasCachedDivisorData = true; - } - - MontgomeryDivisorData divisorData = cachedDivisorData; - ulong divisorCycle; - - if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) - { - divisorCycle = cycleCache.GetCycleLength(candidate); - } - else - { - factorCache ??= new Dictionary(8); - // EvenPerfectBitScanner reaches this path only for divisors above the precomputed cycle snapshot. - // Those candidates follow q = 2 * prime * k + 1 with prime >= 5, so the downstream order calculator - // can rely on an odd divisor and a prime exponent when deriving the cycle length. - if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) + if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) { - divisorCycle = cycleCache.GetCycleLength(candidate); + processedAll = false; + return false; + } + + ulong candidate = (ulong)localDivisor; + candidateSpan[chunkCount] = candidate; + remainder10Span[chunkCount] = localRemainder10; + remainder8Span[chunkCount] = localRemainder8; + remainder5Span[chunkCount] = localRemainder5; + remainder3Span[chunkCount] = localRemainder3; + remainder7Span[chunkCount] = localRemainder7; + remainder11Span[chunkCount] = localRemainder11; + + processedCount++; + lastProcessed = candidate; + iterationCounter++; + chunkCount++; + + localDivisor += step; + localRemainder10 = AddMod(localRemainder10, step10, (byte)10); + localRemainder8 = AddMod(localRemainder8, step8, (byte)8); + localRemainder5 = AddMod(localRemainder5, step5, (byte)5); + localRemainder3 = AddMod(localRemainder3, step3, (byte)3); + localRemainder7 = AddMod(localRemainder7, step7, (byte)7); + localRemainder11 = AddMod(localRemainder11, step11, (byte)11); + } + + if (chunkCount == 0) + { + break; + } + + divisor = localDivisor; + remainder10 = localRemainder10; + remainder8 = localRemainder8; + remainder5 = localRemainder5; + remainder3 = localRemainder3; + remainder7 = localRemainder7; + remainder11 = localRemainder11; + + // Reuse the candidate buffer in-place so the filtered results stay contiguous without additional arrays. + Span chunkCandidates = candidateSpan[..chunkCount]; + int filteredCount = 0; + + for (int i = 0; i < chunkCount; i++) + { + bool admissible = lastIsSeven + ? (remainder10Span[i] == 3 || remainder10Span[i] == 7 || remainder10Span[i] == 9) + : (remainder10Span[i] == 1 || remainder10Span[i] == 3 || remainder10Span[i] == 9); + + if (!admissible) + { + continue; + } + + if ((remainder8Span[i] != 1 && remainder8Span[i] != 7) + || remainder3Span[i] == 0 + || remainder5Span[i] == 0 + || remainder7Span[i] == 0 + || remainder11Span[i] == 0) + { + continue; + } + + // Store the accepted candidate at the next free slot so the front of the span becomes the filtered set. + chunkCandidates[filteredCount++] = chunkCandidates[i]; + } + + if (filteredCount == 0) + { + continue; + } + + Span filteredSpan = chunkCandidates[..filteredCount]; + Span cycleSpan = cycleBuffer.AsSpan(0, filteredCount); + cycleCache.GetCycleLengths(filteredSpan, cycleSpan); + + int admissibleCount = 0; + Span divisorDataSpan = divisorDataBuffer.AsSpan(0, filteredCount); + + for (int i = 0; i < filteredCount; i++) + { + if (cycleSpan[i] != prime) + { + continue; + } + + ulong candidate = filteredSpan[i]; + // Reuse the filtered span's leading segment to accumulate the admissible divisors for Montgomery checks. + filteredSpan[admissibleCount] = candidate; + divisorDataSpan[admissibleCount] = MontgomeryDivisorData.FromModulus(candidate); + admissibleCount++; + } + + if (admissibleCount == 0) + { + continue; + } + + Span admissibleSpan = filteredSpan[..admissibleCount]; + Span admissibleDataSpan = divisorDataSpan[..admissibleCount]; + + for (int i = 0; i < admissibleCount; i++) + { + if (CheckDivisor(prime, prime, in admissibleDataSpan[i]) != 0) + { + lastProcessed = admissibleSpan[i]; + processedAll = true; + return true; } } - - if (divisorCycle == prime) - { - // A cycle equal to the tested exponent (which is prime in this path) - // guarantees that the candidate divides the corresponding Mersenne - // number because the order of 2 modulo the divisor is exactly p. - processedAll = true; - return true; - } - - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } } - - divisor += step; - remainder10 = AddMod(remainder10, step10, (byte)10); - remainder8 = AddMod(remainder8, step8, (byte)8); - remainder5 = AddMod(remainder5, step5, (byte)5); - remainder3 = AddMod(remainder3, step3, (byte)3); - remainder7 = AddMod(remainder7, step7, (byte)7); - remainder11 = AddMod(remainder11, step11, (byte)11); + } + finally + { + ArrayPool.Shared.Return(candidateBuffer, clearArray: false); + ArrayPool.Shared.Return(cycleBuffer, clearArray: false); + ArrayPool.Shared.Return(divisorDataBuffer, clearArray: false); + ArrayPool.Shared.Return(remainderBuffer, clearArray: false); } processedAll = divisor > limit; return false; } - private static byte CheckDivisor(ulong prime, ulong divisorCycle, in MontgomeryDivisorData divisorData) { ulong residue = prime.Pow2MontgomeryModWithCycleCpu(divisorCycle, divisorData); @@ -372,6 +453,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private bool _hasCachedDivisorData; private ulong _cachedDivisor; private MontgomeryDivisorData _cachedDivisorData; + private ulong[]? _primeDeltas; + private ulong[]? _cycleRemainders; + private ulong[]? _montgomeryDeltas; + private int _bufferCapacity; internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) { @@ -412,6 +497,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv cachedData = _cachedDivisorData; } + ulong modulus = cachedData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits.Clear(); + return; + } + if (divisorCycle == 0UL) { divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); @@ -422,35 +514,167 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } } - // Keep these remainder steppers in place so future updates continue reusing the previously computed residues. - // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors. - var exponentStepper = new ExponentRemainderStepper(cachedData); - if (!exponentStepper.IsValidModulus) - { - hits.Clear(); - return; - } + EnsureCapacity(length); - var cycleStepper = new CycleRemainderStepper(divisorCycle); + Span deltaSpan = _primeDeltas!.AsSpan(0, length); + ComputePrimeDeltas(primes, deltaSpan); - ulong remainder = cycleStepper.Initialize(primes[0]); - hits[0] = remainder == 0UL - ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0) - : (byte)0; + Span remainderSpan = _cycleRemainders!.AsSpan(0, length); + ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); + + Span montgomerySpan = _montgomeryDeltas!.AsSpan(0, length); + ComputeMontgomeryDeltas(cachedData, primes, deltaSpan, montgomerySpan); + + ulong nPrime = cachedData.NPrime; + ulong montgomeryOne = cachedData.MontgomeryOne; + ulong currentMontgomery = montgomerySpan[0]; + + hits[0] = remainderSpan[0] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0; for (int i = 1; i < length; i++) { - remainder = cycleStepper.ComputeNext(primes[i]); - if (remainder != 0UL) + currentMontgomery = currentMontgomery.MontgomeryMultiply(montgomerySpan[i], modulus, nPrime); + hits[i] = remainderSpan[i] == 0UL && currentMontgomery == montgomeryOne ? (byte)1 : (byte)0; + } + } + + private void EnsureCapacity(int requiredLength) + { + if (requiredLength <= 0) + { + return; + } + + if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _montgomeryDeltas is not null) + { + return; + } + + int newCapacity = _bufferCapacity; + if (newCapacity == 0) + { + newCapacity = 1; + } + + while (newCapacity < requiredLength) + { + newCapacity <<= 1; + } + + ReturnBuffers(); + + _primeDeltas = ArrayPool.Shared.Rent(newCapacity); + _cycleRemainders = ArrayPool.Shared.Rent(newCapacity); + _montgomeryDeltas = ArrayPool.Shared.Rent(newCapacity); + _bufferCapacity = newCapacity; + } + + private void ReturnBuffers() + { + if (_primeDeltas is not null) + { + ArrayPool.Shared.Return(_primeDeltas, clearArray: false); + _primeDeltas = null; + } + + if (_cycleRemainders is not null) + { + ArrayPool.Shared.Return(_cycleRemainders, clearArray: false); + _cycleRemainders = null; + } + + if (_montgomeryDeltas is not null) + { + ArrayPool.Shared.Return(_montgomeryDeltas, clearArray: false); + _montgomeryDeltas = null; + } + + _bufferCapacity = 0; + } + + private static void ComputePrimeDeltas(ReadOnlySpan primes, Span deltas) + { + if (primes.Length == 0) + { + return; + } + + ulong previous = primes[0]; + deltas[0] = previous; + + for (int i = 1; i < primes.Length; i++) + { + ulong current = primes[i]; + if (current <= previous) { - hits[i] = 0; - continue; + throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing."); } - hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; + deltas[i] = current - previous; + previous = current; } } + private static void ComputeCycleRemainders(ulong cycleLength, ReadOnlySpan primes, ReadOnlySpan deltas, Span remainders) + { + if (cycleLength == 0UL || primes.Length == 0) + { + remainders.Clear(); + return; + } + + ulong remainder = primes[0] % cycleLength; + remainders[0] = remainder; + + for (int i = 1; i < primes.Length; i++) + { + remainder = AddMod(remainder, deltas[i], cycleLength); + remainders[i] = remainder; + } + } + + private static void ComputeMontgomeryDeltas(in MontgomeryDivisorData divisorData, ReadOnlySpan primes, ReadOnlySpan deltas, Span montgomeryValues) + { + if (primes.Length == 0) + { + montgomeryValues.Clear(); + return; + } + + montgomeryValues[0] = primes[0].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true); + + for (int i = 1; i < primes.Length; i++) + { + montgomeryValues[i] = deltas[i].Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: true); + } + } + + private static ulong AddMod(ulong value, ulong delta, ulong modulus) + { + if (modulus == 0UL) + { + return 0UL; + } + + if (ulong.MaxValue - value < delta) + { + UInt128 extended = (UInt128)value + delta; + return (ulong)(extended % modulus); + } + + ulong sum = value + delta; + if (sum >= modulus) + { + sum -= modulus; + if (sum >= modulus) + { + sum %= modulus; + } + } + + return sum; + } + public void Dispose() { if (_disposed)