diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index a4b1cde..cd4d37f 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -15,8 +15,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"; @@ -887,7 +885,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 7201cd0..d1c6410 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -10,290 +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 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>> _kernelCache = new(); - private Action, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); + 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 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 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."); - } + 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; - } + ulong allowedMax = ComputeAllowedMaxDivisor(prime); + if (allowedMax < 3UL) + { + divisorsExhausted = true; + return true; + } - bool composite = CheckDivisors(prime, allowedMax, out ulong lastProcessed, out bool coveredRange); + bool composite = CheckDivisors(prime, allowedMax, out ulong lastProcessed, out bool coveredRange); - if (!composite && lastProcessed != 0UL) - { - ReportStatus(lastProcessed); - } + if (lastProcessed != 0UL) + { + ReportStatus(lastProcessed); + } - if (composite) - { - divisorsExhausted = true; - return false; - } + if (composite) + { + divisorsExhausted = true; + return false; + } - divisorsExhausted = coveredRange; - return true; - } - } + divisorsExhausted = coveredRange; + return true; + } + } - private bool CheckDivisors(ulong prime, ulong allowedMax, out ulong lastProcessed, out bool coveredRange) - { - allowedMax = Math.Min(allowedMax, _divisorLimit); - lastProcessed = 0UL; + 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; - } + if (allowedMax < 3UL) + { + coveredRange = true; + return false; + } - var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); - var accelerator = gpuLease.Accelerator; - var kernel = GetKernel(accelerator); + 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); + 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); + ulong[] divisors = ArrayPool.Shared.Rent(batchCapacity); + byte[] hits = ArrayPool.Shared.Rent(batchCapacity); - bool composite = false; - bool processedAll = false; - ulong divisor = 3UL; + bool composite = false; + bool processedAll = false; + ulong divisor = 3UL; - try - { - while (divisor <= allowedMax) - { - int batchSize = 0; - bool reachedEndInBatch = false; + try + { + while (divisor <= allowedMax) + { + int batchSize = 0; + bool reachedEndInBatch = false; - while (batchSize < batchCapacity && divisor <= allowedMax) - { - ulong current = divisor; - divisors[batchSize++] = current; + while (batchSize < batchCapacity && divisor <= allowedMax) + { + ulong current = divisor; + divisors[batchSize++] = current; - ulong next = current + 2UL; - if (next <= current) - { - reachedEndInBatch = true; - break; - } + ulong next = current + 2UL; + if (next <= current) + { + reachedEndInBatch = true; + break; + } - if (next > allowedMax) - { - divisor = next; - reachedEndInBatch = true; - break; - } + if (next > allowedMax) + { + divisor = next; + reachedEndInBatch = true; + break; + } - divisor = next; - } + divisor = next; + } - if (batchSize == 0) - { - break; - } + if (batchSize == 0) + { + break; + } - var divisorView = divisorsBuffer.View.SubView(0, batchSize); - var hitsView = hitsBuffer.View.SubView(0, batchSize); + 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); + 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(); + divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorSpan), batchSize); + kernel(batchSize, prime, divisorView, hitsView); + accelerator.Synchronize(); - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); + hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); - for (int i = 0; i < batchSize; i++) - { - lastProcessed = divisorSpan[i]; - if (hitsSpan[i] != 0) - { - composite = true; - break; - } - } + for (int i = 0; i < batchSize; i++) + { + lastProcessed = divisorSpan[i]; + if (hitsSpan[i] != 0) + { + composite = true; + break; + } + } - if (composite) - { - break; - } + if (composite) + { + break; + } - if (reachedEndInBatch) - { - processedAll = true; - 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(); - } + 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; - } + return composite; + } - private void ReportStatus(ulong divisor) - { - if (divisor <= _lastStatusDivisor) - { - return; - } + private void ReportStatus(ulong divisor) + { + if (++_lastStatusDivisor == PerfectNumberConstants.ConsoleInterval) + { + Console.WriteLine($"...processed by-divisor candidate = {divisor}"); + _lastStatusDivisor = 0; + } + } - _lastStatusDivisor = divisor; - Console.WriteLine($"...processed by-divisor candidate = {divisor}"); - } + private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) + { + if (maxPrime <= 1UL) + { + return 0UL; + } - private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) - { - if (maxPrime <= 1UL) - { - return 0UL; - } + if (maxPrime - 1UL >= 64UL) + { + return ulong.MaxValue; + } - if (maxPrime - 1UL >= 64UL) - { - return ulong.MaxValue; - } + return (1UL << (int)(maxPrime - 1UL)) - 1UL; + } - return (1UL << (int)(maxPrime - 1UL)) - 1UL; - } + private ulong ComputeAllowedMaxDivisor(ulong prime) + { + if (prime <= 1UL) + { + return 0UL; + } - private ulong ComputeAllowedMaxDivisor(ulong prime) - { - if (prime <= 1UL) - { - return 0UL; - } + if (prime - 1UL >= 64UL) + { + return _divisorLimit; + } - if (prime - 1UL >= 64UL) - { - return _divisorLimit; - } + ulong maxByPrime = (1UL << (int)(prime - 1UL)) - 1UL; + return Math.Min(maxByPrime, _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; + } - 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; + } - 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; + } - 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); + } - 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; + } - exp >>= 1; - if (exp == 0UL) - { - break; - } + baseVal = MulMod(baseVal, baseVal, modulus); + } - baseVal = MulMod(baseVal, baseVal, modulus); - } + return result; + } - return result; - } + private static ulong MulMod(ulong a, ulong b, ulong modulus) + { + if (modulus == 0UL) + { + return 0UL; + } - 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; + } + } - 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; + } - x <<= 1; - if (x >= modulus) - { - x -= modulus; - } + y >>= 1; + } - y >>= 1; - } - - return result; - } + 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; }