From 9b9c6959bec06fabd930d0e255233e5299b6eb8f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 25 Sep 2025 18:51:47 +0200 Subject: [PATCH] Reuse cycle remainder stepper in by-divisor testers --- .../Pow2MontgomeryModBenchmarks.cs | 20 +- EvenPerfectBitScanner/Program.cs | 42 +- .../MersenneNumberDivisorGpuTesterTests.cs | 63 ++- ...MersenneNumberDivisorByDivisorCpuTester.cs | 64 +-- PerfectNumbers.Core/CycleRemainderStepper.cs | 84 ++++ PerfectNumbers.Core/DivisorCycleCache.cs | 154 ++++++- ...MersenneNumberDivisorByDivisorGpuTester.cs | 428 +++++++++--------- .../Gpu/MersenneNumberDivisorGpuTester.cs | 52 +-- PerfectNumbers.Core/MontgomeryDivisorData.cs | 148 ++++++ PerfectNumbers.Core/ULongExtensions.cs | 14 +- 10 files changed, 738 insertions(+), 331 deletions(-) create mode 100644 PerfectNumbers.Core/CycleRemainderStepper.cs create mode 100644 PerfectNumbers.Core/MontgomeryDivisorData.cs diff --git a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs index 299874f..3669e68 100644 --- a/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/Pow2MontgomeryModBenchmarks.cs @@ -1,9 +1,9 @@ using System; using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; +using System.Numerics; using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; -using System.Numerics; namespace EvenPerfectBitScanner.Benchmarks; @@ -13,8 +13,8 @@ public class Pow2MontgomeryModBenchmarks { private const int SampleCount = 256; - private readonly MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[] _smallDivisors = new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[SampleCount]; - private readonly MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[] _largeDivisors = new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[SampleCount]; + private readonly MontgomeryDivisorData[] _smallDivisors = new MontgomeryDivisorData[SampleCount]; + private readonly MontgomeryDivisorData[] _largeDivisors = new MontgomeryDivisorData[SampleCount]; private readonly ulong[] _smallExponents = new ulong[SampleCount]; private readonly ulong[] _largeExponents = new ulong[SampleCount]; private readonly ulong[] _smallCycles = new ulong[SampleCount]; @@ -184,7 +184,7 @@ public class Pow2MontgomeryModBenchmarks return checksum; } - private void GetData(out ulong[] exponents, out MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData[] divisors, out ulong[] cycles) + private void GetData(out ulong[] exponents, out MontgomeryDivisorData[] divisors, out ulong[] cycles) { if (Scale == InputScale.Small) { @@ -217,14 +217,14 @@ public class Pow2MontgomeryModBenchmarks return (modulus, (ulong)bitLength); } - private static MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) + private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) { if (modulus <= 1UL || (modulus & 1UL) == 0UL) { - return new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); + return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); } - return new MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData( + return new MontgomeryDivisorData( modulus, ComputeMontgomeryNPrime(modulus), ComputeMontgomeryResidue(1UL, modulus), @@ -245,7 +245,7 @@ public class Pow2MontgomeryModBenchmarks return unchecked(0UL - inv); } - private static ulong Pow2MontgomeryModLeftToRight(ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor) + private static ulong Pow2MontgomeryModLeftToRight(ulong exponent, in MontgomeryDivisorData divisor) { ulong modulus = divisor.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) @@ -277,7 +277,7 @@ public class Pow2MontgomeryModBenchmarks return result.MontgomeryMultiply(1UL, modulus, nPrime); } - private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor) + private static ulong Pow2MontgomeryModBatched4(ulong exponent, in MontgomeryDivisorData divisor) { ulong modulus = divisor.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) @@ -326,7 +326,7 @@ public class Pow2MontgomeryModBenchmarks return result.MontgomeryMultiply(1UL, modulus, nPrime); } - private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor) + private static ulong Pow2MontgomeryModSlidingWindow(ulong exponent, in MontgomeryDivisorData divisor) { ulong modulus = divisor.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index fada2ac..60b3a3d 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -899,21 +899,29 @@ internal static class Program divisor = AcquireNextDivisor(ref nextDivisor, divisorLimit, ref divisorsExhaustedFlag, ref finalDivisorBits, out exhausted, ref localDivisorCursor, ref localDivisorsRemaining); if (divisor != 0UL) { - divisorCycle = useDivisorCycles - ? DivisorCycleCache.Shared.Acquire(divisor).GetCycle(divisor) - : 0UL; + DivisorCycleCache.CycleBlock? cycleBlock = null; + if (useDivisorCycles) + { + cycleBlock = DivisorCycleCache.Shared.Acquire(divisor); + divisorCycle = cycleBlock.GetCycle(divisor); + } + else + { + divisorCycle = 0UL; + } - activeCount = BuildPrimeBuffer(divisor, primeValues, allowedMaxValues, stateFlags, primeBuffer, indexBuffer, completionsBuffer, ref completionsCount, ref remainingStates, activeStateMask, ref activeStartIndex); + activeCount = BuildPrimeBuffer(divisor, primeValues, allowedMaxValues, stateFlags, primeBuffer, indexBuffer, completionsBuffer, ref completionsCount, ref remainingStates, activeStateMask, ref activeStartIndex); - if (completionsCount > 0) - { - FlushPendingResults(completionsBuffer, ref completionsCount); - } + if (completionsCount > 0) + { + FlushPendingResults(completionsBuffer, ref completionsCount); + } - if (activeCount == 0) - { - continue; - } + if (activeCount == 0) + { + cycleBlock?.Dispose(); + continue; + } hitsSpan = hitsBuffer.AsSpan(0, activeCount); hitsSpan.Clear(); @@ -940,10 +948,12 @@ internal static class Program } } - if (compositesCount > 0) - { - FlushPendingResults(compositesBuffer, ref compositesCount); - } + if (compositesCount > 0) + { + FlushPendingResults(compositesBuffer, ref compositesCount); + } + + cycleBlock?.Dispose(); continue; } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index ed8e6f9..9317c02 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -1,6 +1,7 @@ using FluentAssertions; using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; +using PerfectNumbers.Core.Cpu; using System.Reflection; using Xunit; @@ -37,7 +38,7 @@ public class MersenneNumberDivisorGpuTesterTests var tester = new MersenneNumberDivisorGpuTester(); typeof(MersenneNumberDivisorGpuTester) .GetField("_divisorCandidates", BindingFlags.NonPublic | BindingFlags.Static)! - .SetValue(null, Array.Empty<(ulong, ulong)>()); + .SetValue(null, Array.Empty<(ulong, uint)>()); tester.IsPrime(11UL, UInt128.Zero, 0UL, out bool divisorsExhausted).Should().BeTrue(); divisorsExhausted.Should().BeFalse(); @@ -102,6 +103,22 @@ public class MersenneNumberDivisorGpuTesterTests divisorsExhausted.Should().BeTrue(); } + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_gpu_tester_uses_cycle_remainders_per_divisor() + { + var tester = new MersenneNumberDivisorByDivisorGpuTester + { + UseDivisorCycles = true, + GpuBatchSize = 8, + }; + + tester.ConfigureFromMaxPrime(13UL); + + tester.IsPrime(11UL, out bool divisorsExhausted).Should().BeFalse(); + divisorsExhausted.Should().BeTrue(); + } + [Fact] [Trait("Category", "Fast")] public void ByDivisor_session_checks_divisors_across_primes() @@ -122,5 +139,49 @@ public class MersenneNumberDivisorGpuTesterTests session.CheckDivisor(31UL, cycle31, primes, hits); hits.Should().ContainInOrder(new byte[] { 1, 0, 0, 0 }); } + + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_gpu_session_reuses_cycle_remainders_across_batches() + { + var tester = new MersenneNumberDivisorByDivisorGpuTester + { + UseDivisorCycles = true, + GpuBatchSize = 2, + }; + + tester.ConfigureFromMaxPrime(17UL); + + using var session = tester.CreateDivisorSession(); + ulong[] primes = { 5UL, 7UL, 11UL, 13UL, 17UL }; + 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, 0 }); + } + + [Fact] + [Trait("Category", "Fast")] + public void ByDivisor_cpu_session_reuses_cycle_remainders_across_primes() + { + var tester = new MersenneNumberDivisorByDivisorCpuTester + { + UseDivisorCycles = true, + BatchSize = 2, + }; + + tester.ConfigureFromMaxPrime(17UL); + + using var session = tester.CreateDivisorSession(); + ulong[] primes = { 5UL, 7UL, 11UL, 13UL, 17UL }; + 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, 0 }); + } } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 376a2a7..53a6417 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,9 +1,7 @@ using System; using System.Collections.Concurrent; -using System.Runtime.CompilerServices; +using System.Numerics; using PerfectNumbers.Core; -using PerfectNumbers.Core.Gpu; -using MontgomeryDivisorData = PerfectNumbers.Core.Gpu.MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData; namespace PerfectNumbers.Core.Cpu; @@ -167,12 +165,14 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv processedCount++; lastProcessed = divisor; - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); byte hit; + ulong divisorCycle = 0UL; + DivisorCycleCache.CycleBlock? cycleBlock = null; if (useCycles) { - DivisorCycleCache.CycleBlock cycleBlock = DivisorCycleCache.Shared.Acquire(divisor); - ulong divisorCycle = cycleBlock.GetCycle(divisor); + cycleBlock = DivisorCycleCache.Shared.Acquire(divisor); + divisorCycle = cycleBlock.GetCycle(divisor); hit = CheckDivisor(prime, divisorCycle != 0UL, divisorCycle, divisorData); } else @@ -180,6 +180,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv hit = CheckDivisor(prime, false, 0UL, divisorData); } + cycleBlock?.Dispose(); + if (hit != 0) { processedAll = true; @@ -270,35 +272,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv 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; @@ -335,26 +308,31 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } bool cycleEnabled = _owner._useDivisorCycles && divisorCycle != 0UL; - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); - for (int i = 0; i < length; i++) + if (cycleEnabled) { - ulong prime = primes[i]; + var stepper = new CycleRemainderStepper(divisorCycle); - if (cycleEnabled) + for (int i = 0; i < length; i++) { - ulong remainder = prime % divisorCycle; + ulong remainder = stepper.ComputeNext(primes[i]); if (remainder != 0UL) { hits[i] = 0; continue; } + + ulong residue = remainder.Pow2MontgomeryModFromCycleRemainder(divisorData); + hits[i] = residue == 1UL ? (byte)1 : (byte)0; } - ulong residue = cycleEnabled - ? prime.Pow2MontgomeryModWithCycle(divisorCycle, divisorData) - : prime.Pow2MontgomeryMod(divisorData); + return; + } + for (int i = 0; i < length; i++) + { + ulong residue = primes[i].Pow2MontgomeryMod(divisorData); hits[i] = residue == 1UL ? (byte)1 : (byte)0; } } diff --git a/PerfectNumbers.Core/CycleRemainderStepper.cs b/PerfectNumbers.Core/CycleRemainderStepper.cs new file mode 100644 index 0000000..96f17dc --- /dev/null +++ b/PerfectNumbers.Core/CycleRemainderStepper.cs @@ -0,0 +1,84 @@ +using System.Numerics; + +namespace PerfectNumbers.Core; + +internal struct CycleRemainderStepper +{ + private readonly ulong _cycleLength; + private ulong _previousPrime; + private ulong _currentRemainder; + private bool _hasState; + + public CycleRemainderStepper(ulong cycleLength) + { + _cycleLength = cycleLength; + _previousPrime = 0UL; + _currentRemainder = 0UL; + _hasState = false; + } + + public bool HasCycle => _cycleLength != 0UL; + + public void Reset() + { + _previousPrime = 0UL; + _currentRemainder = 0UL; + _hasState = false; + } + + public void InitializeState(ulong previousPrime, ulong previousRemainder, bool hasState) + { + _previousPrime = previousPrime; + _currentRemainder = previousRemainder; + _hasState = hasState; + } + + public ulong ComputeNext(ulong prime) + { + if (_cycleLength == 0UL) + { + _previousPrime = prime; + _currentRemainder = 0UL; + _hasState = true; + return 0UL; + } + + if (!_hasState || prime <= _previousPrime) + { + _currentRemainder = prime % _cycleLength; + _hasState = true; + } + else + { + ulong delta = prime - _previousPrime; + if (ulong.MaxValue - _currentRemainder < delta) + { + UInt128 extended = (UInt128)_currentRemainder + delta; + _currentRemainder = (ulong)(extended % _cycleLength); + } + else + { + _currentRemainder += delta; + + if (_currentRemainder >= _cycleLength) + { + _currentRemainder -= _cycleLength; + if (_currentRemainder >= _cycleLength) + { + _currentRemainder %= _cycleLength; + } + } + } + } + + _previousPrime = prime; + return _currentRemainder; + } + + public void CaptureState(out ulong previousPrime, out ulong previousRemainder, out bool hasState) + { + previousPrime = _previousPrime; + previousRemainder = _currentRemainder; + hasState = _hasState; + } +} diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index d825534..d0294ce 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -3,6 +3,7 @@ using System.Buffers; using System.Collections.Concurrent; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; +using System.Threading; using System.Threading.Tasks; using ILGPU; using ILGPU.Runtime; @@ -12,16 +13,52 @@ namespace PerfectNumbers.Core; public sealed class DivisorCycleCache { - public sealed class CycleBlock(int index, ulong start, ulong[] cycles) - { - internal readonly int Index = index; - internal readonly ulong Start = start; - internal readonly ulong End = start + (ulong)cycles.Length - 1UL; - internal readonly ulong[] Cycles = cycles; + public sealed class CycleBlock : IDisposable + { + private readonly DivisorCycleCache _owner; + private int _referenceCount; - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ulong GetCycle(ulong divisor) => Cycles[(int)(divisor - Start)]; - } + internal CycleBlock(DivisorCycleCache owner, int index, ulong start, in ulong[] cycles) + { + _owner = owner; + Index = index; + Start = start; + End = start + (ulong)cycles.Length - 1UL; + Cycles = cycles; + } + + internal int Index { get; } + + internal ulong Start { get; } + + internal ulong End { get; } + + internal ulong[] Cycles { get; } + + internal int ReferenceCount => Volatile.Read(ref _referenceCount); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong GetCycle(ulong divisor) => Cycles[(int)(divisor - Start)]; + + internal void Retain() + { + Interlocked.Increment(ref _referenceCount); + } + + public void Dispose() + { + int remaining = Interlocked.Decrement(ref _referenceCount); + if (remaining < 0) + { + throw new InvalidOperationException("Cycle block released more times than it was acquired."); + } + + if (remaining == 0) + { + _owner.OnCycleBlockReleased(this); + } + } + } private const int CycleGenerationBatchSize = 262_144; private const byte ByteZero = 0; @@ -63,24 +100,37 @@ public sealed class DivisorCycleCache private void Initialize(ulong[] snapshot) { + CycleBlock? activeToDispose = null; + CycleBlock? prefetchedToDispose = null; + lock (_sync) { + if (_activeBlock is not null && _activeBlock.Index != 0) + { + activeToDispose = _activeBlock; + } + + if (_prefetchedBlock is not null && _prefetchedBlock.Index != 0) + { + prefetchedToDispose = _prefetchedBlock; + } + _pending.Clear(); - _baseBlock = new CycleBlock(index: 0, start: 0UL, snapshot); + _baseBlock = new CycleBlock(this, index: 0, start: 0UL, snapshot); _activeBlock = null; _prefetchedBlock = null; StartPrefetchLocked(1); } + + activeToDispose?.Dispose(); + prefetchedToDispose?.Dispose(); } public CycleBlock Acquire(ulong divisor) { - if (divisor <= _baseBlock.End) - { - return _baseBlock; - } - - return AcquireDynamicBlock(divisor); + CycleBlock block = divisor <= _baseBlock.End ? _baseBlock : AcquireDynamicBlock(divisor); + block.Retain(); + return block; } private CycleBlock AcquireDynamicBlock(ulong divisor) @@ -149,21 +199,42 @@ public sealed class DivisorCycleCache return; } + CycleBlock? previousActive = _activeBlock; _activeBlock = prefetched; _prefetchedBlock = null; + prefetched.Retain(); + prefetched.Dispose(); + + if (previousActive is not null && previousActive.Index != 0 && previousActive.Index != prefetched.Index) + { + previousActive.Dispose(); + } + StartPrefetchLocked(prefetched.Index + 1); } private void SetActiveBlockLocked(CycleBlock block) { + CycleBlock? previousActive = _activeBlock; _activeBlock = block; + block.Retain(); CycleBlock? prefetched = _prefetchedBlock; if (prefetched is not null && prefetched.Index != block.Index + 1) { + if (prefetched.Index != 0) + { + prefetched.Dispose(); + } + _prefetchedBlock = null; } + if (previousActive is not null && previousActive.Index != 0 && previousActive.Index != block.Index) + { + previousActive.Dispose(); + } + StartPrefetchLocked(block.Index + 1); } @@ -216,12 +287,59 @@ public sealed class DivisorCycleCache CycleBlock? active = _activeBlock; if (active is not null && blockIndex != active.Index + 1) { + MontgomeryDivisorDataCache.ReleaseBlock(block); return; } - if (_prefetchedBlock is null || _prefetchedBlock.Index < blockIndex) + CycleBlock? previousPrefetched = _prefetchedBlock; + if (ReferenceEquals(previousPrefetched, block)) { + return; + } + + if (_prefetchedBlock is null || _prefetchedBlock.Index <= blockIndex) + { + if (previousPrefetched is not null && previousPrefetched.Index != 0) + { + previousPrefetched.Dispose(); + } + _prefetchedBlock = block; + block.Retain(); + } + else + { + MontgomeryDivisorDataCache.ReleaseBlock(block); + } + } + + private void OnCycleBlockReleased(CycleBlock block) + { + if (block.Index == 0) + { + return; + } + + for (int attempt = 0; attempt < 8; attempt++) + { + if (block.ReferenceCount != 0) + { + return; + } + + Thread.SpinWait(1 << attempt); + } + + bool stillTracked; + + lock (_sync) + { + stillTracked = ReferenceEquals(_activeBlock, block) || ReferenceEquals(_prefetchedBlock, block); + } + + if (!stillTracked && block.ReferenceCount == 0) + { + MontgomeryDivisorDataCache.ReleaseBlock(block); } } @@ -245,7 +363,7 @@ public sealed class DivisorCycleCache ComputeCyclesCpu(start, cycles); } - return new CycleBlock(blockIndex, start, cycles); + return new CycleBlock(this, blockIndex, start, cycles); } private int GetBlockIndex(ulong divisor) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 534dddf..bd0c879 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -2,36 +2,33 @@ using System.Buffers; using System.Collections.Concurrent; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; +using System.Numerics; using ILGPU; using ILGPU.Runtime; +using MontgomeryDivisorData = PerfectNumbers.Core.MontgomeryDivisorData; +using MontgomeryDivisorDataCache = PerfectNumbers.Core.MontgomeryDivisorDataCache; +using CycleRemainderStepper = PerfectNumbers.Core.CycleRemainderStepper; namespace PerfectNumbers.Core.Gpu; public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDivisorByDivisorTester { - private int _gpuBatchSize = GpuConstants.ScanBatchSize; - private readonly object _sync = new(); - private ulong _divisorLimit; - private bool _isConfigured; - private bool _useDivisorCycles; + private int _gpuBatchSize = GpuConstants.ScanBatchSize; + private readonly object _sync = new(); + private ulong _divisorLimit; + private bool _isConfigured; + private bool _useDivisorCycles; - public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo) - { - public readonly ulong Modulus = modulus; - public readonly ulong NPrime = nPrime; - public readonly ulong MontgomeryOne = montgomeryOne; - public readonly ulong MontgomeryTwo = montgomeryTwo; - } - - private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); + private readonly ConcurrentDictionary, ArrayView, ArrayView>> _kernelCache = new(); private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _sessionPool = new(); - private Action, ArrayView> GetKernel(Accelerator accelerator) => - _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernel)); - private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputePrimeExponentKernel)); + private Action, ArrayView, ArrayView> GetKernel(Accelerator accelerator) => + _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(CheckKernel)); + + private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => + _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputePrimeExponentKernel)); public int GpuBatchSize { @@ -96,21 +93,23 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv try { - composite = CheckDivisors( - prime, - allowedMax, - useCycles, - accelerator, - kernel, - resources.DivisorsBuffer, - resources.HitsBuffer, - resources.Divisors, - resources.Hits, - resources.DivisorData, - out lastProcessed, - out coveredRange, - out processedCount - ); + composite = CheckDivisors( + prime, + allowedMax, + useCycles, + accelerator, + kernel, + resources.DivisorsBuffer, + resources.ExponentBuffer, + resources.HitsBuffer, + resources.Divisors, + resources.Exponents, + resources.Hits, + resources.DivisorData, + out lastProcessed, + out coveredRange, + out processedCount + ); } finally { @@ -128,17 +127,19 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return true; } - private static bool CheckDivisors( - ulong prime, - ulong allowedMax, - bool useCycles, - Accelerator accelerator, - Action, ArrayView> kernel, - MemoryBuffer1D divisorsBuffer, - MemoryBuffer1D hitsBuffer, - in ulong[] divisors, - byte[] hits, - MontgomeryDivisorData[] divisorData, + private static bool CheckDivisors( + ulong prime, + ulong allowedMax, + bool useCycles, + Accelerator accelerator, + Action, ArrayView, ArrayView> kernel, + MemoryBuffer1D divisorsBuffer, + MemoryBuffer1D exponentBuffer, + MemoryBuffer1D hitsBuffer, + in ulong[] divisors, + ulong[] exponents, + byte[] hits, + MontgomeryDivisorData[] divisorData, out ulong lastProcessed, out bool coveredRange, out ulong processedCount) @@ -154,20 +155,22 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv int batchCapacity = (int)divisorsBuffer.Length; - bool composite = false; - bool processedAll = false; - ulong currentDivisor = 3UL, cycle = 0UL, next, nextDivisor; - int batchSize, i; - bool reachedEndInBatch; + bool composite = false; + bool processedAll = false; + ulong currentDivisor = 3UL, next, nextDivisor; + int batchSize, i; + bool reachedEndInBatch; - Span divisorDataSpan; - Span divisorSpan; - Span hitsSpan; - ArrayView1D divisorView; - ArrayView1D hitsView; - DivisorCycleCache divisorCyclesCache = DivisorCycleCache.Shared; - DivisorCycleCache.CycleBlock? cycleLease = useCycles ? divisorCyclesCache.Acquire(currentDivisor) : null; - ulong cycleLeaseEnd = cycleLease?.End ?? 0UL; + Span divisorDataSpan; + Span divisorSpan; + Span exponentSpan; + Span hitsSpan; + ArrayView1D divisorView; + ArrayView1D exponentView; + ArrayView1D hitsView; + DivisorCycleCache divisorCyclesCache = DivisorCycleCache.Shared; + DivisorCycleCache.CycleBlock? cycleLease = useCycles ? divisorCyclesCache.Acquire(currentDivisor) : null; + ulong cycleLeaseEnd = cycleLease?.End ?? 0UL; while (currentDivisor <= allowedMax) { @@ -179,26 +182,39 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv nextDivisor = currentDivisor + 2UL; processedCount++; - bool includeDivisor = true; - if (useCycles) - { - if (currentDivisor > cycleLeaseEnd) - { - cycleLease = divisorCyclesCache.Acquire(currentDivisor); - cycleLeaseEnd = cycleLease.End; - } + bool includeDivisor = true; + ulong exponent = prime; - cycle = prime % cycleLease!.GetCycle(currentDivisor); - if (cycle != 0UL) - { - includeDivisor = false; - } - } + if (useCycles) + { + if (currentDivisor > cycleLeaseEnd) + { + cycleLease?.Dispose(); + cycleLease = divisorCyclesCache.Acquire(currentDivisor); + cycleLeaseEnd = cycleLease.End; + } - if (includeDivisor) - { - divisors[batchSize++] = currentDivisor; - } + ulong cycleLength = cycleLease!.GetCycle(currentDivisor); + if (cycleLength > 0UL) + { + ulong remainder = prime % cycleLength; + if (remainder != 0UL) + { + includeDivisor = false; + } + else + { + exponent = remainder; + } + } + } + + if (includeDivisor) + { + divisors[batchSize] = currentDivisor; + exponents[batchSize] = exponent; + batchSize++; + } lastProcessed = currentDivisor; next = nextDivisor; @@ -219,32 +235,34 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv currentDivisor = next; } - if (batchSize == 0) - { - if (reachedEndInBatch) - { - processedAll = true; - break; - } + if (batchSize == 0) + { + if (reachedEndInBatch) + { + processedAll = true; + break; + } - continue; - } + continue; + } - divisorDataSpan = divisorData.AsSpan(0, batchSize); - divisorView = divisorsBuffer.View.SubView(0, batchSize); - hitsView = hitsBuffer.View.SubView(0, batchSize); + divisorDataSpan = divisorData.AsSpan(0, batchSize); + divisorView = divisorsBuffer.View.SubView(0, batchSize); + exponentView = exponentBuffer.View.SubView(0, batchSize); + hitsView = hitsBuffer.View.SubView(0, batchSize); - divisorSpan = divisors.AsSpan(0, batchSize); - hitsSpan = hits.AsSpan(0, batchSize); + divisorSpan = divisors.AsSpan(0, batchSize); + exponentSpan = exponents.AsSpan(0, batchSize); + hitsSpan = hits.AsSpan(0, batchSize); - for (i = 0; i < batchSize; i++) - { - divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]); - } + for (i = 0; i < batchSize; i++) + { + divisorDataSpan[i] = MontgomeryDivisorDataCache.Get(divisorSpan[i]); + } - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); - byte useCycleChecks = useCycles ? (byte)1 : (byte)0; - kernel(batchSize, prime, useCycleChecks, cycle, divisorView, hitsView); + divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); + exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSpan), batchSize); + kernel(batchSize, divisorView, exponentView, hitsView); accelerator.Synchronize(); hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); @@ -274,9 +292,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } - coveredRange = composite || processedAll || currentDivisor > allowedMax; + cycleLease?.Dispose(); + coveredRange = composite || processedAll || currentDivisor > allowedMax; - return composite; + return composite; } private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) @@ -309,23 +328,19 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); } - private static void CheckKernel(Index1D index, ulong prime, byte useCycleChecks, ulong cycle, 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, ArrayView divisors, ArrayView exponents, ArrayView hits) + { + MontgomeryDivisorData divisor = divisors[index]; + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + hits[index] = 0; + return; + } - if (useCycleChecks != 0) - { - prime = cycle; - } - - hits[index] = prime.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; - } + ulong exponent = exponents[index]; + hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; + } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession { @@ -364,7 +379,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv int gpuBatchSize = owner._gpuBatchSize, primesLength = primes.Length; - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); ulong modulus = divisorData.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) { @@ -376,52 +391,60 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayView1D exponentView, exponentsView = _exponentsBuffer.View; ArrayView1D resultView, resultsView = _resultsBuffer.View; - Span exponentSlice, hostSpan = new(_hostBuffer); - Span positionSpan = new(_positionBuffer); - Span hitsSlice; - ReadOnlySpan primesSlice; + Span exponentSlice, hostSpan = new(_hostBuffer); + Span positionSpan = new(_positionBuffer); + Span hitsSlice; + ReadOnlySpan primesSlice; - int batchSize, computeCount, i, offset = 0; - bool useDivisorCycles = owner._useDivisorCycles; - ulong residue; - while (offset < primesLength) - { - batchSize = Math.Min(gpuBatchSize, primesLength - offset); - primesSlice = primes.Slice(offset, batchSize); - hitsSlice = hits.Slice(offset, batchSize); + int batchSize, computeCount, i, offset = 0; + bool useDivisorCycles = owner._useDivisorCycles && divisorCycle != 0UL; + CycleRemainderStepper cycleStepper = default; + if (useDivisorCycles) + { + cycleStepper = new CycleRemainderStepper(divisorCycle); + } + while (offset < primesLength) + { + batchSize = Math.Min(gpuBatchSize, primesLength - offset); + primesSlice = primes.Slice(offset, batchSize); + hitsSlice = hits.Slice(offset, batchSize); - computeCount = 0; + computeCount = 0; - if (useDivisorCycles) - { - for (i = 0; i < batchSize; i++) - { - residue = primesSlice[i] % divisorCycle; - if (residue != 0UL) - { - // We don't need to clear hitsSlice[i] because the parent function calls .Clear() on the buffer - continue; - } + if (useDivisorCycles) + { + Span remainderSpan = hostSpan.Slice(0, batchSize); - positionSpan[computeCount] = i; - hostSpan[computeCount] = residue; - computeCount++; - } - } - else - { - for (i = 0; i < batchSize; i++) - { - positionSpan[computeCount] = i; - hostSpan[computeCount] = primesSlice[i]; + for (i = 0; i < batchSize; i++) + { + ulong remainder = cycleStepper.ComputeNext(primesSlice[i]); + remainderSpan[i] = remainder; + hitsSlice[i] = 0; + + if (remainder != 0UL) + { + continue; + } + + positionSpan[computeCount] = i; + hostSpan[computeCount] = remainder; + computeCount++; + } + } + else + { + for (i = 0; i < batchSize; i++) + { + positionSpan[computeCount] = i; + hostSpan[computeCount] = primesSlice[i]; computeCount++; } } - if (computeCount > 0) - { - exponentSlice = hostSpan[..computeCount]; - exponentView = exponentsView.SubView(0, computeCount); + if (computeCount > 0) + { + exponentSlice = hostSpan[..computeCount]; + exponentView = exponentsView.SubView(0, computeCount); exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); resultView = resultsView.SubView(0, computeCount); @@ -429,14 +452,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv accelerator.Synchronize(); resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - for (i = 0; i < computeCount; i++) - { - hitsSlice[positionSpan[i]] = exponentSlice[i] == 1UL ? (byte)1 : (byte)0; - } - } + for (i = 0; i < computeCount; i++) + { + hitsSlice[positionSpan[i]] = exponentSlice[i] == 1UL ? (byte)1 : (byte)0; + } + } - offset += batchSize; - } + offset += batchSize; + } } public void Dispose() @@ -463,35 +486,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv results[index] = exponents[index].Pow2MontgomeryMod(divisor); } - private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) - { - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return new(modulus, 0UL, 0UL, 0UL); - } - - return new( - modulus, - ComputeMontgomeryNPrime(modulus), - ComputeMontgomeryResidue(1UL, modulus), - ComputeMontgomeryResidue(2UL, modulus) - ); - } - - 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); - } - public ulong DivisorLimit { get @@ -523,37 +517,45 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private sealed class BatchResources : IDisposable { - internal BatchResources(Accelerator accelerator, int capacity) - { - DivisorsBuffer = accelerator.Allocate1D(capacity); - HitsBuffer = accelerator.Allocate1D(capacity); - Divisors = ArrayPool.Shared.Rent(capacity); - Hits = ArrayPool.Shared.Rent(capacity); - DivisorData = ArrayPool.Shared.Rent(capacity); - Capacity = capacity; - } + internal BatchResources(Accelerator accelerator, int capacity) + { + DivisorsBuffer = accelerator.Allocate1D(capacity); + ExponentBuffer = accelerator.Allocate1D(capacity); + HitsBuffer = accelerator.Allocate1D(capacity); + Divisors = ArrayPool.Shared.Rent(capacity); + Exponents = ArrayPool.Shared.Rent(capacity); + Hits = ArrayPool.Shared.Rent(capacity); + DivisorData = ArrayPool.Shared.Rent(capacity); + Capacity = capacity; + } - internal readonly MemoryBuffer1D DivisorsBuffer; + internal readonly MemoryBuffer1D DivisorsBuffer; - internal readonly MemoryBuffer1D HitsBuffer; + internal readonly MemoryBuffer1D ExponentBuffer; - internal readonly ulong[] Divisors; + internal readonly MemoryBuffer1D HitsBuffer; - internal readonly byte[] Hits; + internal readonly ulong[] Divisors; - internal readonly MontgomeryDivisorData[] DivisorData; + internal readonly ulong[] Exponents; - internal readonly int Capacity; + internal readonly byte[] Hits; - public void Dispose() - { - DivisorsBuffer.Dispose(); - HitsBuffer.Dispose(); - ArrayPool.Shared.Return(Divisors, clearArray: false); - ArrayPool.Shared.Return(Hits, clearArray: false); - ArrayPool.Shared.Return(DivisorData, clearArray: false); - } - } + internal readonly MontgomeryDivisorData[] DivisorData; + + internal readonly int Capacity; + + public void Dispose() + { + DivisorsBuffer.Dispose(); + ExponentBuffer.Dispose(); + HitsBuffer.Dispose(); + ArrayPool.Shared.Return(Divisors, clearArray: false); + ArrayPool.Shared.Return(Exponents, clearArray: false); + ArrayPool.Shared.Return(Hits, clearArray: false); + ArrayPool.Shared.Return(DivisorData, clearArray: false); + } + } private sealed class AcceleratorReferenceComparer : IEqualityComparer { diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs index 289c442..6f3c2f0 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorGpuTester.cs @@ -13,25 +13,25 @@ public sealed class MersenneNumberDivisorGpuTester private Action> GetKernel(Accelerator accelerator) => _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel>(Kernel)); - public static void BuildDivisorCandidates() - { + public static void BuildDivisorCandidates() + { ulong[] snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - (ulong divisor, ulong cycle)[] list = new (ulong divisor, ulong cycle)[snapshot.Length / 2]; + (ulong divisor, uint cycle)[] list = new (ulong divisor, uint cycle)[snapshot.Length / 2]; ulong cycle; - int count = 0, i, snapshotLength = snapshot.Length; - for (i = 3; i < snapshotLength; i += 2) - { - cycle = snapshot[i]; - if (cycle == 0U) - { - continue; - } + int count = 0, i, snapshotLength = snapshot.Length; + for (i = 3; i < snapshotLength; i += 2) + { + cycle = snapshot[i]; + if (cycle == 0U) + { + continue; + } - list[count++] = ((ulong)i, cycle); - } + list[count++] = ((ulong)i, (uint)cycle); + } _divisorCandidates = count == 0 ? [] : list[..count]; - } + } public bool IsDivisible(ulong exponent, UInt128 divisor) { @@ -52,7 +52,7 @@ public sealed class MersenneNumberDivisorGpuTester return divisible; } - private static (ulong divisor, ulong cycle)[]? _divisorCandidates = Array.Empty<(ulong divisor, ulong cycle)>(); + private static (ulong divisor, uint cycle)[]? _divisorCandidates = Array.Empty<(ulong divisor, uint cycle)>(); public bool IsPrime(ulong p, UInt128 d, ulong divisorCyclesSearchLimit, out bool divisorsExhausted) { @@ -68,17 +68,17 @@ public sealed class MersenneNumberDivisorGpuTester return true; } - if (_divisorCandidates is { Length: > 0 } candidates) - { - int len = candidates.Length; - ulong cycle; - for (int k = 0; k < len; k++) - { - (ulong dSmall, cycle) = candidates[k]; - if (p % cycle != 0UL) - { - continue; - } + if (_divisorCandidates is { Length: > 0 } candidates) + { + int len = candidates.Length; + uint cycle; + for (int k = 0; k < len; k++) + { + (ulong dSmall, cycle) = candidates[k]; + if (p % cycle != 0UL) + { + continue; + } if (IsDivisible(p, dSmall)) { diff --git a/PerfectNumbers.Core/MontgomeryDivisorData.cs b/PerfectNumbers.Core/MontgomeryDivisorData.cs new file mode 100644 index 0000000..708cebd --- /dev/null +++ b/PerfectNumbers.Core/MontgomeryDivisorData.cs @@ -0,0 +1,148 @@ +using System.Collections.Concurrent; +using System.Numerics; +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +public readonly struct MontgomeryDivisorData(ulong modulus, ulong nPrime, ulong montgomeryOne, ulong montgomeryTwo) +{ + public ulong Modulus { get; } = modulus; + + public ulong NPrime { get; } = nPrime; + + public ulong MontgomeryOne { get; } = montgomeryOne; + + public ulong MontgomeryTwo { get; } = montgomeryTwo; +} + +internal static class MontgomeryDivisorDataCache +{ + private const ulong BaseBlockEnd = PerfectNumberConstants.MaxQForDivisorCycles; + + private static readonly ConcurrentDictionary Cache = new(); + private static readonly ConcurrentDictionary> BlockMembership = new(); + + public static MontgomeryDivisorData Get(ulong modulus) + { + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return new MontgomeryDivisorData(modulus, 0UL, 0UL, 0UL); + } + + int blockIndex = ComputeBlockIndex(modulus); + return GetManaged(modulus, blockIndex); + } + + internal static void ReleaseBlock(DivisorCycleCache.CycleBlock block) + { + if (block.Index <= 0) + { + return; + } + + if (BlockMembership.TryRemove(block.Index, out ConcurrentDictionary? members)) + { + foreach (ulong modulus in members.Keys) + { + Cache.TryRemove(modulus, out _); + } + } + } + + private static MontgomeryDivisorData GetManaged(ulong modulus, int blockIndex) + { + while (true) + { + if (Cache.TryGetValue(modulus, out CacheEntry entry)) + { + if (entry.BlockIndex != blockIndex) + { + if (blockIndex > 0 && entry.BlockIndex <= 0) + { + CacheEntry updated = new(entry.Data, blockIndex); + if (Cache.TryUpdate(modulus, updated, entry)) + { + RegisterBlockMembership(modulus, blockIndex); + return entry.Data; + } + + continue; + } + + if (blockIndex > 0) + { + RegisterBlockMembership(modulus, blockIndex); + } + + return entry.Data; + } + return entry.Data; + } + + MontgomeryDivisorData created = Create(modulus); + if (Cache.TryAdd(modulus, new CacheEntry(created, blockIndex))) + { + RegisterBlockMembership(modulus, blockIndex); + return created; + } + } + } + + private static void RegisterBlockMembership(ulong modulus, int blockIndex) + { + if (blockIndex <= 0) + { + return; + } + + ConcurrentDictionary set = BlockMembership.GetOrAdd(blockIndex, static _ => new ConcurrentDictionary()); + set.TryAdd(modulus, 0); + } + + private static int ComputeBlockIndex(ulong modulus) + { + if (modulus <= BaseBlockEnd) + { + return 0; + } + + return (int)((modulus - (BaseBlockEnd + 1UL)) / (ulong)PerfectNumberConstants.MaxQForDivisorCycles) + 1; + } + + private static MontgomeryDivisorData Create(ulong modulus) + { + 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 readonly struct CacheEntry + { + internal CacheEntry(MontgomeryDivisorData data, int blockIndex) + { + Data = data; + BlockIndex = blockIndex; + } + + internal MontgomeryDivisorData Data { get; } + + internal int BlockIndex { get; } + } +} diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index b02b507..839cbe3 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -214,7 +214,7 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryMod(this ulong exponent, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor) + public static ulong Pow2MontgomeryMod(this ulong exponent, in MontgomeryDivisorData divisor) { ulong modulus = divisor.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) @@ -247,7 +247,7 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MersenneNumberDivisorByDivisorGpuTester.MontgomeryDivisorData divisor) + public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) { if (cycleLength == 0UL) { @@ -255,12 +255,18 @@ public static class ULongExtensions } ulong rotationCount = exponent % cycleLength; - if (rotationCount == 0UL) + return rotationCount.Pow2MontgomeryModFromCycleRemainder(divisor); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModFromCycleRemainder(this ulong reducedExponent, in MontgomeryDivisorData divisor) + { + if (reducedExponent == 0UL) { return 1UL; } - return rotationCount.Pow2MontgomeryMod(divisor); + return reducedExponent.Pow2MontgomeryMod(divisor); } [MethodImpl(MethodImplOptions.AggressiveInlining)]