diff --git a/AGENTS.md b/AGENTS.md index 5125129..a32973e 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -56,7 +56,9 @@ Never ever run tests marked with trait Category = Slow. Only run tests with cate There is only one agents file - this one. Don't search for another one. -Always use indentation with 4 spaces and / or 1 tab in source code files. +Always use indentation with 4 spaces in source code files. Do not use 8-space indentation or tabs. + +Keep the divisor cycle cache limited to three blocks in memory (the base snapshot plus at most two additional blocks) and always start background computation of the next block immediately after loading the first block from disk and whenever work begins on the most recently generated block. All code, comments, commit messages, branch names and PR descriptions must be written in American English. diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index fcbd93c..6abf920 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -457,7 +457,8 @@ internal static class Program MersenneDivisorCycles.Shared.LoadFrom(cyclesPath); } - DivisorCycleCache.Shared.ReloadFromCurrentSnapshot(); + DivisorCycleCache.Shared.ConfigureGeneratorDevice(orderOnGpu); + DivisorCycleCache.Shared.ReloadFromCurrentSnapshot(); Console.WriteLine("Divisor cycles are ready"); diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 1965ad1..92df4ca 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,5 +1,4 @@ using System; -using System.Buffers; using System.Collections.Concurrent; using System.Runtime.CompilerServices; using PerfectNumbers.Core; @@ -8,556 +7,368 @@ using MontgomeryDivisorData = PerfectNumbers.Core.Gpu.MersenneNumberDivisorByDiv namespace PerfectNumbers.Core.Cpu; -// TODO: The current implementation is over-complicated for CPU. We shouldn't need all the hassle with intermediate buffers, checking the status there etc. 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 bool _useDivisorCycles; - private int _batchSize = 1_024; - - public bool UseDivisorCycles - { - get => _useDivisorCycles; - set => _useDivisorCycles = value; - } - - public int BatchSize - { - get => _batchSize; - set => _batchSize = Math.Max(1, value); - } - - public void ConfigureFromMaxPrime(ulong maxPrime) - { - lock (_sync) - { - _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) - { - ulong allowedMax; - bool useCycles; - int batchCapacity; - - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); - useCycles = _useDivisorCycles; - batchCapacity = _batchSize; - } - - if (allowedMax < 3UL) - { - divisorsExhausted = true; - return true; - } - - ulong[] divisors = ArrayPool.Shared.Rent(batchCapacity); - byte[] hits = ArrayPool.Shared.Rent(batchCapacity); - MontgomeryDivisorData[] divisorData = ArrayPool.Shared.Rent(batchCapacity); - - ulong processedCount = 0UL; - ulong lastProcessed = 0UL; - bool composite = false; - bool processedAll = false; - - try - { - composite = CheckDivisors( - prime, - allowedMax, - useCycles, - divisors, - hits, - divisorData, - out lastProcessed, - out processedAll, - out processedCount); - } - finally - { - ArrayPool.Shared.Return(divisors, clearArray: false); - ArrayPool.Shared.Return(hits, clearArray: false); - ArrayPool.Shared.Return(divisorData, clearArray: false); - } - - if (processedCount > 0UL) - { - lock (_sync) - { - UpdateStatusUnsafe(lastProcessed, processedCount); - } - } - - if (composite) - { - divisorsExhausted = true; - return false; - } - - divisorsExhausted = processedAll || composite; - return true; - } - - 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, - bool useCycles, - ulong[] divisors, - byte[] hits, - MontgomeryDivisorData[] divisorData, - out ulong lastProcessed, - out bool processedAll, - out ulong processedCount) - { - lastProcessed = 0UL; - processedCount = 0UL; - processedAll = false; - - if (allowedMax < 3UL) - { - return false; - } - - int batchCapacity = _batchSize; - ulong divisor = 3UL; - bool composite = false; - - while (divisor <= allowedMax) - { - int batchSize = 0; - bool reachedEndInBatch = false; - - while (batchSize < batchCapacity && divisor <= allowedMax) - { - ulong currentDivisor = divisor; - ulong nextDivisor = currentDivisor + 2UL; - - divisors[batchSize++] = currentDivisor; - processedCount++; - lastProcessed = currentDivisor; - - if (nextDivisor <= currentDivisor) - { - reachedEndInBatch = true; - break; - } - - if (nextDivisor > allowedMax) - { - divisor = nextDivisor; - reachedEndInBatch = true; - break; - } - - divisor = nextDivisor; - } - - if (batchSize == 0) - { - if (reachedEndInBatch) - { - processedAll = true; - } - - break; - } - - Span divisorDataSpan = divisorData.AsSpan(0, batchSize); - Span divisorSpan = divisors.AsSpan(0, batchSize); - Span hitsSpan = hits.AsSpan(0, batchSize); - - for (int i = 0; i < batchSize; i++) - { - divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]); - } - - for (int i = 0; i < batchSize; i++) - { - hitsSpan[i] = CheckDivisor(prime, useCycles, divisorSpan[i], divisorDataSpan[i]); - if (hitsSpan[i] != 0) - { - composite = true; - lastProcessed = divisorSpan[i]; - break; - } - } - - if (!composite) - { - lastProcessed = divisorSpan[batchSize - 1]; - } - - if (composite) - { - break; - } - - if (reachedEndInBatch) - { - processedAll = true; - break; - } - } - - processedAll = processedAll || divisor > allowedMax; - return composite; - } - - private static byte CheckDivisor(ulong prime, bool useCycles, ulong divisor, in MontgomeryDivisorData divisorData) - { - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return 0; - } - - ulong exponent = prime; - - if (useCycles) - { - ulong cycle = MersenneDivisorCycles.CalculateCycleLength(divisor); - if (cycle == 0UL) - { - return 0; - } - - ulong remainder = prime % cycle; - if (remainder != 0UL) - { - return 0; - } - - exponent = remainder; - } - - return exponent.Pow2MontgomeryMod(divisorData) == 1UL ? (byte)1 : (byte)0; - } - - private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount) - { - if (processedCount == 0UL) - { - return; - } - - ulong interval = PerfectNumberConstants.ConsoleInterval; - if (interval == 0UL) - { - _lastStatusDivisor = 0UL; - return; - } - - ulong total = _lastStatusDivisor + processedCount; - _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 static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) - { - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); - } - - return new MontgomeryDivisorData( - modulus, - ComputeMontgomeryNPrime(modulus), - ComputeMontgomeryResidue(1UL, modulus), - ComputeMontgomeryResidue(2UL, modulus)); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); - - private static ulong ComputeMontgomeryNPrime(ulong modulus) - { - ulong inv = modulus; - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - return unchecked(0UL - inv); - } - - private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession - { - private readonly MersenneNumberDivisorByDivisorCpuTester _owner; - private ulong[] _primeBuffer; - private int[] _positionBuffer; - private ulong[] _exponentBuffer; - private int _capacity; - private bool _disposed; - - internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) - { - _owner = owner; - _capacity = Math.Max(1, owner._batchSize); - _primeBuffer = ArrayPool.Shared.Rent(_capacity); - _positionBuffer = ArrayPool.Shared.Rent(_capacity); - _exponentBuffer = ArrayPool.Shared.Rent(_capacity); - } - - internal void Reset() - { - _disposed = false; - } - - public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) - { - if (_disposed) - { - throw new ObjectDisposedException(nameof(DivisorScanSession)); - } - - int primesLength = primes.Length; - if (primesLength == 0) - { - return; - } - - MersenneNumberDivisorByDivisorCpuTester owner = _owner; - int batchSize = owner._batchSize; - EnsureCapacity(batchSize); - - ulong modulus = divisor; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return; - } - - bool cycleEnabled = owner._useDivisorCycles && divisorCycle != 0UL; - ulong cycle = divisorCycle; - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); - - int offset = 0; - Span primeSpan = _primeBuffer.AsSpan(); - Span positionSpan = _positionBuffer.AsSpan(); - Span exponentSpan = _exponentBuffer.AsSpan(); - ulong previousResidue = 0UL; - ulong previousPrime = 0UL; - bool hasPreviousResidue = false; - - while (offset < primesLength) - { - int sliceLength = Math.Min(batchSize, primesLength - offset); - ReadOnlySpan primesSlice = primes.Slice(offset, sliceLength); - Span hitsSlice = hits.Slice(offset, sliceLength); - hitsSlice.Clear(); - - int computeCount = 0; - - if (cycleEnabled) - { - for (int i = 0; i < sliceLength; i++) - { - ulong prime = primesSlice[i]; - if (prime % cycle != 0UL) - { - continue; - } - - primeSpan[computeCount] = prime; - positionSpan[computeCount] = i; - computeCount++; - } - } - else - { - for (int i = 0; i < sliceLength; i++) - { - primeSpan[computeCount] = primesSlice[i]; - positionSpan[computeCount] = i; - computeCount++; - } - } - - if (computeCount > 0) - { - Span exponentSlice = exponentSpan[..computeCount]; - ulong currentPrime = hasPreviousResidue ? previousPrime : 0UL; - bool hasDeltaSource = hasPreviousResidue; - - for (int i = 0; i < computeCount; i++) - { - ulong prime = primeSpan[i]; - ulong exponentValue; - if (hasDeltaSource) - { - exponentValue = prime - currentPrime; - } - else - { - exponentValue = prime; - hasDeltaSource = true; - } - - if (cycleEnabled) - { - exponentValue %= cycle; - } - - exponentSlice[i] = exponentValue; - currentPrime = prime; - } - - for (int i = 0; i < computeCount; i++) - { - exponentSlice[i] = exponentSlice[i].Pow2MontgomeryMod(divisorData); - } - - for (int i = 0; i < computeCount; i++) - { - int position = positionSpan[i]; - ulong stepResidue = exponentSlice[i]; - ulong prime = primeSpan[i]; - ulong residue = hasPreviousResidue ? MultiplyMod(previousResidue, stepResidue, modulus) : stepResidue; - hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0; - previousResidue = residue; - previousPrime = prime; - hasPreviousResidue = true; - } - } - - offset += sliceLength; - } - } - - public void Dispose() - { - if (_disposed) - { - return; - } - - _disposed = true; - _owner.ReturnSession(this); - } - - private void EnsureCapacity(int requiredCapacity) - { - if (requiredCapacity <= _capacity) - { - return; - } - - ArrayPool.Shared.Return(_primeBuffer, clearArray: false); - ArrayPool.Shared.Return(_positionBuffer, clearArray: false); - ArrayPool.Shared.Return(_exponentBuffer, clearArray: false); - - _capacity = requiredCapacity; - _primeBuffer = ArrayPool.Shared.Rent(_capacity); - _positionBuffer = ArrayPool.Shared.Rent(_capacity); - _exponentBuffer = ArrayPool.Shared.Rent(_capacity); - } - - private static ulong MultiplyMod(ulong left, ulong right, ulong modulus) - { - if (modulus == 0UL) - { - return 0UL; - } - - UInt128 product = (UInt128)left * right; - return (ulong)(product % modulus); - } - } + private readonly object _sync = new(); + private readonly ConcurrentBag _sessionPool = new(); + private ulong _divisorLimit; + private ulong _lastStatusDivisor; + private bool _isConfigured; + private bool _useDivisorCycles; + private int _batchSize = 1_024; + + public bool UseDivisorCycles + { + get => _useDivisorCycles; + set => _useDivisorCycles = value; + } + + public int BatchSize + { + get => _batchSize; + set => _batchSize = Math.Max(1, value); + } + + public void ConfigureFromMaxPrime(ulong maxPrime) + { + lock (_sync) + { + _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) + { + ulong allowedMax; + bool useCycles; + + lock (_sync) + { + if (!_isConfigured) + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } + + allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); + useCycles = _useDivisorCycles; + } + + if (allowedMax < 3UL) + { + divisorsExhausted = true; + return true; + } + + ulong processedCount = 0UL; + ulong lastProcessed = 0UL; + bool processedAll = false; + + bool composite = CheckDivisors( + prime, + allowedMax, + useCycles, + out lastProcessed, + out processedAll, + out processedCount); + + if (processedCount > 0UL) + { + lock (_sync) + { + UpdateStatusUnsafe(lastProcessed, processedCount); + } + } + + if (composite) + { + divisorsExhausted = true; + return false; + } + + divisorsExhausted = processedAll || composite; + return true; + } + + 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, + bool useCycles, + out ulong lastProcessed, + out bool processedAll, + out ulong processedCount) + { + lastProcessed = 0UL; + processedCount = 0UL; + processedAll = false; + + if (allowedMax < 3UL) + { + return false; + } + + ulong divisor = 3UL; + while (divisor <= allowedMax) + { + processedCount++; + lastProcessed = divisor; + + MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + byte hit; + if (useCycles) + { + using DivisorCycleCache.Lease cycleLease = DivisorCycleCache.Shared.Acquire(divisor); + ulong divisorCycle = cycleLease.GetCycle(divisor); + hit = CheckDivisor(prime, divisorCycle != 0UL, divisorCycle, divisorData); + } + else + { + hit = CheckDivisor(prime, false, 0UL, divisorData); + } + + if (hit != 0) + { + processedAll = true; + return true; + } + + if (divisor >= allowedMax - 1UL) + { + processedAll = true; + break; + } + + if (divisor > ulong.MaxValue - 2UL) + { + break; + } + + divisor += 2UL; + } + + if (!processedAll) + { + processedAll = divisor > allowedMax; + } + return false; + } + + private static byte CheckDivisor(ulong prime, bool useCycles, ulong divisorCycle, in MontgomeryDivisorData divisorData) + { + ulong modulus = divisorData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return 0; + } + + ulong residue = useCycles && divisorCycle != 0UL + ? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData) + : prime.Pow2MontgomeryMod(divisorData); + + return residue == 1UL ? (byte)1 : (byte)0; + } + + private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount) + { + if (processedCount == 0UL) + { + return; + } + + ulong interval = PerfectNumberConstants.ConsoleInterval; + if (interval == 0UL) + { + _lastStatusDivisor = 0UL; + return; + } + + ulong total = _lastStatusDivisor + processedCount; + _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 static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) + { + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); + } + + return new MontgomeryDivisorData( + modulus, + ComputeMontgomeryNPrime(modulus), + ComputeMontgomeryResidue(1UL, modulus), + ComputeMontgomeryResidue(2UL, modulus)); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); + + private static ulong ComputeMontgomeryNPrime(ulong modulus) + { + ulong inv = modulus; + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + return unchecked(0UL - inv); + } + + private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession + { + private readonly MersenneNumberDivisorByDivisorCpuTester _owner; + private bool _disposed; + + internal DivisorScanSession(MersenneNumberDivisorByDivisorCpuTester owner) + { + _owner = owner; + } + + internal void Reset() + { + _disposed = false; + } + + public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(DivisorScanSession)); + } + + int length = primes.Length; + if (length == 0) + { + return; + } + + ulong modulus = divisor; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits.Clear(); + return; + } + + bool cycleEnabled = _owner._useDivisorCycles && divisorCycle != 0UL; + MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + + for (int i = 0; i < length; i++) + { + ulong prime = primes[i]; + + if (cycleEnabled) + { + ulong remainder = prime % divisorCycle; + if (remainder != 0UL) + { + hits[i] = 0; + continue; + } + } + + ulong residue = cycleEnabled + ? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData) + : prime.Pow2MontgomeryMod(divisorData); + + hits[i] = residue == 1UL ? (byte)1 : (byte)0; + } + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _owner.ReturnSession(this); + } + } } diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index a99d608..15e6316 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -1,428 +1,378 @@ +using System; +using System.Buffers; using System.Collections.Concurrent; using System.Runtime.CompilerServices; +using System.Threading.Tasks; +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; public sealed class DivisorCycleCache { - internal sealed class CycleBlock - { - // TODO: This can be removed. See other comments. - private int _referenceCount; - - internal CycleBlock(int index, ulong start, ulong[] cycles) - { - Index = index; - Start = start; - Cycles = cycles; - End = start + (ulong)cycles.Length - 1UL; - } - - // TODO: Let's use fields instead of the properties below for better performance - internal int Index { get; } - - internal ulong Start { get; } - - internal ulong End { get; } - - internal ulong[] Cycles { get; } - - // TODO: We don't need full reference tracking. It's enough to assign divisor block to a task, instead. We can remove excessive blocks as required. - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal void AddRef() - { - Interlocked.Increment(ref _referenceCount); - } - - // TODO: We don't need full reference tracking. It's enough to assign divisor block to a task, instead. We can remove excessive blocks as required. - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal void Release() - { - Interlocked.Decrement(ref _referenceCount); - } - - // TODO: We don't need full reference tracking. It's enough to assign divisor block to a task, instead. This can be removed - internal bool IsInUse => Volatile.Read(ref _referenceCount) > 0; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal ulong GetCycle(ulong divisor) - { - // TODO: The production code should guarantee that this will never occur and we'll use divisors with their corresponding blocks, only. This can be removed and we can simplify this method to => - if (divisor < Start || divisor > End) - { - return 0UL; - } - - return Cycles[(int)(divisor - Start)]; - } - } - - public struct Lease : IDisposable - { - // TODO: We can probably remove _owner, if we don't need full cycle management which will help GC to clear these out - private readonly DivisorCycleCache? _owner; - private CycleBlock? _block; - - internal Lease(DivisorCycleCache owner, CycleBlock block) - { - _owner = owner; - _block = block; - _block.AddRef(); - } - - // TODO: We're highly over-complicating this. _block is always assigned by the constructor. We don't need to make it nullable. This can be removed. - public readonly bool IsValid => _block is not null; - - // TODO: Once we make _block non-nullable we can simplify these to just fields assigned during the creation for the best performance. - public readonly ulong Start => _block?.Start ?? 0UL; - - public readonly ulong End => _block?.End ?? 0UL; - - public readonly ulong[]? Values => _block?.Cycles; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public readonly ulong GetCycle(ulong divisor) => _block?.GetCycle(divisor) ?? 0UL; - - public void Dispose() - { - // TODO: We don't need full cycle block management. We can rely on GC to do its job. Thanks to that we could possibly make the struct read-only. - CycleBlock? block = _block; - if (block is null) - { - return; - } - - block.Release(); - _owner?.Release(block); - _block = null; - } - } - - private readonly object _sync = new(); - private readonly ConcurrentDictionary> _pending = new(); - // TODO: _baseBlock should be like no other block. Why do we need a separate property for it? It should become current upon initialization - private CycleBlock _baseBlock = null!; - private CycleBlock? _previous; - private CycleBlock? _current; - private CycleBlock? _next; - private bool _initialized; - - // TODO: Let's not use lazy but directly create the instance, instead for better performance. Then we can remove _lazy and just assign the instance to Shared. - private static readonly Lazy _lazy = new(() => new DivisorCycleCache()); - - public static DivisorCycleCache Shared => _lazy.Value; - - private DivisorCycleCache() - { - ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - Initialize(snapshot); - } - - public void ReloadFromCurrentSnapshot() - { - ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - Initialize(snapshot); - } - - private void Initialize(ulong[] snapshot) - { - lock (_sync) - { - _pending.Clear(); - _baseBlock = new CycleBlock(index: 0, start: 0UL, snapshot); - _previous = _baseBlock; - _current = null; - _next = null; - _initialized = true; - StartPrefetchLocked(1); - } - } - - public Lease Acquire(ulong divisor) - { - // TODO: We don't need _initialized. The constructor calls Initialize and makes sure it's always initialized, before use. This check can be removed, too. - if (!_initialized) - { - throw new InvalidOperationException("DivisorCycleCache must be initialized before use."); - } - - // TODO: Do we need Lease instances, if we don't care who is using our blocks once they were given and just remove them internally from future use? - if (divisor <= _baseBlock.End) - { - return new Lease(this, _baseBlock); - } - - CycleBlock block = AcquireDynamicBlock(divisor); - return new Lease(this, block); - } - - private CycleBlock AcquireDynamicBlock(ulong divisor) - { - int blockIndex = GetBlockIndex(divisor); - while (true) - { - Task? taskToAwait = null; - lock (_sync) - { - if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cached)) - { - return cached!; - } - - if (!_pending.TryGetValue(blockIndex, out taskToAwait)) - { - taskToAwait = Task.Run(() => GenerateBlock(blockIndex)); - _pending[blockIndex] = taskToAwait; - } - } - - CycleBlock block = taskToAwait.GetAwaiter().GetResult(); - lock (_sync) - { - _pending.TryRemove(blockIndex, out _); - - if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cachedAfterAwait)) - { - return cachedAfterAwait!; - } - - PromoteBlockLocked(block); - return block; - } - } - } - - private bool TryGetCachedBlockLocked(int blockIndex, out CycleBlock? block) - { - if (_current is { Index: var currentIndex } current && currentIndex == blockIndex) - { - block = current; - return true; - } - - if (_previous is { Index: var previousIndex } previous && previousIndex == blockIndex) - { - block = previous; - return true; - } - - if (_next is { Index: var nextIndex } next && nextIndex == blockIndex) - { - PromoteNextLocked(); - block = _current!; - return true; - } - - block = null; - return false; - } - - private void PromoteNextLocked() - { - if (_next is null) - { - return; - } - - CycleBlock? oldPrevious = _previous; - _previous = _current ?? _baseBlock; - _current = _next; - _next = null; - StartPrefetchLocked(_current!.Index + 1); - TryReleaseBlockLocked(oldPrevious); - } - - private void PromoteBlockLocked(CycleBlock block) - { - if (_current is null) - { - _previous = _baseBlock; - _current = block; - StartPrefetchLocked(block.Index + 1); - return; - } - - if (block.Index > _current.Index) - { - CycleBlock? oldPrevious = _previous; - _previous = _current; - _current = block; - StartPrefetchLocked(block.Index + 1); - TryReleaseBlockLocked(oldPrevious); - return; - } - - if (_previous is null || block.Index >= _previous.Index) - { - CycleBlock? oldPrevious = _previous; - _previous = block; - TryReleaseBlockLocked(oldPrevious); - } - } - - private void StartPrefetchLocked(int blockIndex) - { - // TODO: Why do we need this check? This should never happen in production code. - if (blockIndex <= 0) - { - return; - } - - // TODO: Why do we need this check? This should never happen in production code. - if (_next is { Index: var nextIndex } next && nextIndex == blockIndex) - { - return; - } - - if (_pending.TryGetValue(blockIndex, out Task? existing)) - { - if (existing.IsCompletedSuccessfully) - { - if (_next is null || _next.Index < blockIndex) - { - _next = existing.Result; - } - - _pending.TryRemove(blockIndex, out _); - } - - return; - } - - Task task = Task.Run(() => GenerateBlock(blockIndex)); - _pending[blockIndex] = task; - task.ContinueWith(t => OnPrefetchCompleted(blockIndex, t), TaskScheduler.Default); - } - - private void OnPrefetchCompleted(int blockIndex, Task task) - { - if (!task.IsCompletedSuccessfully) - { - _pending.TryRemove(blockIndex, out _); - return; - } - - lock (_sync) - { - if (_next is null || _next.Index < blockIndex) - { - _next = task.Result; - } - - _pending.TryRemove(blockIndex, out _); - TrimCacheLocked(); - } - } - - private void TrimCacheLocked() - { - if (_previous is null) - { - return; - } - - if (ReferenceEquals(_previous, _baseBlock)) - { - return; - } - - if (_previous.IsInUse) - { - return; - } - - if (_next is not null && _current is not null) - { - _previous = null; - } - } - - private void TryReleaseBlockLocked(CycleBlock? block) - { - if (block is null || ReferenceEquals(block, _baseBlock)) - { - return; - } - - if (ReferenceEquals(block, _current) || ReferenceEquals(block, _next)) - { - return; - } - - if (block.IsInUse) - { - return; - } - - if (ReferenceEquals(block, _previous)) - { - _previous = null; - } - } - - private void Release(CycleBlock block) - { - if (!block.IsInUse) - { - lock (_sync) - { - TryReleaseBlockLocked(block); - } - } - } - - private CycleBlock GenerateBlock(int blockIndex) - { - if (blockIndex == 0) - { - return _baseBlock; - } - - ulong start = GetBlockStart(blockIndex); - int length = GetBlockLength(blockIndex); - ulong[] cycles = new ulong[length]; - ulong divisor = start; - for (int i = 0; i < length; i++) - { - cycles[i] = MersenneDivisorCycles.CalculateCycleLength(divisor); - divisor++; - } - - return new CycleBlock(blockIndex, start, cycles); - } - - private int GetBlockLength(int blockIndex) - { - // TODO: We're over-complicating this again. We should just remember and use the _baseBlock.Cycles.Length. Remember it in a field for better performance during initialization / creation - if (blockIndex == 0) - { - return _baseBlock.Cycles.Length; - } - - return PerfectNumberConstants.MaxQForDivisorCycles; - } - - private ulong GetBlockStart(int blockIndex) - { - // TODO: We're over-complicating this. There should be only 3 blocks kept. We should simply check all 3 to return the correct one, instead of expensive multiplications. - if (blockIndex == 0) - { - return 0UL; - } - - ulong baseEnd = _baseBlock.End; - return baseEnd + 1UL + (ulong)(blockIndex - 1) * (ulong)PerfectNumberConstants.MaxQForDivisorCycles; - } - - private int GetBlockIndex(ulong divisor) - { - // TODO: We're over-complicating this. There should be only 3 blocks kept. We should simply check all 3 to return the correct one, instead of expensive multiplications. - if (divisor <= _baseBlock.End) - { - return 0; - } - - ulong offset = divisor - (_baseBlock.End + 1UL); - return (int)(offset / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1; - } + internal sealed class CycleBlock + { + internal CycleBlock(int index, ulong start, ulong[] cycles) + { + Index = index; + Start = start; + Cycles = cycles; + End = start + (ulong)cycles.Length - 1UL; + } + + internal int Index { get; } + + internal ulong Start { get; } + + internal ulong End { get; } + + internal ulong[] Cycles { get; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal ulong GetCycle(ulong divisor) + { + if (divisor < Start || divisor > End) + { + return 0UL; + } + + return Cycles[(int)(divisor - Start)]; + } + } + + public struct Lease : IDisposable + { + private CycleBlock? _block; + + internal Lease(CycleBlock block) + { + _block = block; + } + + public readonly bool IsValid => _block is not null; + + public readonly ulong Start => _block?.Start ?? 0UL; + + public readonly ulong End => _block?.End ?? 0UL; + + public readonly ulong[]? Values => _block?.Cycles; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public readonly ulong GetCycle(ulong divisor) => _block?.GetCycle(divisor) ?? 0UL; + + public void Dispose() + { + _block = null; + } + } + + private readonly object _sync = new(); + private readonly ConcurrentDictionary> _pending = new(); + private readonly ConcurrentDictionary, ArrayView1D>> _gpuKernelCache = new(AcceleratorReferenceComparer.Instance); + private CycleBlock _baseBlock = null!; + private CycleBlock? _activeBlock; + private CycleBlock? _prefetchedBlock; + private bool _useGpuGeneration = true; + + public static DivisorCycleCache Shared { get; } = new DivisorCycleCache(); + + private DivisorCycleCache() + { + ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + Initialize(snapshot); + } + + public void ReloadFromCurrentSnapshot() + { + ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + Initialize(snapshot); + } + + public void ConfigureGeneratorDevice(bool useGpu) + { + lock (_sync) + { + _useGpuGeneration = useGpu; + } + } + + private void Initialize(ulong[] snapshot) + { + lock (_sync) + { + _pending.Clear(); + _baseBlock = new CycleBlock(index: 0, start: 0UL, snapshot); + _activeBlock = null; + _prefetchedBlock = null; + StartPrefetchLocked(1); + } + } + + public Lease Acquire(ulong divisor) + { + if (divisor <= _baseBlock.End) + { + return new Lease(_baseBlock); + } + + CycleBlock block = AcquireDynamicBlock(divisor); + return new Lease(block); + } + + private CycleBlock AcquireDynamicBlock(ulong divisor) + { + int blockIndex = GetBlockIndex(divisor); + + while (true) + { + Task? pendingTask = null; + lock (_sync) + { + if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cached)) + { + return cached!; + } + + if (!_pending.TryGetValue(blockIndex, out pendingTask)) + { + pendingTask = Task.Run(() => GenerateBlock(blockIndex)); + _pending[blockIndex] = pendingTask; + } + } + + CycleBlock block = pendingTask.GetAwaiter().GetResult(); + lock (_sync) + { + _pending.TryRemove(blockIndex, out _); + + if (TryGetCachedBlockLocked(blockIndex, out CycleBlock? cachedAfterAwait)) + { + return cachedAfterAwait!; + } + + SetActiveBlockLocked(block); + return block; + } + } + } + + private bool TryGetCachedBlockLocked(int blockIndex, out CycleBlock? block) + { + CycleBlock? active = _activeBlock; + if (active is not null && active.Index == blockIndex) + { + block = active; + return true; + } + + CycleBlock? prefetched = _prefetchedBlock; + if (prefetched is not null && prefetched.Index == blockIndex) + { + PromotePrefetchedLocked(); + block = _activeBlock!; + return true; + } + + block = null; + return false; + } + + private void PromotePrefetchedLocked() + { + CycleBlock? prefetched = _prefetchedBlock; + if (prefetched is null) + { + return; + } + + _activeBlock = prefetched; + _prefetchedBlock = null; + StartPrefetchLocked(prefetched.Index + 1); + } + + private void SetActiveBlockLocked(CycleBlock block) + { + _activeBlock = block; + + CycleBlock? prefetched = _prefetchedBlock; + if (prefetched is not null && prefetched.Index != block.Index + 1) + { + _prefetchedBlock = null; + } + + StartPrefetchLocked(block.Index + 1); + } + + private void StartPrefetchLocked(int blockIndex) + { + if (blockIndex <= 0) + { + return; + } + + CycleBlock? prefetched = _prefetchedBlock; + if (prefetched is not null && prefetched.Index == blockIndex) + { + return; + } + + if (_pending.TryGetValue(blockIndex, out Task? existing)) + { + if (existing.IsCompletedSuccessfully) + { + _pending.TryRemove(blockIndex, out _); + AssignPrefetchedBlockLocked(blockIndex, existing.Result); + } + + return; + } + + Task task = Task.Run(() => GenerateBlock(blockIndex)); + _pending[blockIndex] = task; + task.ContinueWith(t => OnPrefetchCompleted(blockIndex, t), TaskScheduler.Default); + } + + private void OnPrefetchCompleted(int blockIndex, Task task) + { + if (!task.IsCompletedSuccessfully) + { + _pending.TryRemove(blockIndex, out _); + return; + } + + lock (_sync) + { + _pending.TryRemove(blockIndex, out _); + AssignPrefetchedBlockLocked(blockIndex, task.Result); + } + } + + private void AssignPrefetchedBlockLocked(int blockIndex, CycleBlock block) + { + CycleBlock? active = _activeBlock; + if (active is not null && blockIndex != active.Index + 1) + { + return; + } + + if (_prefetchedBlock is null || _prefetchedBlock.Index < blockIndex) + { + _prefetchedBlock = block; + } + } + + private CycleBlock GenerateBlock(int blockIndex) + { + if (blockIndex == 0) + { + return _baseBlock; + } + + ulong start = GetBlockStart(blockIndex); + int length = _baseBlock.Cycles.Length; + ulong[] cycles = new ulong[length]; + + if (_useGpuGeneration) + { + ComputeCyclesGpu(start, cycles); + } + else + { + ComputeCyclesCpu(start, cycles); + } + + return new CycleBlock(blockIndex, start, cycles); + } + + private int GetBlockIndex(ulong divisor) + { + if (divisor <= _baseBlock.End) + { + return 0; + } + + ulong offset = divisor - (_baseBlock.End + 1UL); + return (int)(offset / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1; + } + + private ulong GetBlockStart(int blockIndex) + { + if (blockIndex == 0) + { + return 0UL; + } + + ulong baseEnd = _baseBlock.End; + return baseEnd + 1UL + (ulong)(blockIndex - 1) * (ulong)PerfectNumberConstants.MaxQForDivisorCycles; + } + + private void ComputeCyclesCpu(ulong start, ulong[] destination) + { + ulong divisor = start; + for (int i = 0; i < destination.Length; i++) + { + destination[i] = MersenneDivisorCycles.CalculateCycleLength(divisor); + divisor++; + } + } + + private void ComputeCyclesGpu(ulong start, ulong[] destination) + { + var gpuLease = GpuKernelPool.GetKernel(_useGpuGeneration); + var execution = gpuLease.EnterExecutionScope(); + + try + { + Accelerator accelerator = gpuLease.Accelerator; + Action, ArrayView1D> kernel = _gpuKernelCache.GetOrAdd(accelerator, LoadKernel); + + int length = destination.Length; + ulong[] divisors = ArrayPool.Shared.Rent(length); + + try + { + for (int i = 0; i < length; i++) + { + divisors[i] = start + (ulong)i; + } + + using MemoryBuffer1D divisorBuffer = accelerator.Allocate1D(length); + using MemoryBuffer1D resultBuffer = accelerator.Allocate1D(length); + + divisorBuffer.View.CopyFromCPU(ref divisors[0], length); + kernel(length, divisorBuffer.View, resultBuffer.View); + accelerator.Synchronize(); + resultBuffer.View.CopyToCPU(ref destination[0], length); + } + finally + { + ArrayPool.Shared.Return(divisors, clearArray: false); + } + } + finally + { + execution.Dispose(); + gpuLease.Dispose(); + } + } + + private static Action, ArrayView1D> LoadKernel(Accelerator accelerator) + { + return accelerator.LoadAutoGroupedStreamKernel, ArrayView1D>(GpuDivisorCycleKernel); + } + + private static void GpuDivisorCycleKernel(Index1D index, ArrayView1D divisors, ArrayView1D cycles) + { + cycles[index] = MersenneDivisorCycles.CalculateCycleLengthGpu(divisors[index]); + } + + private sealed class AcceleratorReferenceComparer : IEqualityComparer + { + internal static AcceleratorReferenceComparer Instance { get; } = new(); + + public bool Equals(Accelerator? x, Accelerator? y) => ReferenceEquals(x, y); + + public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj); + } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index d295418..9598590 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -11,780 +11,753 @@ namespace PerfectNumbers.Core.Gpu; public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDivisorByDivisorTester { - private int _gpuBatchSize = GpuConstants.ScanBatchSize; - private readonly object _sync = new(); - private ulong _divisorLimit; - private bool _isConfigured; - private ulong _lastStatusDivisor; - private bool _useDivisorCycles; - - public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo) - { - public readonly ulong Modulus = modulus; - public readonly ulong NPrime = nPrime; - public readonly ulong MontgomeryOne = montgomeryOne; - public readonly ulong MontgomeryTwo = montgomeryTwo; - } - - private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); - private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView>> _cycleKernelCache = new(); - private readonly ConcurrentBag _sessionPool = new(); - - private Action, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); - private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputePrimeExponentKernel)); - private Action, ArrayView> GetCycleKernel(Accelerator accelerator) => - _cycleKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckPrimeCycleKernel)); - - public int GpuBatchSize - { - get => _gpuBatchSize; - set => _gpuBatchSize = Math.Max(1, value); - } - - int IMersenneNumberDivisorByDivisorTester.BatchSize - { - get => GpuBatchSize; - set => GpuBatchSize = value; - } - - public bool UseDivisorCycles - { - get => _useDivisorCycles; - set => _useDivisorCycles = value; - } - - public void ConfigureFromMaxPrime(ulong maxPrime) - { - lock (_sync) - { - _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); - _lastStatusDivisor = 0UL; - _isConfigured = true; - } - } - - public bool IsPrime(ulong prime, out bool divisorsExhausted) - { - ulong allowedMax; - bool useCycles; - int batchCapacity; - - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); - useCycles = _useDivisorCycles; - batchCapacity = _gpuBatchSize; - } - - if (allowedMax < 3UL) - { - divisorsExhausted = true; - return true; - } - - bool composite; - bool coveredRange; - ulong processedCount; - ulong lastProcessed; - - var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); - var accelerator = gpuLease.Accelerator; - var kernel = GetKernel(accelerator); - BatchResources resources = RentBatchResources(accelerator, batchCapacity); - - try - { - composite = CheckDivisors( - prime, - allowedMax, - useCycles, - accelerator, - kernel, - resources.DivisorsBuffer, - resources.HitsBuffer, - resources.Divisors, - resources.Hits, - resources.DivisorData, - out lastProcessed, - out coveredRange, - out processedCount - ); - } - finally - { - ReturnBatchResources(accelerator, resources); - gpuLease.Dispose(); - } - - if (processedCount > 0UL) - { - lock (_sync) - { - UpdateStatusUnsafe(lastProcessed, processedCount); - } - } - - if (composite) - { - divisorsExhausted = true; - return false; - } - - divisorsExhausted = coveredRange; - return true; - } - - private static bool CheckDivisors( - ulong prime, - ulong allowedMax, - bool useCycles, - Accelerator accelerator, - Action, ArrayView> kernel, - MemoryBuffer1D divisorsBuffer, - MemoryBuffer1D hitsBuffer, - ulong[] divisors, - byte[] hits, - MontgomeryDivisorData[] divisorData, - out ulong lastProcessed, - out bool coveredRange, - out ulong processedCount) - { - lastProcessed = 0UL; - processedCount = 0UL; - - if (allowedMax < 3UL) - { - coveredRange = true; - return false; - } - - int batchCapacity = (int)divisorsBuffer.Length; - - bool composite = false; - bool processedAll = false; - ulong currentDivisor, divisor = 3UL, next, nextDivisor; - int batchSize, i; - bool reachedEndInBatch; - - Span divisorDataSpan; - Span divisorSpan; - Span hitsSpan; - ArrayView1D divisorView; - ArrayView1D hitsView; - - while (divisor <= allowedMax) - { - batchSize = 0; - reachedEndInBatch = false; - - while (batchSize < batchCapacity && divisor <= allowedMax) - { - currentDivisor = divisor; - nextDivisor = currentDivisor + 2UL; - processedCount++; - - divisors[batchSize++] = currentDivisor; - lastProcessed = currentDivisor; - next = nextDivisor; - - // This will be only true when we exceed ulong.MaxValue - if (next <= divisor) - { - reachedEndInBatch = true; - break; - } - - if (next > allowedMax) - { - divisor = next; - reachedEndInBatch = true; - break; - } - - divisor = next; - } - - if (batchSize == 0) - { - if (reachedEndInBatch) - { - processedAll = true; - } - - break; - } - - divisorDataSpan = divisorData.AsSpan(0, batchSize); - divisorView = divisorsBuffer.View.SubView(0, batchSize); - hitsView = hitsBuffer.View.SubView(0, batchSize); - - divisorSpan = divisors.AsSpan(0, batchSize); - hitsSpan = hits.AsSpan(0, batchSize); - - for (i = 0; i < batchSize; i++) - { - divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]); - } - - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); - kernel(batchSize, prime, useCycles ? (byte)1 : (byte)0, divisorView, hitsView); - accelerator.Synchronize(); - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); - - for (i = 0; i < batchSize; i++) - { - if (hitsSpan[i] != 0) - { - composite = true; - lastProcessed = divisorSpan[i]; - break; - } - } - - if (!composite) - { - lastProcessed = divisorSpan[batchSize - 1]; - } - if (composite) - { - break; - } - - if (reachedEndInBatch) - { - processedAll = true; - break; - } - } - - coveredRange = composite || processedAll || divisor > allowedMax; - - return composite; - } - - 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); - } - - // TODO: Since lastProcessed is unused, let's remove it. - private void UpdateStatusUnsafe(ulong lastProcessed, ulong processedCount) - { - if (processedCount == 0UL) - { - return; - } - - ulong interval = PerfectNumberConstants.ConsoleInterval; - if (interval == 0UL) - { - _lastStatusDivisor = 0UL; - return; - } - - // Can we avoid the modulo calculation here? For the higher divisors it may become very expensive. Do we need to calculate the total, too? - ulong total = _lastStatusDivisor + processedCount; - _lastStatusDivisor = total % interval; - } - - private static void CheckKernel(Index1D index, ulong prime, byte useCycleChecks, ArrayView divisors, ArrayView hits) - { - MontgomeryDivisorData divisor = divisors[index]; - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } - - ulong exponent = prime; - if (useCycleChecks != 0) - { - ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); - if (cycle == 0UL) - { - hits[index] = 0; - return; - } - // TODO: Can we move this back before the GPU call and let CPU handle it? We significantly slowed down after this was moved to GPU. - ulong remainder = prime % cycle; - if (remainder != 0UL) - { - hits[index] = 0; - return; - } - - exponent = remainder; - } - - hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; - } - - public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession - { - private readonly MersenneNumberDivisorByDivisorGpuTester _owner; - private readonly GpuContextPool.GpuContextLease _lease; - private readonly Accelerator _accelerator; - private readonly Action, ArrayView> _kernel; - private readonly Action, ArrayView> _cycleKernel; - private MemoryBuffer1D _exponentsBuffer; - private MemoryBuffer1D _resultsBuffer; - private MemoryBuffer1D _cycleMaskBuffer; - private ulong[] _hostBuffer; - private ulong[] _primeBuffer; - private int[] _positionBuffer; - private byte[] _cycleMaskHost; - private int _capacity; - private bool _disposed; - - internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) - { - _owner = owner; - _lease = GpuContextPool.RentPreferred(preferCpu: false); - _accelerator = _lease.Accelerator; - _kernel = owner.GetKernelByPrimeExponent(_accelerator); - _cycleKernel = owner.GetCycleKernel(_accelerator); - _capacity = Math.Max(1, owner._gpuBatchSize); - _exponentsBuffer = _accelerator.Allocate1D(_capacity); - _resultsBuffer = _accelerator.Allocate1D(_capacity); - _cycleMaskBuffer = _accelerator.Allocate1D(_capacity); - _hostBuffer = ArrayPool.Shared.Rent(_capacity); - _primeBuffer = ArrayPool.Shared.Rent(_capacity); - _positionBuffer = ArrayPool.Shared.Rent(_capacity); - _cycleMaskHost = ArrayPool.Shared.Rent(_capacity); - } - - internal void Reset() - { - _disposed = false; - } - - private void EnsureCapacity(int requiredCapacity) - { - if (requiredCapacity <= _capacity) - { - return; - } - - _exponentsBuffer.Dispose(); - _resultsBuffer.Dispose(); - _cycleMaskBuffer.Dispose(); - ArrayPool.Shared.Return(_hostBuffer, clearArray: false); - ArrayPool.Shared.Return(_primeBuffer, clearArray: false); - ArrayPool.Shared.Return(_positionBuffer, clearArray: false); - ArrayPool.Shared.Return(_cycleMaskHost, clearArray: false); - - _capacity = requiredCapacity; - _exponentsBuffer = _accelerator.Allocate1D(_capacity); - _resultsBuffer = _accelerator.Allocate1D(_capacity); - _cycleMaskBuffer = _accelerator.Allocate1D(_capacity); - _hostBuffer = ArrayPool.Shared.Rent(_capacity); - _primeBuffer = ArrayPool.Shared.Rent(_capacity); - _positionBuffer = ArrayPool.Shared.Rent(_capacity); - _cycleMaskHost = ArrayPool.Shared.Rent(_capacity); - } - - public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) - { - if (_disposed) - { - throw new ObjectDisposedException(nameof(DivisorScanSession)); - } - - int primesLength = primes.Length; - if (primesLength == 0) - { - return; - } - - MersenneNumberDivisorByDivisorGpuTester owner = _owner; - int gpuBatchSize = owner._gpuBatchSize; - EnsureCapacity(gpuBatchSize); - - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return; - } - - ulong cycle = divisorCycle; - bool cycleEnabled = owner._useDivisorCycles && cycle != 0UL; - - Accelerator accelerator = _accelerator; - Action, ArrayView> kernel = _kernel; - Action, ArrayView> cycleKernel = _cycleKernel; - ArrayView1D exponentsView = _exponentsBuffer.View; - ArrayView1D resultsView = _resultsBuffer.View; - ArrayView1D cycleMaskViewFull = _cycleMaskBuffer.View; - - int offset = 0; - Span hostSpan = _hostBuffer.AsSpan(); - Span primeSpan = _primeBuffer.AsSpan(); - Span positionSpan = _positionBuffer.AsSpan(); - bool hasPreviousResidue = false; - ulong previousResidue = 0UL; - ulong previousPrime = 0UL; - - while (offset < primesLength) - { - int batchSize = Math.Min(gpuBatchSize, primesLength - offset); - ReadOnlySpan primesSlice = primes.Slice(offset, batchSize); - Span hitsSlice = hits.Slice(offset, batchSize); - - int computeCount = 0; - Span cycleMaskSpan = default; - - if (cycleEnabled) - { - // TODO: Can we move this back before the GPU call and let CPU handle it? We significantly slowed down after this was moved to GPU. - cycleMaskSpan = _cycleMaskHost.AsSpan(0, batchSize); - ArrayView1D primeInputView = _resultsBuffer.View.SubView(0, batchSize); - primeInputView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize); - ArrayView1D cycleMaskView = cycleMaskViewFull.SubView(0, batchSize); - cycleKernel(batchSize, cycle, primeInputView, cycleMaskView); - accelerator.Synchronize(); - cycleMaskView.CopyToCPU(ref MemoryMarshal.GetReference(cycleMaskSpan), batchSize); - } - - for (int i = 0; i < batchSize; i++) - { - if (cycleEnabled && cycleMaskSpan[i] == 0) - { - hitsSlice[i] = 0; - continue; - } - - // TODO: Avoid one-time use variables like prime here. - ulong prime = primesSlice[i]; - primeSpan[computeCount] = prime; - positionSpan[computeCount] = i; - computeCount++; - } - - if (computeCount > 0) - { - Span exponentSlice = hostSpan[..computeCount]; - ulong currentPrime = hasPreviousResidue ? previousPrime : 0UL; - bool hasDeltaSource = hasPreviousResidue; - - for (int i = 0; i < computeCount; i++) - { - ulong prime = primeSpan[i]; - ulong exponentValue; - if (hasDeltaSource) - { - exponentValue = prime - currentPrime; - } - else - { - exponentValue = prime; - hasDeltaSource = true; - } - - if (cycleEnabled) - { - exponentValue %= cycle; - } - - exponentSlice[i] = exponentValue; - - currentPrime = prime; - } - - ArrayView1D exponentView = exponentsView.SubView(0, computeCount); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - ArrayView1D resultView = resultsView.SubView(0, computeCount); - kernel(computeCount, divisorData, exponentView, resultView); - accelerator.Synchronize(); - resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - - ulong residueAccumulator = previousResidue; - ulong lastPrime = previousPrime; - bool hasAccumulator = hasPreviousResidue; - - for (int i = 0; i < computeCount; i++) - { - // TODO: Avoid one-time use variables like prime and position here. - int position = positionSpan[i]; - ulong stepResidue = exponentSlice[i]; - ulong prime = primeSpan[i]; - ulong residue = hasAccumulator ? MultiplyMod(residueAccumulator, stepResidue, modulus) : stepResidue; - hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0; - exponentSlice[i] = residue; - residueAccumulator = residue; - lastPrime = prime; - hasAccumulator = true; - } - - if (hasAccumulator) - { - previousResidue = residueAccumulator; - previousPrime = lastPrime; - hasPreviousResidue = true; - } - } - - offset += batchSize; - } - } - - public void Dispose() - { - if (_disposed) - { - return; - } - - _disposed = true; - _owner.ReturnSession(this); - } - - private static ulong MultiplyMod(ulong left, ulong right, ulong modulus) - { - if (modulus == 0UL) - { - return 0UL; - } - - UInt128 product = (UInt128)left * right; - return (ulong)(product % modulus); - } - - } - - internal void ReturnSession(DivisorScanSession session) - { - _sessionPool.Add(session); - } - - private static void CheckPrimeCycleKernel(Index1D index, ulong cycle, ArrayView primes, ArrayView mask) - { - if (cycle == 0UL) - { - mask[index] = 1; - return; - } - - mask[index] = primes[index] % cycle == 0UL ? (byte)1 : (byte)0; - } - - private static void ComputePrimeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) - { - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[index] = 0UL; - return; - } - - ulong exponent = exponents[index]; - results[index] = exponent.Pow2MontgomeryMod(divisor); - } - - private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) - { - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return new(modulus, 0UL, 0UL, 0UL); - } - - return new( - modulus, - ComputeMontgomeryNPrime(modulus), - ComputeMontgomeryResidue(1UL, modulus), - ComputeMontgomeryResidue(2UL, modulus) - ); - } - - private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); - - private static ulong ComputeMontgomeryNPrime(ulong modulus) - { - ulong inv = modulus; - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - inv *= unchecked(2UL - modulus * inv); - return unchecked(0UL - inv); - } - - public ulong DivisorLimit - { - get - { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - return _divisorLimit; - } - } - } - - private BatchResources RentBatchResources(Accelerator accelerator, int capacity) - { - var bag = _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()); - if (bag.TryTake(out BatchResources? resources)) - { - resources.EnsureCapacity(capacity); - return resources; - } - - return new BatchResources(accelerator, capacity); - } - - private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) - { - _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()).Add(resources); - } - - private sealed class BatchResources : IDisposable - { - private readonly Accelerator _accelerator; - - internal BatchResources(Accelerator accelerator, int capacity) - { - _accelerator = accelerator; - DivisorsBuffer = accelerator.Allocate1D(capacity); - HitsBuffer = accelerator.Allocate1D(capacity); - Divisors = ArrayPool.Shared.Rent(capacity); - Hits = ArrayPool.Shared.Rent(capacity); - DivisorData = ArrayPool.Shared.Rent(capacity); - Capacity = capacity; - } - - internal MemoryBuffer1D DivisorsBuffer { get; private set; } - - internal MemoryBuffer1D HitsBuffer { get; private set; } - - internal ulong[] Divisors { get; private set; } - - internal byte[] Hits { get; private set; } - - internal MontgomeryDivisorData[] DivisorData { get; private set; } - - internal int Capacity { get; private set; } - - public void EnsureCapacity(int requiredCapacity) - { - if (requiredCapacity <= Capacity) - { - return; - } - - Resize(requiredCapacity); - } - - private void Resize(int newCapacity) - { - DivisorsBuffer.Dispose(); - HitsBuffer.Dispose(); - ArrayPool.Shared.Return(Divisors, clearArray: false); - ArrayPool.Shared.Return(Hits, clearArray: false); - ArrayPool.Shared.Return(DivisorData, clearArray: false); - - DivisorsBuffer = _accelerator.Allocate1D(newCapacity); - HitsBuffer = _accelerator.Allocate1D(newCapacity); - Divisors = ArrayPool.Shared.Rent(newCapacity); - Hits = ArrayPool.Shared.Rent(newCapacity); - DivisorData = ArrayPool.Shared.Rent(newCapacity); - Capacity = newCapacity; - } - - public void Dispose() - { - DivisorsBuffer.Dispose(); - HitsBuffer.Dispose(); - ArrayPool.Shared.Return(Divisors, clearArray: false); - ArrayPool.Shared.Return(Hits, clearArray: false); - ArrayPool.Shared.Return(DivisorData, clearArray: false); - } - } - - private sealed class AcceleratorReferenceComparer : IEqualityComparer - { - public static AcceleratorReferenceComparer Instance { get; } = new(); - - public bool Equals(Accelerator? x, Accelerator? y) => ReferenceEquals(x, y); - - public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj); - } - - 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 DivisorScanSession 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); - } - } - - IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() - { - return CreateDivisorSession(); - } + private int _gpuBatchSize = GpuConstants.ScanBatchSize; + private readonly object _sync = new(); + private ulong _divisorLimit; + private bool _isConfigured; + private ulong _lastStatusDivisor; + private bool _useDivisorCycles; + + public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo) + { + public readonly ulong Modulus = modulus; + public readonly ulong NPrime = nPrime; + public readonly ulong MontgomeryOne = montgomeryOne; + public readonly ulong MontgomeryTwo = montgomeryTwo; + } + + private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); + private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView>> _cycleKernelCache = new(); + private readonly ConcurrentBag _sessionPool = new(); + + private Action, ArrayView> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); + private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => + _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputePrimeExponentKernel)); + private Action, ArrayView> GetCycleKernel(Accelerator accelerator) => + _cycleKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckPrimeCycleKernel)); + + public int GpuBatchSize + { + get => _gpuBatchSize; + set => _gpuBatchSize = Math.Max(1, value); + } + + int IMersenneNumberDivisorByDivisorTester.BatchSize + { + get => GpuBatchSize; + set => GpuBatchSize = value; + } + + public bool UseDivisorCycles + { + get => _useDivisorCycles; + set => _useDivisorCycles = value; + } + + public void ConfigureFromMaxPrime(ulong maxPrime) + { + lock (_sync) + { + _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); + _lastStatusDivisor = 0UL; + _isConfigured = true; + } + } + + public bool IsPrime(ulong prime, out bool divisorsExhausted) + { + ulong allowedMax; + bool useCycles; + int batchCapacity; + + lock (_sync) + { + if (!_isConfigured) + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } + + allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); + useCycles = _useDivisorCycles; + batchCapacity = _gpuBatchSize; + } + + if (allowedMax < 3UL) + { + divisorsExhausted = true; + return true; + } + + bool composite; + bool coveredRange; + ulong processedCount; + ulong lastProcessed; + + var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); + var accelerator = gpuLease.Accelerator; + var kernel = GetKernel(accelerator); + BatchResources resources = RentBatchResources(accelerator, batchCapacity); + + try + { + composite = CheckDivisors( + prime, + allowedMax, + useCycles, + accelerator, + kernel, + resources.DivisorsBuffer, + resources.HitsBuffer, + resources.Divisors, + resources.Hits, + resources.DivisorData, + out lastProcessed, + out coveredRange, + out processedCount + ); + } + finally + { + ReturnBatchResources(accelerator, resources); + gpuLease.Dispose(); + } + + if (processedCount > 0UL) + { + lock (_sync) + { + UpdateStatusUnsafe(processedCount); + } + } + + if (composite) + { + divisorsExhausted = true; + return false; + } + + divisorsExhausted = coveredRange; + return true; + } + + private static bool CheckDivisors( + ulong prime, + ulong allowedMax, + bool useCycles, + Accelerator accelerator, + Action, ArrayView> kernel, + MemoryBuffer1D divisorsBuffer, + MemoryBuffer1D hitsBuffer, + ulong[] divisors, + byte[] hits, + MontgomeryDivisorData[] divisorData, + out ulong lastProcessed, + out bool coveredRange, + out ulong processedCount) + { + lastProcessed = 0UL; + processedCount = 0UL; + + if (allowedMax < 3UL) + { + coveredRange = true; + return false; + } + + int batchCapacity = (int)divisorsBuffer.Length; + + bool composite = false; + bool processedAll = false; + ulong currentDivisor, divisor = 3UL, next, nextDivisor; + int batchSize, i; + bool reachedEndInBatch; + bool requiresCycleKernelCheck; + + Span divisorDataSpan; + Span divisorSpan; + Span hitsSpan; + ArrayView1D divisorView; + ArrayView1D hitsView; + + while (divisor <= allowedMax) + { + batchSize = 0; + reachedEndInBatch = false; + requiresCycleKernelCheck = false; + + while (batchSize < batchCapacity && divisor <= allowedMax) + { + currentDivisor = divisor; + nextDivisor = currentDivisor + 2UL; + processedCount++; + + bool includeDivisor = true; + if (useCycles) + { + using DivisorCycleCache.Lease cycleLease = DivisorCycleCache.Shared.Acquire(currentDivisor); + ulong cycle = cycleLease.GetCycle(currentDivisor); + if (cycle == 0UL) + { + requiresCycleKernelCheck = true; + } + else if (prime % cycle != 0UL) + { + includeDivisor = false; + } + } + + if (includeDivisor) + { + divisors[batchSize++] = currentDivisor; + } + + lastProcessed = currentDivisor; + next = nextDivisor; + + if (next <= divisor) + { + reachedEndInBatch = true; + break; + } + + if (next > allowedMax) + { + divisor = next; + reachedEndInBatch = true; + break; + } + + divisor = next; + } + + if (batchSize == 0) + { + if (reachedEndInBatch) + { + processedAll = true; + break; + } + + continue; + } + + divisorDataSpan = divisorData.AsSpan(0, batchSize); + divisorView = divisorsBuffer.View.SubView(0, batchSize); + hitsView = hitsBuffer.View.SubView(0, batchSize); + + divisorSpan = divisors.AsSpan(0, batchSize); + hitsSpan = hits.AsSpan(0, batchSize); + + for (i = 0; i < batchSize; i++) + { + divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]); + } + + divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); + byte useCycleChecks = useCycles && requiresCycleKernelCheck ? (byte)1 : (byte)0; + kernel(batchSize, prime, useCycleChecks, divisorView, hitsView); + accelerator.Synchronize(); + hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); + + for (i = 0; i < batchSize; i++) + { + if (hitsSpan[i] != 0) + { + composite = true; + lastProcessed = divisorSpan[i]; + break; + } + } + + if (!composite) + { + lastProcessed = divisorSpan[batchSize - 1]; + } + if (composite) + { + break; + } + + if (reachedEndInBatch) + { + processedAll = true; + break; + } + } + + coveredRange = composite || processedAll || divisor > allowedMax; + + return composite; + } + + 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 void UpdateStatusUnsafe(ulong processedCount) + { + if (processedCount == 0UL) + { + return; + } + + ulong interval = PerfectNumberConstants.ConsoleInterval; + if (interval == 0UL) + { + _lastStatusDivisor = 0UL; + return; + } + + // Can we avoid the modulo calculation here? For the higher divisors it may become very expensive. Do we need to calculate the total, too? + ulong total = _lastStatusDivisor + processedCount; + _lastStatusDivisor = total % interval; + } + + private static void CheckKernel(Index1D index, ulong prime, byte useCycleChecks, ArrayView divisors, ArrayView hits) + { + MontgomeryDivisorData divisor = divisors[index]; + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits[index] = 0; + return; + } + + ulong exponent = prime; + if (useCycleChecks != 0) + { + ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); + if (cycle == 0UL) + { + hits[index] = 0; + return; + } + ulong remainder = prime % cycle; + if (remainder != 0UL) + { + hits[index] = 0; + return; + } + + exponent = remainder; + } + + hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; + } + + public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession + { + private readonly MersenneNumberDivisorByDivisorGpuTester _owner; + private readonly GpuContextPool.GpuContextLease _lease; + private readonly Accelerator _accelerator; + private readonly Action, ArrayView> _kernel; + private readonly Action, ArrayView> _cycleKernel; + private MemoryBuffer1D _exponentsBuffer; + private MemoryBuffer1D _resultsBuffer; + private MemoryBuffer1D _cycleMaskBuffer; + private ulong[] _hostBuffer; + private ulong[] _primeBuffer; + private int[] _positionBuffer; + private byte[] _cycleMaskHost; + private int _capacity; + private bool _disposed; + + internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) + { + _owner = owner; + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = owner.GetKernelByPrimeExponent(_accelerator); + _cycleKernel = owner.GetCycleKernel(_accelerator); + _capacity = Math.Max(1, owner._gpuBatchSize); + _exponentsBuffer = _accelerator.Allocate1D(_capacity); + _resultsBuffer = _accelerator.Allocate1D(_capacity); + _cycleMaskBuffer = _accelerator.Allocate1D(_capacity); + _hostBuffer = ArrayPool.Shared.Rent(_capacity); + _primeBuffer = ArrayPool.Shared.Rent(_capacity); + _positionBuffer = ArrayPool.Shared.Rent(_capacity); + _cycleMaskHost = ArrayPool.Shared.Rent(_capacity); + } + + internal void Reset() + { + _disposed = false; + } + + private void EnsureCapacity(int requiredCapacity) + { + if (requiredCapacity <= _capacity) + { + return; + } + + _exponentsBuffer.Dispose(); + _resultsBuffer.Dispose(); + _cycleMaskBuffer.Dispose(); + ArrayPool.Shared.Return(_hostBuffer, clearArray: false); + ArrayPool.Shared.Return(_primeBuffer, clearArray: false); + ArrayPool.Shared.Return(_positionBuffer, clearArray: false); + ArrayPool.Shared.Return(_cycleMaskHost, clearArray: false); + + _capacity = requiredCapacity; + _exponentsBuffer = _accelerator.Allocate1D(_capacity); + _resultsBuffer = _accelerator.Allocate1D(_capacity); + _cycleMaskBuffer = _accelerator.Allocate1D(_capacity); + _hostBuffer = ArrayPool.Shared.Rent(_capacity); + _primeBuffer = ArrayPool.Shared.Rent(_capacity); + _positionBuffer = ArrayPool.Shared.Rent(_capacity); + _cycleMaskHost = ArrayPool.Shared.Rent(_capacity); + } + + public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(DivisorScanSession)); + } + + int primesLength = primes.Length; + if (primesLength == 0) + { + return; + } + + MersenneNumberDivisorByDivisorGpuTester owner = _owner; + int gpuBatchSize = owner._gpuBatchSize; + EnsureCapacity(gpuBatchSize); + + MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + ulong modulus = divisorData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return; + } + + ulong cycle = divisorCycle; + bool cycleEnabled = owner._useDivisorCycles && cycle != 0UL; + + Accelerator accelerator = _accelerator; + Action, ArrayView> kernel = _kernel; + Action, ArrayView> cycleKernel = _cycleKernel; + ArrayView1D exponentsView = _exponentsBuffer.View; + ArrayView1D resultsView = _resultsBuffer.View; + ArrayView1D cycleMaskViewFull = _cycleMaskBuffer.View; + + int offset = 0; + Span hostSpan = _hostBuffer.AsSpan(); + Span primeSpan = _primeBuffer.AsSpan(); + Span positionSpan = _positionBuffer.AsSpan(); + + while (offset < primesLength) + { + int batchSize = Math.Min(gpuBatchSize, primesLength - offset); + ReadOnlySpan primesSlice = primes.Slice(offset, batchSize); + Span hitsSlice = hits.Slice(offset, batchSize); + + int computeCount = 0; + + if (cycleEnabled) + { + Span cycleMaskSpan = _cycleMaskHost.AsSpan(0, batchSize); + ArrayView1D primeInputView = _resultsBuffer.View.SubView(0, batchSize); + primeInputView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize); + ArrayView1D cycleMaskView = cycleMaskViewFull.SubView(0, batchSize); + cycleKernel(batchSize, cycle, primeInputView, cycleMaskView); + accelerator.Synchronize(); + cycleMaskView.CopyToCPU(ref MemoryMarshal.GetReference(cycleMaskSpan), batchSize); + + for (int i = 0; i < batchSize; i++) + { + if (cycleMaskSpan[i] == 0) + { + hitsSlice[i] = 0; + continue; + } + + primeSpan[computeCount] = primesSlice[i]; + positionSpan[computeCount] = i; + computeCount++; + } + } + else + { + for (int i = 0; i < batchSize; i++) + { + primeSpan[computeCount] = primesSlice[i]; + positionSpan[computeCount] = i; + computeCount++; + } + } + + if (computeCount > 0) + { + Span exponentSlice = hostSpan[..computeCount]; + for (int i = 0; i < computeCount; i++) + { + ulong prime = primeSpan[i]; + exponentSlice[i] = cycleEnabled ? prime % cycle : prime; + } + + ArrayView1D exponentView = exponentsView.SubView(0, computeCount); + exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); + ArrayView1D resultView = resultsView.SubView(0, computeCount); + kernel(computeCount, divisorData, exponentView, resultView); + accelerator.Synchronize(); + resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); + + for (int i = 0; i < computeCount; i++) + { + int position = positionSpan[i]; + ulong residue = exponentSlice[i]; + hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0; + } + } + + offset += batchSize; + } + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _owner.ReturnSession(this); + } + + + } + + internal void ReturnSession(DivisorScanSession session) + { + _sessionPool.Add(session); + } + + private static void CheckPrimeCycleKernel(Index1D index, ulong cycle, ArrayView primes, ArrayView mask) + { + if (cycle == 0UL) + { + mask[index] = 1; + return; + } + + mask[index] = primes[index] % cycle == 0UL ? (byte)1 : (byte)0; + } + + private static void ComputePrimeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) + { + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + results[index] = 0UL; + return; + } + + ulong exponent = exponents[index]; + results[index] = exponent.Pow2MontgomeryMod(divisor); + } + + private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) + { + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return new(modulus, 0UL, 0UL, 0UL); + } + + return new( + modulus, + ComputeMontgomeryNPrime(modulus), + ComputeMontgomeryResidue(1UL, modulus), + ComputeMontgomeryResidue(2UL, modulus) + ); + } + + private static ulong ComputeMontgomeryResidue(ulong value, ulong modulus) => (ulong)((UInt128)value * (UInt128.One << 64) % modulus); + + private static ulong ComputeMontgomeryNPrime(ulong modulus) + { + ulong inv = modulus; + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + inv *= unchecked(2UL - modulus * inv); + return unchecked(0UL - inv); + } + + public ulong DivisorLimit + { + get + { + lock (_sync) + { + if (!_isConfigured) + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } + + return _divisorLimit; + } + } + } + + private BatchResources RentBatchResources(Accelerator accelerator, int capacity) + { + var bag = _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()); + if (bag.TryTake(out BatchResources? resources)) + { + resources.EnsureCapacity(capacity); + return resources; + } + + return new BatchResources(accelerator, capacity); + } + + private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) + { + _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()).Add(resources); + } + + private sealed class BatchResources : IDisposable + { + private readonly Accelerator _accelerator; + + internal BatchResources(Accelerator accelerator, int capacity) + { + _accelerator = accelerator; + DivisorsBuffer = accelerator.Allocate1D(capacity); + HitsBuffer = accelerator.Allocate1D(capacity); + Divisors = ArrayPool.Shared.Rent(capacity); + Hits = ArrayPool.Shared.Rent(capacity); + DivisorData = ArrayPool.Shared.Rent(capacity); + Capacity = capacity; + } + + internal MemoryBuffer1D DivisorsBuffer { get; private set; } + + internal MemoryBuffer1D HitsBuffer { get; private set; } + + internal ulong[] Divisors { get; private set; } + + internal byte[] Hits { get; private set; } + + internal MontgomeryDivisorData[] DivisorData { get; private set; } + + internal int Capacity { get; private set; } + + public void EnsureCapacity(int requiredCapacity) + { + if (requiredCapacity <= Capacity) + { + return; + } + + Resize(requiredCapacity); + } + + private void Resize(int newCapacity) + { + DivisorsBuffer.Dispose(); + HitsBuffer.Dispose(); + ArrayPool.Shared.Return(Divisors, clearArray: false); + ArrayPool.Shared.Return(Hits, clearArray: false); + ArrayPool.Shared.Return(DivisorData, clearArray: false); + + DivisorsBuffer = _accelerator.Allocate1D(newCapacity); + HitsBuffer = _accelerator.Allocate1D(newCapacity); + Divisors = ArrayPool.Shared.Rent(newCapacity); + Hits = ArrayPool.Shared.Rent(newCapacity); + DivisorData = ArrayPool.Shared.Rent(newCapacity); + Capacity = newCapacity; + } + + public void Dispose() + { + DivisorsBuffer.Dispose(); + HitsBuffer.Dispose(); + ArrayPool.Shared.Return(Divisors, clearArray: false); + ArrayPool.Shared.Return(Hits, clearArray: false); + ArrayPool.Shared.Return(DivisorData, clearArray: false); + } + } + + private sealed class AcceleratorReferenceComparer : IEqualityComparer + { + public static AcceleratorReferenceComparer Instance { get; } = new(); + + public bool Equals(Accelerator? x, Accelerator? y) => ReferenceEquals(x, y); + + public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj); + } + + 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 DivisorScanSession 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); + } + } + + IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() + { + return CreateDivisorSession(); + } } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 3f2579c..b02b507 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -246,6 +246,23 @@ public static class ULongExtensions return result.MontgomeryMultiply(1UL, modulus, nPrime); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor) + { + if (cycleLength == 0UL) + { + return exponent.Pow2MontgomeryMod(divisor); + } + + ulong rotationCount = exponent % cycleLength; + if (rotationCount == 0UL) + { + return 1UL; + } + + return rotationCount.Pow2MontgomeryMod(divisor); + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static UInt128 PowMod(this ulong exponent, UInt128 modulus) {