From f91f526121bc09f10e88a7346f56035ab501acc7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 21 Sep 2025 13:00:17 +0200 Subject: [PATCH] Rework by-divisor scanning order --- EvenPerfectBitScanner/Program.cs | 224 +++++++++++++++++- .../MersenneNumberDivisorGpuTesterTests.cs | 18 ++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 132 +++++++++++ 3 files changed, 363 insertions(+), 11 deletions(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index a4b1cde..08e84eb 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -1,3 +1,4 @@ +using System.Buffers; using System.Globalization; using System.Numerics; using System.Runtime.CompilerServices; @@ -11,8 +12,8 @@ namespace EvenPerfectBitScanner; internal static class Program { - private static ThreadLocal PrimeTesters = null!; - private static ThreadLocal MersenneTesters = null!; + private static ThreadLocal PrimeTesters = null!; + private static ThreadLocal MersenneTesters = null!; private static ThreadLocal PResidue = null!; // p mod d tracker (per-thread) private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; @@ -31,6 +32,7 @@ internal static class Program private static bool _useDivisor; private static bool _useResidueMode; private static bool _useByDivisorMode; + private static bool _byDivisorPrecheckOnly; private static UInt128 _divisor; private static MersenneNumberDivisorGpuTester? _divisorTester; private static MersenneNumberDivisorByDivisorGpuTester? _byDivisorTester; @@ -54,9 +56,18 @@ internal static class Program private static double _zeroFracConj = -1.0; // disabled when < 0 private static int _maxZeroConj = -1; // disabled when < 0 - private static unsafe void Main(string[] args) - { - ulong currentP = InitialP; + private struct ByDivisorPrimeState + { + internal ulong Prime; + internal ulong AllowedMax; + internal bool Completed; + internal bool Composite; + internal bool DetailedCheck; + } + + private static unsafe void Main(string[] args) + { + ulong currentP = InitialP; int threadCount = Environment.ProcessorCount; int blockSize = 1; int gpuPrimeThreads = 1; @@ -577,6 +588,7 @@ internal static class Program bool useFilter = !string.IsNullOrEmpty(filterFile); HashSet filter = []; + List filterList = []; ulong maxP = 0UL; if (useFilter) { @@ -595,7 +607,7 @@ internal static class Program if (count == 1024) { - filter.AddRange(localFilter[..count]); + filter.AddRange(localFilter[..count]); count = 0; Console.WriteLine($"Added {p}"); } @@ -604,10 +616,9 @@ internal static class Program if (count > 0) { - filter.AddRange(localFilter[..count]); + filter.AddRange(localFilter[..count]); } - // Restore this if you want to use List instead of HashSet - // filter.Sort(); + filterList = [..filter]; } if (useByDivisor) @@ -633,6 +644,14 @@ internal static class Program // Configure batch size for GPU primality sieve PrimeTester.GpuBatchSize = gpuPrimeBatch; + if (_useByDivisorMode) + { + RunByDivisorMode(filterList, divisorCyclesSearchLimit); + FlushBuffer(); + StringBuilderPool.Return(_outputBuilder!); + return; + } + if (!useDivisor && !useByDivisor) { Console.WriteLine("Warming up orders..."); @@ -713,8 +732,184 @@ internal static class Program Task.WaitAll(tasks); FlushBuffer(); - StringBuilderPool.Return(_outputBuilder!); - } + StringBuilderPool.Return(_outputBuilder!); + } + + private static void RunByDivisorMode(List primes, ulong divisorCyclesSearchLimit) + { + if (primes.Count == 0) + { + return; + } + + primes.Sort(); + List states = new(primes.Count); + + _byDivisorPrecheckOnly = true; + foreach (ulong prime in primes) + { + bool passed = IsEvenPerfectCandidate(prime, divisorCyclesSearchLimit, out bool searched, out bool detailed); + if (!passed) + { + PrintResult(prime, searched, detailed, false); + continue; + } + + ulong allowedMax = _byDivisorTester!.GetAllowedMaxDivisor(prime); + if (allowedMax < 3UL) + { + PrintResult(prime, searchedMersenne: true, detailedCheck: true, passedAllTests: true); + continue; + } + + states.Add(new ByDivisorPrimeState + { + Prime = prime, + AllowedMax = allowedMax, + Completed = false, + Composite = false, + DetailedCheck = false, + }); + } + + _byDivisorPrecheckOnly = false; + + if (states.Count == 0) + { + return; + } + + using var session = _byDivisorTester!.CreateDivisorSession(); + byte[] hitsBuffer = ArrayPool.Shared.Rent(states.Count); + ulong[] primeBuffer = ArrayPool.Shared.Rent(states.Count); + int[] indexBuffer = ArrayPool.Shared.Rent(states.Count); + + try + { + ulong divisor = 3UL; + ulong divisorLimit = _byDivisorTester!.DivisorLimit; + + while (true) + { + bool anyPending = false; + int activeCount = 0; + + for (int i = 0; i < states.Count; i++) + { + var state = states[i]; + if (state.Completed) + { + continue; + } + + if (state.AllowedMax < divisor) + { + state.Completed = true; + state.DetailedCheck = true; + states[i] = state; + PrintResult(state.Prime, searchedMersenne: true, detailedCheck: true, passedAllTests: true); + continue; + } + + anyPending = true; + + if (divisor > divisorLimit) + { + states[i] = state; + continue; + } + + primeBuffer[activeCount] = state.Prime; + indexBuffer[activeCount] = i; + activeCount++; + states[i] = state; + } + + if (!anyPending || divisor > divisorLimit) + { + break; + } + + if (activeCount == 0) + { + ulong nextDivisor = divisor + 2UL; + if (nextDivisor <= divisor) + { + break; + } + + divisor = nextDivisor; + continue; + } + + session.CheckDivisor(divisor, primeBuffer.AsSpan(0, activeCount), hitsBuffer.AsSpan(0, activeCount)); + + for (int i = 0; i < activeCount; i++) + { + if (hitsBuffer[i] == 0) + { + continue; + } + + int index = indexBuffer[i]; + var state = states[index]; + if (state.Completed) + { + continue; + } + + state.Completed = true; + state.Composite = true; + state.DetailedCheck = true; + states[index] = state; + PrintResult(state.Prime, searchedMersenne: true, detailedCheck: true, passedAllTests: false); + } + + bool remaining = false; + for (int i = 0; i < states.Count; i++) + { + if (!states[i].Completed) + { + remaining = true; + break; + } + } + + if (!remaining) + { + break; + } + + ulong next = divisor + 2UL; + if (next <= divisor) + { + break; + } + + divisor = next; + } + + for (int i = 0; i < states.Count; i++) + { + var state = states[i]; + if (state.Completed) + { + continue; + } + + state.Completed = true; + state.DetailedCheck = divisor > state.AllowedMax; + states[i] = state; + PrintResult(state.Prime, searchedMersenne: true, detailedCheck: state.DetailedCheck, passedAllTests: true); + } + } + finally + { + ArrayPool.Shared.Return(indexBuffer, clearArray: true); + ArrayPool.Shared.Return(primeBuffer, clearArray: true); + ArrayPool.Shared.Return(hitsBuffer, clearArray: true); + } + } private static unsafe void LoadResultsFile(string resultsFileName, Action lineProcessorAction) { @@ -1163,6 +1358,13 @@ internal static class Program searchedMersenne = true; if (_useByDivisorMode) { + if (_byDivisorPrecheckOnly) + { + searchedMersenne = false; + detailedCheck = false; + return true; + } + return _byDivisorTester!.IsPrime(p, out detailedCheck); } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index e41c4d4..0d767ee 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -103,5 +103,23 @@ public class MersenneNumberDivisorGpuTesterTests divisorsExhausted.Should().BeTrue(); } + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_session_checks_divisors_across_primes() + { + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + tester.ConfigureFromMaxPrime(13UL); + + using var session = tester.CreateDivisorSession(); + ulong[] primes = { 5UL, 7UL, 11UL, 13UL }; + byte[] hits = new byte[primes.Length]; + + session.CheckDivisor(23UL, primes, hits); + hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); + + Array.Fill(hits, (byte)0); + session.CheckDivisor(31UL, primes, hits); + hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 }); + } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 7201cd0..f2464b5 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -17,10 +17,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester 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; @@ -37,6 +41,48 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester } } + 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) @@ -179,6 +225,79 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester 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) @@ -234,6 +353,19 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester 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)