diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 08e84eb..36dc115 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -16,8 +16,6 @@ internal static class Program private static ThreadLocal MersenneTesters = null!; private static ThreadLocal PResidue = null!; // p mod d tracker (per-thread) private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; - - private const int ConsoleInterval = 100_000; private const int WriteBatchSize = 100; private const string DefaultCyclesPath = "divisor_cycles.bin"; private static string ResultsFileName = "even_perfect_bit_scan_results.csv"; @@ -1082,7 +1080,7 @@ internal static class Program lock (Sync) { - if (_consoleCounter >= ConsoleInterval) + if (_consoleCounter >= PerfectNumberConstants.ConsoleInterval) { printToConsole = true; primeFlag = passedAllTests && !lastWasComposite && _primeFoundAfterInit; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index f2464b5..d1c6410 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -10,422 +10,288 @@ namespace PerfectNumbers.Core.Gpu; public sealed class MersenneNumberDivisorByDivisorGpuTester { - private int _gpuBatchSize = GpuConstants.ScanBatchSize; - private readonly object _sync = new(); - private ulong _divisorLimit; - private bool _isConfigured; - private ulong _lastStatusDivisor; - - private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelByPrimeCache = new(); - - private Action, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); - - private Action, ArrayView> GetKernelByPrime(Accelerator accelerator) => - _kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernelByPrime)); - - public int GpuBatchSize - { - get => _gpuBatchSize; - set => _gpuBatchSize = 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); - } - } - - public DivisorScanSession CreateDivisorSession() - { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - return new DivisorScanSession(this); - } - } - - public bool IsPrime(ulong prime, out bool divisorsExhausted) - { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - ulong allowedMax = ComputeAllowedMaxDivisor(prime); - if (allowedMax < 3UL) - { - divisorsExhausted = true; - return true; - } - - bool composite = CheckDivisors(prime, allowedMax, out ulong lastProcessed, out bool coveredRange); - - if (!composite && lastProcessed != 0UL) - { - ReportStatus(lastProcessed); - } - - if (composite) - { - divisorsExhausted = true; - return false; - } - - divisorsExhausted = coveredRange; - return true; - } - } - - private bool CheckDivisors(ulong prime, ulong allowedMax, out ulong lastProcessed, out bool coveredRange) - { - allowedMax = Math.Min(allowedMax, _divisorLimit); - lastProcessed = 0UL; - - if (allowedMax < 3UL) - { - coveredRange = true; - return false; - } - - var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); - var accelerator = gpuLease.Accelerator; - var kernel = GetKernel(accelerator); - - int batchCapacity = _gpuBatchSize; - var divisorsBuffer = accelerator.Allocate1D(batchCapacity); - var hitsBuffer = accelerator.Allocate1D(batchCapacity); - - ulong[] divisors = ArrayPool.Shared.Rent(batchCapacity); - byte[] hits = ArrayPool.Shared.Rent(batchCapacity); - - bool composite = false; - bool processedAll = false; - ulong divisor = 3UL; - - try - { - while (divisor <= allowedMax) - { - int batchSize = 0; - bool reachedEndInBatch = false; - - while (batchSize < batchCapacity && divisor <= allowedMax) - { - ulong current = divisor; - divisors[batchSize++] = current; - - ulong next = current + 2UL; - if (next <= current) - { - reachedEndInBatch = true; - break; - } - - if (next > allowedMax) - { - divisor = next; - reachedEndInBatch = true; - break; - } - - divisor = next; - } - - if (batchSize == 0) - { - break; - } - - var divisorView = divisorsBuffer.View.SubView(0, batchSize); - var hitsView = hitsBuffer.View.SubView(0, batchSize); - - Span divisorSpan = divisors.AsSpan(0, batchSize); - Span hitsSpan = hits.AsSpan(0, batchSize); - - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorSpan), batchSize); - kernel(batchSize, prime, divisorView, hitsView); - accelerator.Synchronize(); - - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); - - for (int i = 0; i < batchSize; i++) - { - lastProcessed = divisorSpan[i]; - if (hitsSpan[i] != 0) - { - composite = true; - break; - } - } - - if (composite) - { - break; - } - - if (reachedEndInBatch) - { - processedAll = true; - break; - } - } - - coveredRange = composite || processedAll || divisor > allowedMax; - } - finally - { - ArrayPool.Shared.Return(hits, clearArray: true); - ArrayPool.Shared.Return(divisors, clearArray: true); - hitsBuffer.Dispose(); - divisorsBuffer.Dispose(); - gpuLease.Dispose(); - } - - return composite; - } - - public sealed class DivisorScanSession : IDisposable - { - private readonly MersenneNumberDivisorByDivisorGpuTester _owner; - private readonly GpuContextPool.GpuContextLease _lease; - private readonly Accelerator _accelerator; - private readonly Action, ArrayView> _kernel; - private readonly MemoryBuffer1D _primesBuffer; - private readonly MemoryBuffer1D _hitsBuffer; - private readonly ulong[] _primesHost; - private readonly byte[] _hitsHost; - private bool _disposed; - - internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) - { - _owner = owner; - _lease = GpuContextPool.RentPreferred(preferCpu: false); - _accelerator = _lease.Accelerator; - _kernel = owner.GetKernelByPrime(_accelerator); - _primesBuffer = _accelerator.Allocate1D(owner._gpuBatchSize); - _hitsBuffer = _accelerator.Allocate1D(owner._gpuBatchSize); - _primesHost = ArrayPool.Shared.Rent(owner._gpuBatchSize); - _hitsHost = ArrayPool.Shared.Rent(owner._gpuBatchSize); - } - - public void CheckDivisor(ulong divisor, ReadOnlySpan primes, Span hits) - { - if (_disposed) - { - throw new ObjectDisposedException(nameof(DivisorScanSession)); - } - - if (primes.Length == 0) - { - return; - } - - int offset = 0; - while (offset < primes.Length) - { - int batchSize = Math.Min(_owner._gpuBatchSize, primes.Length - offset); - Span primeSpan = _primesHost.AsSpan(0, batchSize); - primes.Slice(offset, batchSize).CopyTo(primeSpan); - - var primesView = _primesBuffer.View.SubView(0, batchSize); - var hitsView = _hitsBuffer.View.SubView(0, batchSize); - - primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primeSpan), batchSize); - _kernel(batchSize, divisor, primesView, hitsView); - _accelerator.Synchronize(); - - Span hitSpan = _hitsHost.AsSpan(0, batchSize); - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitSpan), batchSize); - hitSpan.CopyTo(hits.Slice(offset, batchSize)); - offset += batchSize; - } - } - - public void Dispose() - { - if (_disposed) - { - return; - } - - _disposed = true; - ArrayPool.Shared.Return(_primesHost, clearArray: true); - ArrayPool.Shared.Return(_hitsHost, clearArray: true); - _primesBuffer.Dispose(); - _hitsBuffer.Dispose(); - _lease.Dispose(); - } - } - - private void ReportStatus(ulong divisor) - { - if (divisor <= _lastStatusDivisor) - { - return; - } - - _lastStatusDivisor = divisor; - Console.WriteLine($"...processed by-divisor candidate = {divisor}"); - } - - 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 ulong ComputeAllowedMaxDivisor(ulong prime) - { - if (prime <= 1UL) - { - return 0UL; - } - - if (prime - 1UL >= 64UL) - { - return _divisorLimit; - } - - ulong maxByPrime = (1UL << (int)(prime - 1UL)) - 1UL; - return Math.Min(maxByPrime, _divisorLimit); - } - - private static void CheckKernel(Index1D index, ulong prime, ArrayView divisors, ArrayView hits) - { - ulong divisor = divisors[index]; - if (divisor <= 1UL || (divisor & 1UL) == 0UL) - { - hits[index] = 0; - return; - } - - ulong pow = Pow2Mod(prime, divisor); - hits[index] = pow == 1UL ? (byte)1 : (byte)0; - } - - private static void CheckKernelByPrime(Index1D index, ulong divisor, ArrayView primes, ArrayView hits) - { - ulong prime = primes[index]; - if (divisor <= 1UL || (divisor & 1UL) == 0UL) - { - hits[index] = 0; - return; - } - - ulong pow = Pow2Mod(prime, divisor); - hits[index] = pow == 1UL ? (byte)1 : (byte)0; - } - - private static ulong Pow2Mod(ulong exponent, ulong modulus) - { - if (modulus <= 1UL) - { - return 0UL; - } - - ulong result = 1UL % modulus; - ulong baseVal = 2UL % modulus; - ulong exp = exponent; - while (exp > 0UL) - { - if ((exp & 1UL) != 0UL) - { - result = MulMod(result, baseVal, modulus); - } - - exp >>= 1; - if (exp == 0UL) - { - break; - } - - baseVal = MulMod(baseVal, baseVal, modulus); - } - - return result; - } - - private static ulong MulMod(ulong a, ulong b, ulong modulus) - { - if (modulus == 0UL) - { - return 0UL; - } - - ulong result = 0UL; - ulong x = a % modulus; - ulong y = b; - while (y > 0UL) - { - if ((y & 1UL) != 0UL) - { - result += x; - if (result >= modulus) - { - result -= modulus; - } - } - - x <<= 1; - if (x >= modulus) - { - x -= modulus; - } - - y >>= 1; - } - - return result; - } + private int _gpuBatchSize = GpuConstants.ScanBatchSize; + private readonly object _sync = new(); + private ulong _divisorLimit; + private bool _isConfigured; + private ulong _lastStatusDivisor; + + private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); + + private Action, ArrayView> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); + + public int GpuBatchSize + { + get => _gpuBatchSize; + set => _gpuBatchSize = Math.Max(1, value); + } + + public void ConfigureFromMaxPrime(ulong maxPrime) + { + lock (_sync) + { + _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); + _lastStatusDivisor = 0UL; + _isConfigured = true; + } + } + + public bool IsPrime(ulong prime, out bool divisorsExhausted) + { + lock (_sync) + { + if (!_isConfigured) + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } + + ulong allowedMax = ComputeAllowedMaxDivisor(prime); + if (allowedMax < 3UL) + { + divisorsExhausted = true; + return true; + } + + bool composite = CheckDivisors(prime, allowedMax, out ulong lastProcessed, out bool coveredRange); + + if (lastProcessed != 0UL) + { + ReportStatus(lastProcessed); + } + + if (composite) + { + divisorsExhausted = true; + return false; + } + + divisorsExhausted = coveredRange; + return true; + } + } + + private bool CheckDivisors(ulong prime, ulong allowedMax, out ulong lastProcessed, out bool coveredRange) + { + allowedMax = Math.Min(allowedMax, _divisorLimit); + lastProcessed = 0UL; + + if (allowedMax < 3UL) + { + coveredRange = true; + return false; + } + + var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); + var accelerator = gpuLease.Accelerator; + var kernel = GetKernel(accelerator); + + int batchCapacity = _gpuBatchSize; + var divisorsBuffer = accelerator.Allocate1D(batchCapacity); + var hitsBuffer = accelerator.Allocate1D(batchCapacity); + + ulong[] divisors = ArrayPool.Shared.Rent(batchCapacity); + byte[] hits = ArrayPool.Shared.Rent(batchCapacity); + + bool composite = false; + bool processedAll = false; + ulong divisor = 3UL; + + try + { + while (divisor <= allowedMax) + { + int batchSize = 0; + bool reachedEndInBatch = false; + + while (batchSize < batchCapacity && divisor <= allowedMax) + { + ulong current = divisor; + divisors[batchSize++] = current; + + ulong next = current + 2UL; + if (next <= current) + { + reachedEndInBatch = true; + break; + } + + if (next > allowedMax) + { + divisor = next; + reachedEndInBatch = true; + break; + } + + divisor = next; + } + + if (batchSize == 0) + { + break; + } + + var divisorView = divisorsBuffer.View.SubView(0, batchSize); + var hitsView = hitsBuffer.View.SubView(0, batchSize); + + Span divisorSpan = divisors.AsSpan(0, batchSize); + Span hitsSpan = hits.AsSpan(0, batchSize); + + divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorSpan), batchSize); + kernel(batchSize, prime, divisorView, hitsView); + accelerator.Synchronize(); + + hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); + + for (int i = 0; i < batchSize; i++) + { + lastProcessed = divisorSpan[i]; + if (hitsSpan[i] != 0) + { + composite = true; + break; + } + } + + if (composite) + { + break; + } + + if (reachedEndInBatch) + { + processedAll = true; + break; + } + } + + coveredRange = composite || processedAll || divisor > allowedMax; + } + finally + { + ArrayPool.Shared.Return(hits, clearArray: true); + ArrayPool.Shared.Return(divisors, clearArray: true); + hitsBuffer.Dispose(); + divisorsBuffer.Dispose(); + gpuLease.Dispose(); + } + + return composite; + } + + private void ReportStatus(ulong divisor) + { + if (++_lastStatusDivisor == PerfectNumberConstants.ConsoleInterval) + { + Console.WriteLine($"...processed by-divisor candidate = {divisor}"); + _lastStatusDivisor = 0; + } + } + + 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 ulong ComputeAllowedMaxDivisor(ulong prime) + { + if (prime <= 1UL) + { + return 0UL; + } + + if (prime - 1UL >= 64UL) + { + return _divisorLimit; + } + + ulong maxByPrime = (1UL << (int)(prime - 1UL)) - 1UL; + return Math.Min(maxByPrime, _divisorLimit); + } + + private static void CheckKernel(Index1D index, ulong prime, ArrayView divisors, ArrayView hits) + { + ulong divisor = divisors[index]; + if (divisor <= 1UL || (divisor & 1UL) == 0UL) + { + hits[index] = 0; + return; + } + + ulong pow = Pow2Mod(prime, divisor); + hits[index] = pow == 1UL ? (byte)1 : (byte)0; + } + + private static ulong Pow2Mod(ulong exponent, ulong modulus) + { + if (modulus <= 1UL) + { + return 0UL; + } + + ulong result = 1UL % modulus; + ulong baseVal = 2UL % modulus; + ulong exp = exponent; + while (exp > 0UL) + { + if ((exp & 1UL) != 0UL) + { + result = MulMod(result, baseVal, modulus); + } + + exp >>= 1; + if (exp == 0UL) + { + break; + } + + baseVal = MulMod(baseVal, baseVal, modulus); + } + + return result; + } + + private static ulong MulMod(ulong a, ulong b, ulong modulus) + { + if (modulus == 0UL) + { + return 0UL; + } + + ulong result = 0UL; + ulong x = a % modulus; + ulong y = b; + while (y > 0UL) + { + if ((y & 1UL) != 0UL) + { + result += x; + if (result >= modulus) + { + result -= modulus; + } + } + + x <<= 1; + if (x >= modulus) + { + x -= modulus; + } + + y >>= 1; + } + + return result; + } } diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index eb53ca1..f28077c 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -3,8 +3,9 @@ namespace PerfectNumbers.Core; public static class PerfectNumberConstants { public const ulong BiggestKnownEvenPerfectP = 136279841UL; - public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; - public const int MaxQForDivisorCycles = 4_000_000; + public const int ConsoleInterval = 100_000; public const ulong ExtraDivisorCycleSearchLimit = 64UL; + public const int MaxQForDivisorCycles = 4_000_000; + public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; }