using System; using System.Buffers; using System.Collections.Concurrent; using System.Collections.Generic; using System.Numerics; using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Cpu; public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDivisorByDivisorTester { private readonly object _sync = new(); private readonly ConcurrentBag _sessionPool = new(); private ulong _divisorLimit; private ulong _lastStatusDivisor; private bool _isConfigured; private int _batchSize = 1_024; private ulong _configuredMaxPrime; private ByDivisorDeltasDevice _deltasDevice = ByDivisorDeltasDevice.Cpu; public int BatchSize { get => _batchSize; set { int sanitized = Math.Max(1, value); _batchSize = sanitized; } } public ByDivisorDeltasDevice DeltasDevice { get => _deltasDevice; set { lock (_sync) { _deltasDevice = value; } } } public void ConfigureFromMaxPrime(ulong maxPrime) { lock (_sync) { _configuredMaxPrime = maxPrime; _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); _lastStatusDivisor = 0UL; _isConfigured = true; } } public ulong DivisorLimit { get { lock (_sync) { if (!_isConfigured) { throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } return _divisorLimit; } } } public ulong GetAllowedMaxDivisor(ulong prime) { lock (_sync) { if (!_isConfigured) { throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } return ComputeAllowedMaxDivisor(prime, _divisorLimit); } } public bool IsPrime(ulong prime, out bool divisorsExhausted, TimeSpan? timeLimit = null) { ulong allowedMax; lock (_sync) { if (!_isConfigured) { throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); } if (allowedMax < 3UL) { divisorsExhausted = true; return true; } ulong processedCount; ulong lastProcessed; bool processedAll; bool composite = CheckDivisors( prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); if (processedCount > 0UL) { lock (_sync) { UpdateStatusUnsafe(processedCount); } } if (composite) { divisorsExhausted = true; return false; } divisorsExhausted = processedAll || composite; return true; } public void PrepareCandidates(ReadOnlySpan primes, Span allowedMaxValues) { if (allowedMaxValues.Length < primes.Length) { throw new ArgumentException("allowedMaxValues span must be at least as long as primes span.", nameof(allowedMaxValues)); } ulong divisorLimit; lock (_sync) { if (!_isConfigured) { throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } divisorLimit = _divisorLimit; } for (int index = 0; index < primes.Length; index++) { allowedMaxValues[index] = ComputeAllowedMaxDivisor(primes[index], divisorLimit); } } public IMersenneNumberDivisorByDivisorTester.IDivisorScanSession CreateDivisorSession() { lock (_sync) { if (!_isConfigured) { throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } if (_sessionPool.TryTake(out DivisorScanSession? session)) { session.Reset(); return session; } return new DivisorScanSession(this); } } private void ReturnSession(DivisorScanSession session) { _sessionPool.Add(session); } private bool CheckDivisors( ulong prime, ulong allowedMax, out ulong lastProcessed, out bool processedAll, out ulong processedCount, TimeSpan? timeLimit) { if (_deltasDevice == ByDivisorDeltasDevice.Cpu) { return CheckDivisorsCpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); } return CheckDivisorsGpu(prime, allowedMax, out lastProcessed, out processedAll, out processedCount, timeLimit); } private bool CheckDivisorsCpu( ulong prime, ulong allowedMax, out ulong lastProcessed, out bool processedAll, out ulong processedCount, TimeSpan? timeLimit) { lastProcessed = 0UL; processedCount = 0UL; processedAll = false; PrimeTestTimeLimit limitGuard; if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard)) { return false; } bool enforceLimit = limitGuard.IsActive; ulong iterationCounter = 0UL; if (allowedMax < 3UL) { return false; } UInt128 step = (UInt128)prime << 1; if (step == UInt128.Zero) { processedAll = true; return false; } UInt128 limit = allowedMax; UInt128 divisor = step + UInt128.One; if (divisor > limit) { processedAll = true; return false; } Dictionary? factorCache = null; DivisorCycleCache cycleCache = DivisorCycleCache.Shared; ulong cachedModulus = 0UL; MontgomeryDivisorData cachedDivisorData = default; bool hasCachedDivisorData = false; byte step10 = (byte)(step % 10UL); byte step8 = (byte)(step % 8UL); byte step5 = (byte)(step % 5UL); byte step3 = (byte)(step % 3UL); byte step7 = (byte)(step % 7UL); byte step11 = (byte)(step % 11UL); byte remainder10 = (byte)(divisor % 10UL); byte remainder8 = (byte)(divisor % 8UL); byte remainder5 = (byte)(divisor % 5UL); byte remainder3 = (byte)(divisor % 3UL); 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) { if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) { processedAll = false; return false; } iterationCounter++; ulong candidate = (ulong)divisor; processedCount++; lastProcessed = candidate; 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) { 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) { divisorCycle = cycleCache.GetCycleLength(candidate); } } 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); } processedAll = divisor > limit; return false; } private bool CheckDivisorsGpu( ulong prime, ulong allowedMax, out ulong lastProcessed, out bool processedAll, out ulong processedCount, TimeSpan? timeLimit) { lastProcessed = 0UL; processedCount = 0UL; processedAll = false; PrimeTestTimeLimit limitGuard; if (!PrimeTestTimeLimit.TryCreate(timeLimit, out limitGuard)) { return false; } bool enforceLimit = limitGuard.IsActive; ulong iterationCounter = 0UL; if (allowedMax < 3UL) { return false; } UInt128 step = (UInt128)prime << 1; if (step == UInt128.Zero) { processedAll = true; return false; } UInt128 limit = allowedMax; UInt128 divisor = step + UInt128.One; if (divisor > limit) { processedAll = true; return false; } int batchCapacity = Math.Max(1, _batchSize); DivisorCycleCache cycleCache = DivisorCycleCache.Shared; int preferredCycleBatch = cycleCache.PreferredBatchSize; if (preferredCycleBatch > 0 && preferredCycleBatch < batchCapacity) { batchCapacity = preferredCycleBatch; } 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); byte step5 = (byte)(step % 5UL); byte step3 = (byte)(step % 3UL); byte step7 = (byte)(step % 7UL); byte step11 = (byte)(step % 11UL); byte remainder10 = (byte)(divisor % 10UL); byte remainder8 = (byte)(divisor % 8UL); byte remainder5 = (byte)(divisor % 5UL); byte remainder3 = (byte)(divisor % 3UL); byte remainder7 = (byte)(divisor % 7UL); byte remainder11 = (byte)(divisor % 11UL); bool lastIsSeven = (prime & 3UL) == 3UL; try { while (divisor <= limit) { Span candidateSpan = candidateBuffer.AsSpan(0, batchCapacity); int chunkCount = 0; UInt128 localDivisor = divisor; byte localRemainder10 = remainder10; byte localRemainder8 = remainder8; byte localRemainder5 = remainder5; byte localRemainder3 = remainder3; byte localRemainder7 = remainder7; byte localRemainder11 = remainder11; while (chunkCount < batchCapacity && localDivisor <= limit) { if (enforceLimit && (iterationCounter == 0UL || (iterationCounter & 63UL) == 0UL) && limitGuard.HasExpired()) { 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; } } } } 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); return residue == 1UL ? (byte)1 : (byte)0; } private static byte AddMod(byte value, byte delta, byte modulus) { int sum = value + delta; if (sum >= modulus) { sum -= modulus; } return (byte)sum; } private void UpdateStatusUnsafe(ulong processedCount) { if (processedCount == 0UL) { return; } ulong interval = PerfectNumberConstants.ConsoleInterval; if (interval == 0UL) { _lastStatusDivisor = 0UL; return; } ulong total = _lastStatusDivisor + processedCount; // TODO: Replace this modulo with the ring-buffer style counter (subtract loop) used in the fast CLI // status benchmarks so we avoid `%` in this hot loop while still wrapping progress correctly. _lastStatusDivisor = total % interval; } private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) { if (maxPrime <= 1UL) { return 0UL; } if (maxPrime - 1UL >= 64UL) { return ulong.MaxValue; } return (1UL << (int)(maxPrime - 1UL)) - 1UL; } private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { if (prime <= 1UL) { return 0UL; } if (prime - 1UL >= 64UL) { return divisorLimit; } return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); } private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession { private readonly MersenneNumberDivisorByDivisorCpuTester _owner; private bool _disposed; private bool _hasCachedDivisorData; private ulong _cachedDivisor; private MontgomeryDivisorData _cachedDivisorData; private ulong[]? _primeDeltas; private ulong[]? _cycleRemainders; private ulong[]? _montgomeryDeltas; private int _bufferCapacity; private MersenneNumberDivisorResidueGpuEvaluator? _gpuResidueEvaluator; internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) { _owner = owner; } internal void Reset() { _disposed = false; _hasCachedDivisorData = false; _cachedDivisor = 0UL; _cachedDivisorData = default; } public void CheckDivisor(ulong divisor, in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits) { if (_disposed) { throw new ObjectDisposedException(nameof(DivisorScanSession)); } int length = primes.Length; if (length == 0) { return; } MontgomeryDivisorData cachedData = divisorData; if (cachedData.Modulus != divisor) { if (!_hasCachedDivisorData || _cachedDivisor != divisor) { _cachedDivisorData = MontgomeryDivisorData.FromModulus(divisor); _cachedDivisor = divisor; _hasCachedDivisorData = true; } cachedData = _cachedDivisorData; } ulong modulus = cachedData.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) { hits.Clear(); return; } if (divisorCycle == 0UL) { divisorCycle = DivisorCycleCache.Shared.GetCycleLength(divisor); if (divisorCycle == 0UL) { hits.Clear(); return; } } if (_owner._deltasDevice == ByDivisorDeltasDevice.Cpu) { CheckDivisorCpu(cachedData, divisorCycle, primes, hits); return; } CheckDivisorGpu(cachedData, divisorCycle, primes, hits, length); } private static void CheckDivisorCpu( in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits) { // 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(divisorData); if (!exponentStepper.IsValidModulus) { hits.Clear(); return; } var cycleStepper = new CycleRemainderStepper(divisorCycle); ulong remainder = cycleStepper.Initialize(primes[0]); hits[0] = remainder == 0UL ? (exponentStepper.ComputeNextIsUnity(primes[0]) ? (byte)1 : (byte)0) : (byte)0; int length = primes.Length; for (int i = 1; i < length; i++) { remainder = cycleStepper.ComputeNext(primes[i]); if (remainder != 0UL) { hits[i] = 0; continue; } hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; } } private void CheckDivisorGpu( in MontgomeryDivisorData divisorData, ulong divisorCycle, ReadOnlySpan primes, Span hits, int length) { EnsureCapacity(length); Span deltaSpan = _primeDeltas!.AsSpan(0, length); ComputePrimeDeltas(primes, deltaSpan); Span remainderSpan = _cycleRemainders!.AsSpan(0, length); ComputeCycleRemainders(divisorCycle, primes, deltaSpan, remainderSpan); EnsureGpuEvaluator(); Span residueSpan = _montgomeryDeltas!.AsSpan(0, length); _gpuResidueEvaluator!.ComputeResidues(primes, divisorData, residueSpan); ulong montgomeryOne = divisorData.MontgomeryOne; for (int i = 0; i < length; i++) { hits[i] = remainderSpan[i] == 0UL && residueSpan[i] == montgomeryOne ? (byte)1 : (byte)0; } } private void EnsureGpuEvaluator() { if (_gpuResidueEvaluator is null) { _gpuResidueEvaluator = new MersenneNumberDivisorResidueGpuEvaluator(_owner._batchSize); } } private void ReleaseGpuEvaluator() { if (_gpuResidueEvaluator is null) { return; } _gpuResidueEvaluator.Dispose(); _gpuResidueEvaluator = null; } 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) { throw new ArgumentOutOfRangeException(nameof(primes), "Primes must be strictly increasing."); } 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 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) { return; } _disposed = true; ReleaseGpuEvaluator(); _owner.ReturnSession(this); } } }