From bc4e5269467b57b27a22ca7f55a2f85079aa2cbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 25 Sep 2025 01:07:58 +0200 Subject: [PATCH] TODOs added --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 1097 +++++++++-------- ...MersenneNumberDivisorByDivisorGpuTester.cs | 185 +-- 2 files changed, 643 insertions(+), 639 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index e9b2cea..1965ad1 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -8,555 +8,556 @@ 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; + 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); + } + } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index b34d441..d295418 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -39,17 +39,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv 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); - } + public int GpuBatchSize + { + get => _gpuBatchSize; + set => _gpuBatchSize = Math.Max(1, value); + } - int IMersenneNumberDivisorByDivisorTester.BatchSize - { - get => GpuBatchSize; - set => GpuBatchSize = value; - } + int IMersenneNumberDivisorByDivisorTester.BatchSize + { + get => GpuBatchSize; + set => GpuBatchSize = value; + } public bool UseDivisorCycles { @@ -301,6 +301,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv 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) @@ -315,47 +316,46 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv 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; - - // Removed noisy console output; status is tracked internally only. } - 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; - } + 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 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; + } - ulong remainder = prime % cycle; - if (remainder != 0UL) - { - hits[index] = 0; - return; - } + exponent = remainder; + } - exponent = remainder; - } + hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; + } - hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; - } - - public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession - { + public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession + { private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; @@ -442,8 +442,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return; } - ulong cycle = divisorCycle; - bool cycleEnabled = owner._useDivisorCycles && cycle != 0UL; + ulong cycle = divisorCycle; + bool cycleEnabled = owner._useDivisorCycles && cycle != 0UL; Accelerator accelerator = _accelerator; Action, ArrayView> kernel = _kernel; @@ -469,18 +469,19 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv int computeCount = 0; Span cycleMaskSpan = default; - if (cycleEnabled) - { - 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); - } + 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++) + for (int i = 0; i < batchSize; i++) { if (cycleEnabled && cycleMaskSpan[i] == 0) { @@ -488,6 +489,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv continue; } + // TODO: Avoid one-time use variables like prime here. ulong prime = primesSlice[i]; primeSpan[computeCount] = prime; positionSpan[computeCount] = i; @@ -497,32 +499,32 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv if (computeCount > 0) { Span exponentSlice = hostSpan[..computeCount]; - ulong currentPrime = hasPreviousResidue ? previousPrime : 0UL; - bool hasDeltaSource = hasPreviousResidue; + 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; - } + 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; - } + if (cycleEnabled) + { + exponentValue %= cycle; + } - exponentSlice[i] = exponentValue; + exponentSlice[i] = exponentValue; - currentPrime = prime; - } + currentPrime = prime; + } ArrayView1D exponentView = exponentsView.SubView(0, computeCount); exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); @@ -537,6 +539,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv 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]; @@ -760,11 +763,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } - public DivisorScanSession CreateDivisorSession() - { - lock (_sync) - { - if (!_isConfigured) + public DivisorScanSession CreateDivisorSession() + { + lock (_sync) + { + if (!_isConfigured) { throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } @@ -773,15 +776,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { session.Reset(); return session; - } + } - return new DivisorScanSession(this); - } - } + return new DivisorScanSession(this); + } + } - IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() - { - return CreateDivisorSession(); - } + IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() + { + return CreateDivisorSession(); + } }