diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 949d5ae..195f355 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -6,6 +6,7 @@ using System.Text; using Open.Collections; using Open.Numeric.Primes; using PerfectNumbers.Core; +using PerfectNumbers.Core.Cpu; using PerfectNumbers.Core.Gpu; namespace EvenPerfectBitScanner; @@ -33,8 +34,8 @@ internal static class Program private static bool _useByDivisorMode; private static bool _byDivisorPrecheckOnly; private static UInt128 _divisor; - private static MersenneNumberDivisorGpuTester? _divisorTester; - private static MersenneNumberDivisorByDivisorGpuTester? _byDivisorTester; + private static MersenneNumberDivisorGpuTester? _divisorTester; + private static IMersenneNumberDivisorByDivisorTester? _byDivisorTester; private static ulong? _orderWarmupLimitOverride; private static unsafe delegate* _transformP; private static long _state; @@ -518,14 +519,14 @@ internal static class Program _divisorTester = new MersenneNumberDivisorGpuTester(); } - else if (useByDivisor) - { - _byDivisorTester = new MersenneNumberDivisorByDivisorGpuTester - { - GpuBatchSize = scanBatchSize, - UseDivisorCycles = useDivisorCycles, - }; - } + else if (useByDivisor) + { + _byDivisorTester = mersenneOnGpu + ? new MersenneNumberDivisorByDivisorGpuTester() + : new MersenneNumberDivisorByDivisorCpuTester(); + _byDivisorTester.BatchSize = scanBatchSize; + _byDivisorTester.UseDivisorCycles = useDivisorCycles; + } // Load RLE blacklist (optional) if (!string.IsNullOrEmpty(_rleBlacklistPath)) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs new file mode 100644 index 0000000..eaf7f71 --- /dev/null +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -0,0 +1,48 @@ +using System; +using FluentAssertions; +using PerfectNumbers.Core.Cpu; +using Xunit; + +namespace PerfectNumbers.Core.Tests; + +public class MersenneNumberDivisorCpuTesterTests +{ + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_tester_tracks_divisors_across_primes() + { + var tester = new MersenneNumberDivisorByDivisorCpuTester(); + tester.ConfigureFromMaxPrime(11UL); + + tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue(); + divisorsExhausted.Should().BeTrue(); + + tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue(); + divisorsExhausted.Should().BeTrue(); + + tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse(); + divisorsExhausted.Should().BeTrue(); + } + + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_session_checks_divisors_across_primes() + { + var tester = new MersenneNumberDivisorByDivisorCpuTester(); + tester.ConfigureFromMaxPrime(13UL); + + using var session = tester.CreateDivisorSession(); + ulong[] primes = { 5UL, 7UL, 11UL, 13UL }; + byte[] hits = new byte[primes.Length]; + + ulong cycle23 = MersenneDivisorCycles.CalculateCycleLength(23UL); + session.CheckDivisor(23UL, cycle23, primes, hits); + hits.Should().ContainInOrder(new byte[] { 0, 0, 1, 0 }); + + Array.Fill(hits, (byte)0); + ulong cycle31 = MersenneDivisorCycles.CalculateCycleLength(31UL); + session.CheckDivisor(31UL, cycle31, primes, hits); + hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 }); + } +} + diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs new file mode 100644 index 0000000..e9b2cea --- /dev/null +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -0,0 +1,562 @@ +using System; +using System.Buffers; +using System.Collections.Concurrent; +using System.Runtime.CompilerServices; +using PerfectNumbers.Core; +using PerfectNumbers.Core.Gpu; +using MontgomeryDivisorData = PerfectNumbers.Core.Gpu.MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData; + +namespace PerfectNumbers.Core.Cpu; + +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); + } + } +} + diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 45480a1..b34d441 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -9,7 +9,7 @@ using PerfectNumbers.Core; namespace PerfectNumbers.Core.Gpu; -public sealed class MersenneNumberDivisorByDivisorGpuTester +public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDivisorByDivisorTester { private int _gpuBatchSize = GpuConstants.ScanBatchSize; private readonly object _sync = new(); @@ -39,11 +39,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester 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; + } public bool UseDivisorCycles { @@ -348,8 +354,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; } - public sealed class DivisorScanSession : IDisposable - { + public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession + { private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; @@ -754,11 +760,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester } } - 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."); } @@ -767,10 +773,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester { session.Reset(); return session; - } + } - return new DivisorScanSession(this); - } - } + return new DivisorScanSession(this); + } + } + + IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() + { + return CreateDivisorSession(); + } } diff --git a/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs new file mode 100644 index 0000000..93413b1 --- /dev/null +++ b/PerfectNumbers.Core/IMersenneNumberDivisorByDivisorTester.cs @@ -0,0 +1,26 @@ +using System; + +namespace PerfectNumbers.Core; + +public interface IMersenneNumberDivisorByDivisorTester +{ + bool UseDivisorCycles { get; set; } + + int BatchSize { get; set; } + + void ConfigureFromMaxPrime(ulong maxPrime); + + ulong DivisorLimit { get; } + + ulong GetAllowedMaxDivisor(ulong prime); + + bool IsPrime(ulong prime, out bool divisorsExhausted); + + IDivisorScanSession CreateDivisorSession(); + + public interface IDivisorScanSession : IDisposable + { + void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits); + } +} +