From 9b7611fe62748a579cf013b2183ac87fde29a01d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 25 Sep 2025 20:23:25 +0200 Subject: [PATCH] Trim GPU exponent batches using cached divisor cycles --- EvenPerfectBitScanner/Program.cs | 44 +-- ...MersenneNumberDivisorByDivisorCpuTester.cs | 13 +- .../ExponentRemainderStepper.cs | 124 +++++++++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 251 ++++++++++++------ PerfectNumbers.Core/ULongExtensions.cs | 32 +++ 5 files changed, 353 insertions(+), 111 deletions(-) create mode 100644 PerfectNumbers.Core/ExponentRemainderStepper.cs diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 60b3a3d..07bfa60 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -811,30 +811,32 @@ internal static class Program _byDivisorPrecheckOnly = false; - if (states.Count == 0) - { - if (applyStartPrime) - { - Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); - } + if (states.Count == 0) + { + if (applyStartPrime) + { + Console.WriteLine($"No primes greater than or equal to {startPrime.ToString(CultureInfo.InvariantCulture)} were found for --mersenne=bydivisor."); + } - return; - } + return; + } - int stateCount = states.Count; - ulong[] primeValues = new ulong[stateCount]; - ulong[] allowedMaxValues = new ulong[stateCount]; - int[] stateFlags = new int[stateCount]; + states.Sort(static (left, right) => left.Prime.CompareTo(right.Prime)); - int stateIndex = 0; - ByDivisorPrimeState currentState; - for (; stateIndex < stateCount; stateIndex++) - { - currentState = states[stateIndex]; - primeValues[stateIndex] = currentState.Prime; - allowedMaxValues[stateIndex] = currentState.AllowedMax; - stateFlags[stateIndex] = ByDivisorStateActive; - } + int stateCount = states.Count; + ulong[] primeValues = new ulong[stateCount]; + ulong[] allowedMaxValues = new ulong[stateCount]; + int[] stateFlags = new int[stateCount]; + + int stateIndex = 0; + ByDivisorPrimeState currentState; + for (; stateIndex < stateCount; stateIndex++) + { + currentState = states[stateIndex]; + primeValues[stateIndex] = currentState.Prime; + allowedMaxValues[stateIndex] = currentState.AllowedMax; + stateFlags[stateIndex] = ByDivisorStateActive; + } states.Clear(); diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 53a6417..965161f 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -310,6 +310,13 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv bool cycleEnabled = _owner._useDivisorCycles && divisorCycle != 0UL; MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); + var exponentStepper = new ExponentRemainderStepper(divisorData); + if (!exponentStepper.IsValidModulus) + { + hits.Clear(); + return; + } + if (cycleEnabled) { var stepper = new CycleRemainderStepper(divisorCycle); @@ -323,8 +330,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv continue; } - ulong residue = remainder.Pow2MontgomeryModFromCycleRemainder(divisorData); - hits[i] = residue == 1UL ? (byte)1 : (byte)0; + hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; } return; @@ -332,8 +338,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv for (int i = 0; i < length; i++) { - ulong residue = primes[i].Pow2MontgomeryMod(divisorData); - hits[i] = residue == 1UL ? (byte)1 : (byte)0; + hits[i] = exponentStepper.ComputeNextIsUnity(primes[i]) ? (byte)1 : (byte)0; } } diff --git a/PerfectNumbers.Core/ExponentRemainderStepper.cs b/PerfectNumbers.Core/ExponentRemainderStepper.cs new file mode 100644 index 0000000..2a16d7e --- /dev/null +++ b/PerfectNumbers.Core/ExponentRemainderStepper.cs @@ -0,0 +1,124 @@ +namespace PerfectNumbers.Core; + +internal struct ExponentRemainderStepper +{ + private readonly MontgomeryDivisorData _divisor; + private readonly ulong _modulus; + private readonly ulong _nPrime; + private readonly ulong _montgomeryOne; + private ulong _previousExponent; + private ulong _currentMontgomery; + private bool _hasState; + + public ExponentRemainderStepper(in MontgomeryDivisorData divisor) + { + _divisor = divisor; + _modulus = divisor.Modulus; + _nPrime = divisor.NPrime; + _montgomeryOne = divisor.MontgomeryOne; + _previousExponent = 0UL; + _currentMontgomery = divisor.MontgomeryOne; + _hasState = false; + } + + public bool IsValidModulus => _modulus > 1UL && (_modulus & 1UL) != 0UL; + + public bool HasState => _hasState; + + public ulong PreviousExponent => _previousExponent; + + public void Reset() + { + _previousExponent = 0UL; + _currentMontgomery = _montgomeryOne; + _hasState = false; + } + + public ulong ComputeNext(ulong exponent) + { + if (!IsValidModulus) + { + _hasState = false; + return 0UL; + } + + if (!_hasState || exponent <= _previousExponent) + { + _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); + _previousExponent = exponent; + _hasState = true; + return ReduceCurrent(); + } + + ulong delta = exponent - _previousExponent; + ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor); + _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); + _previousExponent = exponent; + return ReduceCurrent(); + } + + public bool ComputeNextIsUnity(ulong exponent) + { + if (!IsValidModulus) + { + _hasState = false; + return false; + } + + if (!_hasState || exponent <= _previousExponent) + { + _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); + _previousExponent = exponent; + _hasState = true; + return _currentMontgomery == _montgomeryOne; + } + + ulong delta = exponent - _previousExponent; + ulong multiplier = delta.Pow2MontgomeryModMontgomery(_divisor); + _currentMontgomery = _currentMontgomery.MontgomeryMultiply(multiplier, _modulus, _nPrime); + _previousExponent = exponent; + return _currentMontgomery == _montgomeryOne; + } + + public bool TryInitializeFromMontgomeryResult(ulong exponent, ulong montgomeryResult, out bool isUnity) + { + if (!IsValidModulus) + { + _hasState = false; + isUnity = false; + return false; + } + + _currentMontgomery = montgomeryResult; + _previousExponent = exponent; + _hasState = true; + isUnity = _currentMontgomery == _montgomeryOne; + return true; + } + + public bool TryAdvanceWithMontgomeryDelta(ulong exponent, ulong montgomeryDelta, out bool isUnity) + { + if (!IsValidModulus) + { + _hasState = false; + isUnity = false; + return false; + } + + if (!_hasState || exponent <= _previousExponent) + { + _currentMontgomery = exponent.Pow2MontgomeryModMontgomery(_divisor); + _previousExponent = exponent; + _hasState = true; + isUnity = _currentMontgomery == _montgomeryOne; + return true; + } + + _currentMontgomery = _currentMontgomery.MontgomeryMultiply(montgomeryDelta, _modulus, _nPrime); + _previousExponent = exponent; + isUnity = _currentMontgomery == _montgomeryOne; + return true; + } + + private ulong ReduceCurrent() => _currentMontgomery.MontgomeryMultiply(1UL, _modulus, _nPrime); +} diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index bd0c879..eba49b2 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -1,3 +1,4 @@ +using System; using System.Buffers; using System.Collections.Concurrent; using System.Runtime.CompilerServices; @@ -20,7 +21,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private bool _useDivisorCycles; private readonly ConcurrentDictionary, ArrayView, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); private readonly ConcurrentBag _sessionPool = new(); @@ -28,7 +29,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView, ArrayView>(CheckKernel)); private Action, ArrayView> GetKernelByPrimeExponent(Accelerator accelerator) => - _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputePrimeExponentKernel)); + _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputeMontgomeryExponentKernel)); public int GpuBatchSize { @@ -347,25 +348,27 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly GpuContextPool.GpuContextLease _lease; private readonly Accelerator _accelerator; - private readonly Action, ArrayView> _kernel; - private MemoryBuffer1D _exponentsBuffer; - private MemoryBuffer1D _resultsBuffer; - private ulong[] _hostBuffer; - private int[] _positionBuffer; - private int _capacity; - private bool _disposed; + private readonly Action, ArrayView> _kernel; + private MemoryBuffer1D _exponentsBuffer; + private MemoryBuffer1D _resultsBuffer; + private ulong[] _hostBuffer; + private int[] _positionBuffer; + private byte[] _modeBuffer; + private int _capacity; + private bool _disposed; internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) { _owner = owner; _lease = GpuContextPool.RentPreferred(preferCpu: false); _accelerator = _lease.Accelerator; - _kernel = owner.GetKernelByPrimeExponent(_accelerator); + _kernel = owner.GetKernelByPrimeExponent(_accelerator); _capacity = Math.Max(1, owner._gpuBatchSize); _exponentsBuffer = _accelerator.Allocate1D(_capacity); _resultsBuffer = _accelerator.Allocate1D(_capacity); - _hostBuffer = ArrayPool.Shared.Rent(_capacity); - _positionBuffer = ArrayPool.Shared.Rent(_capacity); + _hostBuffer = ArrayPool.Shared.Rent(_capacity); + _positionBuffer = ArrayPool.Shared.Rent(_capacity); + _modeBuffer = ArrayPool.Shared.Rent(_capacity); } internal void Reset() @@ -373,94 +376,168 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv _disposed = false; } - public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) - { - MersenneNumberDivisorByDivisorGpuTester owner = _owner; - int gpuBatchSize = owner._gpuBatchSize, - primesLength = primes.Length; + public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(DivisorScanSession)); + } + + MersenneNumberDivisorByDivisorGpuTester owner = _owner; + int gpuBatchSize = owner._gpuBatchSize; + int primesLength = primes.Length; + if (primesLength == 0) + { + return; + } + + EnsureCapacity(gpuBatchSize); MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(divisor); - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return; - } - - Accelerator accelerator = _accelerator; - Action, ArrayView> kernel = _kernel; - ArrayView1D exponentView, exponentsView = _exponentsBuffer.View; - ArrayView1D resultView, resultsView = _resultsBuffer.View; - - Span exponentSlice, hostSpan = new(_hostBuffer); - Span positionSpan = new(_positionBuffer); - Span hitsSlice; - ReadOnlySpan primesSlice; - - int batchSize, computeCount, i, offset = 0; - bool useDivisorCycles = owner._useDivisorCycles && divisorCycle != 0UL; - CycleRemainderStepper cycleStepper = default; - if (useDivisorCycles) + var exponentStepper = new ExponentRemainderStepper(divisorData); + if (!exponentStepper.IsValidModulus) { - cycleStepper = new CycleRemainderStepper(divisorCycle); + hits.Clear(); + return; } + + bool useDivisorCycles = owner._useDivisorCycles; + if (useDivisorCycles && divisorCycle == 0UL) + { + throw new InvalidOperationException($"Missing divisor cycle for divisor {divisor}."); + } + + CycleRemainderStepper cycleStepper = useDivisorCycles ? new CycleRemainderStepper(divisorCycle) : default; + + Accelerator accelerator = _accelerator; + Action, ArrayView> kernel = _kernel; + ArrayView1D exponentsView = _exponentsBuffer.View; + ArrayView1D resultsView = _resultsBuffer.View; + + Span hostSpan = _hostBuffer.AsSpan(0, _capacity); + Span positionSpan = _positionBuffer.AsSpan(0, _capacity); + Span modeSpan = _modeBuffer.AsSpan(0, _capacity); + + int offset = 0; + while (offset < primesLength) { - batchSize = Math.Min(gpuBatchSize, primesLength - offset); - primesSlice = primes.Slice(offset, batchSize); - hitsSlice = hits.Slice(offset, batchSize); + int batchSize = Math.Min(gpuBatchSize, primesLength - offset); + ReadOnlySpan primesSlice = primes.Slice(offset, batchSize); + Span hitsSlice = hits.Slice(offset, batchSize); + hitsSlice.Clear(); - computeCount = 0; + int computeCount = 0; + bool localHasState = exponentStepper.HasState; + ulong localPreviousExponent = exponentStepper.PreviousExponent; - if (useDivisorCycles) + for (int i = 0; i < batchSize; i++) { - Span remainderSpan = hostSpan.Slice(0, batchSize); + ulong primeValue = primesSlice[i]; - for (i = 0; i < batchSize; i++) + if (useDivisorCycles) { - ulong remainder = cycleStepper.ComputeNext(primesSlice[i]); - remainderSpan[i] = remainder; - hitsSlice[i] = 0; - + ulong remainder = cycleStepper.ComputeNext(primeValue); if (remainder != 0UL) { continue; } - - positionSpan[computeCount] = i; - hostSpan[computeCount] = remainder; - computeCount++; } - } - else - { - for (i = 0; i < batchSize; i++) + + ulong workValue; + byte mode; + + if (!localHasState || primeValue <= localPreviousExponent) { - positionSpan[computeCount] = i; - hostSpan[computeCount] = primesSlice[i]; - computeCount++; - } - } + mode = 0; + workValue = primeValue; + localHasState = true; + } + else + { + mode = 1; + workValue = primeValue - localPreviousExponent; + } + + localPreviousExponent = primeValue; + + if (useDivisorCycles) + { + workValue %= divisorCycle; + } + + positionSpan[computeCount] = i; + hostSpan[computeCount] = workValue; + modeSpan[computeCount] = mode; + computeCount++; + } if (computeCount > 0) { - exponentSlice = hostSpan[..computeCount]; - exponentView = exponentsView.SubView(0, computeCount); + Span exponentSlice = hostSpan[..computeCount]; + ArrayView1D exponentView = exponentsView.SubView(0, computeCount); + exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - resultView = resultsView.SubView(0, computeCount); - kernel(computeCount, divisorData, exponentView, resultView); - accelerator.Synchronize(); - resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); + ArrayView1D resultView = resultsView.SubView(0, computeCount); + kernel(computeCount, divisorData, exponentView, resultView); + accelerator.Synchronize(); + resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - for (i = 0; i < computeCount; i++) + for (int i = 0; i < computeCount; i++) { - hitsSlice[positionSpan[i]] = exponentSlice[i] == 1UL ? (byte)1 : (byte)0; + int position = positionSpan[i]; + ulong primeValue = primesSlice[position]; + ulong montgomeryResult = exponentSlice[i]; + bool isUnity; + + if (modeSpan[i] == 0) + { + exponentStepper.TryInitializeFromMontgomeryResult(primeValue, montgomeryResult, out isUnity); + } + else + { + exponentStepper.TryAdvanceWithMontgomeryDelta(primeValue, montgomeryResult, out isUnity); + } + + hitsSlice[position] = isUnity ? (byte)1 : (byte)0; } } offset += batchSize; } - } + } + + private void EnsureCapacity(int requiredCapacity) + { + if (requiredCapacity <= _capacity) + { + return; + } + + _exponentsBuffer.Dispose(); + _resultsBuffer.Dispose(); + ArrayPool.Shared.Return(_hostBuffer, clearArray: false); + ArrayPool.Shared.Return(_positionBuffer, clearArray: false); + ArrayPool.Shared.Return(_modeBuffer, clearArray: false); + + int newCapacity = _capacity; + if (newCapacity < 1) + { + newCapacity = 1; + } + + while (newCapacity < requiredCapacity) + { + newCapacity <<= 1; + } + + _exponentsBuffer = _accelerator.Allocate1D(newCapacity); + _resultsBuffer = _accelerator.Allocate1D(newCapacity); + _hostBuffer = ArrayPool.Shared.Rent(newCapacity); + _positionBuffer = ArrayPool.Shared.Rent(newCapacity); + _modeBuffer = ArrayPool.Shared.Rent(newCapacity); + _capacity = newCapacity; + } public void Dispose() { @@ -469,22 +546,24 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return; } - _disposed = true; - _owner._sessionPool.Add(this); - } - } + _disposed = true; + _owner._sessionPool.Add(this); + } + } - private static void ComputePrimeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) - { - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[index] = 0UL; - return; - } + private static void ComputeMontgomeryExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) + { + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + results[index] = 0UL; + return; + } - results[index] = exponents[index].Pow2MontgomeryMod(divisor); - } + ulong exponent = exponents[index]; + + results[index] = exponent.Pow2MontgomeryModMontgomery(divisor); + } public ulong DivisorLimit { diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 839cbe3..2c31787 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -246,6 +246,38 @@ public static class ULongExtensions return result.MontgomeryMultiply(1UL, modulus, nPrime); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static ulong Pow2MontgomeryModMontgomery(this ulong exponent, in MontgomeryDivisorData divisor) + { + ulong modulus = divisor.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return 0UL; + } + + ulong result = divisor.MontgomeryOne; + ulong baseVal = divisor.MontgomeryTwo; + ulong remainingExponent = exponent; + + while (remainingExponent > 0UL) + { + if ((remainingExponent & 1UL) != 0UL) + { + result = result.MontgomeryMultiply(baseVal, modulus, divisor.NPrime); + } + + remainingExponent >>= 1; + if (remainingExponent == 0UL) + { + break; + } + + baseVal = baseVal.MontgomeryMultiply(baseVal, modulus, divisor.NPrime); + } + + return result; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Pow2MontgomeryModWithCycle(this ulong exponent, ulong cycleLength, in MontgomeryDivisorData divisor) {