From df69dcbadb37686289af4d2a6d3d8698d7061020 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 24 Sep 2025 15:44:21 +0200 Subject: [PATCH] Improve by-divisor residue reuse and fix Montgomery mul --- EvenPerfectBitScanner/Program.cs | 11 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 259 ++++++++++++------ PerfectNumbers.Core/ULongExtensions.cs | 17 +- 3 files changed, 190 insertions(+), 97 deletions(-) diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index e08bb18..43b0526 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -813,7 +813,16 @@ internal static class Program return; } - states.Sort((left, right) => left.AllowedMax.CompareTo(right.AllowedMax)); + states.Sort(static (left, right) => + { + int compare = left.AllowedMax.CompareTo(right.AllowedMax); + if (compare != 0) + { + return compare; + } + + return left.Prime.CompareTo(right.Prime); + }); int stateCount = states.Count; ulong[] primeValues = new ulong[stateCount]; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 44c9ebc..5fab5ba 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; @@ -26,14 +27,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester } private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); - private readonly ConcurrentDictionary, ArrayView>> _kernelByPrimeCache = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelDeltaCache = 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> GetKernelByPrime(Accelerator accelerator) => - _kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckKernelByPrime)); + private Action, ArrayView> GetKernelByPrimeDelta(Accelerator accelerator) => + _kernelDeltaCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(ComputePrimeDeltaKernel)); public int GpuBatchSize { @@ -115,11 +116,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester gpuLease.Dispose(); } - if (processedCount > 0UL) - { - lock (_sync) - { - UpdateStatusUnsafe(lastProcessed, processedCount); + if (processedCount > 0UL) + { + lock (_sync) + { + UpdateStatusUnsafe(lastProcessed, processedCount); } } @@ -230,14 +231,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester accelerator.Synchronize(); hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize); - for (i = 0; i < batchSize; i++) - { - if (hitsSpan[i] != 0) - { - composite = true; - lastProcessed = divisorSpan[i]; - break; - } + for (i = 0; i < batchSize; i++) + { + if (hitsSpan[i] != 0) + { + composite = true; + lastProcessed = divisorSpan[i]; + break; + } } if (!composite) @@ -334,112 +335,198 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester hits[index] = prime.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; } - 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; + public sealed class DivisorScanSession : IDisposable + { + 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; - 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); + internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) + { + _owner = owner; + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = owner.GetKernelByPrimeDelta(_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); } - internal void Reset() - { - _disposed = false; - } + internal void Reset() + { + _disposed = false; + } - public void CheckDivisor(ulong divisor, ReadOnlySpan primes, Span hits) - { - if (_disposed) - { - throw new ObjectDisposedException(nameof(DivisorScanSession)); - } + 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); + + _capacity = requiredCapacity; + _exponentsBuffer = _accelerator.Allocate1D(_capacity); + _resultsBuffer = _accelerator.Allocate1D(_capacity); + _hostBuffer = ArrayPool.Shared.Rent(_capacity); + _positionBuffer = ArrayPool.Shared.Rent(_capacity); + } + + public void CheckDivisor(ulong divisor, ReadOnlySpan primes, Span hits) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(DivisorScanSession)); + } + + int primesLength = primes.Length; + if (primesLength == 0) + { + return; + } + + MersenneNumberDivisorByDivisorGpuTester owner = _owner; + int gpuBatchSize = owner._gpuBatchSize; + EnsureCapacity(gpuBatchSize); - int offset = 0; - MersenneNumberDivisorByDivisorGpuTester owner = _owner; - int gpuBatchSize = owner._gpuBatchSize; - bool useCycles = owner._useDivisorCycles; - Accelerator accelerator = _accelerator; - MemoryBuffer1D hitsBuffer = _hitsBuffer; - MemoryBuffer1D primesBuffer = _primesBuffer; - ArrayView1D primesBufferView = primesBuffer.View; - ArrayView1D hitsBufferView = hitsBuffer.View; - Action, ArrayView> kernel = _kernel; - ArrayView1D primesView; - ArrayView1D hitsView; - byte useCyclesFlag = useCycles ? (byte)1 : (byte)0; MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + ulong modulus = divisorData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return; + } + + bool useCycles = owner._useDivisorCycles; + ulong cycle = 0UL; + bool cycleEnabled = false; + if (useCycles) + { + cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); + cycleEnabled = cycle != 0UL; + } + + Accelerator accelerator = _accelerator; + Action, ArrayView> kernel = _kernel; + ArrayView1D exponentsView = _exponentsBuffer.View; + ArrayView1D resultsView = _resultsBuffer.View; + + int offset = 0; + int batchSize; ReadOnlySpan primesSlice; Span hitsSlice; - int batchSize, primesLength = primes.Length; + Span hostSpan = _hostBuffer.AsSpan(); + Span positionSpan = _positionBuffer.AsSpan(); + bool hasResidue = false; + ulong lastPrimeWithResidue = 0UL; + ulong lastResidue = 0UL; while (offset < primesLength) { batchSize = Math.Min(gpuBatchSize, primesLength - offset); primesSlice = primes.Slice(offset, batchSize); - hitsSlice = hits.Slice(offset, batchSize); + hitsSlice = hits.Slice(offset, batchSize); - primesView = primesBufferView.SubView(0, batchSize); - hitsView = hitsBufferView.SubView(0, batchSize); + int computeCount = 0; + ulong deltaBasePrime = lastPrimeWithResidue; + bool deltaHasResidue = hasResidue; - primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize); - kernel(batchSize, divisorData, useCyclesFlag, primesView, hitsView); - accelerator.Synchronize(); + for (int i = 0; i < batchSize; i++) + { + ulong prime = primesSlice[i]; + if (cycleEnabled && prime % cycle != 0UL) + { + hitsSlice[i] = 0; + continue; + } - hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSlice), batchSize); + ulong delta; + if (deltaHasResidue) + { + delta = prime >= deltaBasePrime ? prime - deltaBasePrime : prime; + } + else + { + delta = prime; + } + hostSpan[computeCount] = delta; + positionSpan[computeCount] = i; + computeCount++; + deltaBasePrime = prime; + deltaHasResidue = true; + } + + if (computeCount > 0) + { + Span exponentSlice = hostSpan.Slice(0, computeCount); + ArrayView1D exponentView = exponentsView.SubView(0, computeCount); + exponentView.CopyFromCPU(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 (int i = 0; i < computeCount; i++) + { + int position = positionSpan[i]; + ulong prime = primesSlice[position]; + ulong deltaResult = exponentSlice[i]; + ulong residue = hasResidue + ? lastResidue.MulMod64(deltaResult, modulus) + : deltaResult % modulus; + + hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0; + hasResidue = true; + lastResidue = residue; + lastPrimeWithResidue = prime; + } + } offset += batchSize; } } - public void Dispose() - { - if (_disposed) - { - return; - } + public void Dispose() + { + if (_disposed) + { + return; + } - _disposed = true; - _owner.ReturnSession(this); - } - } + _disposed = true; + _owner.ReturnSession(this); + } + + } internal void ReturnSession(DivisorScanSession session) { _sessionPool.Add(session); } - private static void CheckKernelByPrime(Index1D index, MontgomeryDivisorData divisor, byte useCycleChecks, ArrayView primes, ArrayView hits) + private static void ComputePrimeDeltaKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView deltas, ArrayView results) { ulong modulus = divisor.Modulus; if (modulus <= 1UL || (modulus & 1UL) == 0UL) { - hits[index] = 0; + results[index] = 0UL; return; } - if (useCycleChecks != 0) - { - ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); - if (cycle == 0UL || primes[index] % cycle != 0UL) - { - hits[index] = 0; - return; - } - } - - hits[index] = primes[index].Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; + ulong exponent = deltas[index]; + results[index] = exponent.Pow2MontgomeryMod(divisor); } private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 726dc35..3f2579c 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -163,18 +163,15 @@ public static class ULongExtensions { ulong xLow = (uint)x; ulong xHigh = x >> 32; - ulong yLow = (uint)y; - ulong yHigh = y >> 32; + ulong yLow = (uint)y; + ulong yHigh = y >> 32; - ulong w1 = xLow * yHigh; - ulong w2 = xHigh * yLow; + ulong w1 = xLow * yHigh; + ulong w2 = xHigh * yLow; - return (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32) + - ( - ((xLow * yLow) >> 32) + - (uint)w1 + - (uint)w2 - ) >> 32; + ulong result = (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32); + result += (((xLow * yLow) >> 32) + (uint)w1 + (uint)w2) >> 32; + return result; } [MethodImpl(MethodImplOptions.AggressiveInlining)]