diff --git a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs index 139a513..357f250 100644 --- a/PerfectNumbers.Core/Gpu/GpuKernelPool.cs +++ b/PerfectNumbers.Core/Gpu/GpuKernelPool.cs @@ -275,7 +275,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } - ulong div = phi64.FastDiv64(exponent, divMul); + ulong div = FastDiv64Gpu(phi64, exponent, divMul); GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, q) - GpuUInt128.One; if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) { @@ -311,7 +311,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } - ulong div = phi64.FastDiv64(exponent, divMul); + ulong div = FastDiv64Gpu(phi64, exponent, divMul); GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, q) - GpuUInt128.One; if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) { @@ -432,6 +432,22 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne orders[index] = exponent; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong FastDiv64Gpu(ulong value, ulong divisor, ulong mul) + { + ulong quotient = GpuUInt128.MulHigh(value, mul); + GpuUInt128 remainder = new(0UL, value); + GpuUInt128 product = new(GpuUInt128.MulHigh(quotient, divisor), quotient * divisor); + remainder.Sub(product); + + if (remainder.High != 0UL || remainder.Low >= divisor) + { + quotient++; + } + + return quotient; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong Mod128By64(GpuUInt128 value, ulong modulus) { @@ -552,7 +568,7 @@ public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, Kerne return; } - ulong div = phi64.FastDiv64(exponent, divMul); + ulong div = FastDiv64Gpu(phi64, exponent, divMul); GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, q) - GpuUInt128.One; if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One) { diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 9762ff9..a820e13 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -28,13 +28,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); - private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); + private readonly ConcurrentDictionary, ArrayView>> _cycleKernelCache = new(); 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> GetCycleKernel(Accelerator accelerator) => + _cycleKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel, ArrayView>(CheckPrimeCycleKernel)); public int GpuBatchSize { @@ -337,56 +340,66 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester 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 MersenneNumberDivisorByDivisorGpuTester _owner; + private readonly GpuContextPool.GpuContextLease _lease; + private readonly Accelerator _accelerator; + private readonly Action, ArrayView> _kernel; + private readonly Action, ArrayView> _cycleKernel; private MemoryBuffer1D _exponentsBuffer; private MemoryBuffer1D _resultsBuffer; + private MemoryBuffer1D _cycleMaskBuffer; private ulong[] _hostBuffer; private ulong[] _primeBuffer; private int[] _positionBuffer; - private int _capacity; - private bool _disposed; + private byte[] _cycleMaskHost; + private int _capacity; + private bool _disposed; - internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) - { - _owner = owner; - _lease = GpuContextPool.RentPreferred(preferCpu: false); - _accelerator = _lease.Accelerator; + internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) + { + _owner = owner; + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; _kernel = owner.GetKernelByPrimeExponent(_accelerator); - _capacity = Math.Max(1, owner._gpuBatchSize); - _exponentsBuffer = _accelerator.Allocate1D(_capacity); - _resultsBuffer = _accelerator.Allocate1D(_capacity); + _cycleKernel = owner.GetCycleKernel(_accelerator); + _capacity = Math.Max(1, owner._gpuBatchSize); + _exponentsBuffer = _accelerator.Allocate1D(_capacity); + _resultsBuffer = _accelerator.Allocate1D(_capacity); + _cycleMaskBuffer = _accelerator.Allocate1D(_capacity); _hostBuffer = ArrayPool.Shared.Rent(_capacity); _primeBuffer = ArrayPool.Shared.Rent(_capacity); _positionBuffer = ArrayPool.Shared.Rent(_capacity); + _cycleMaskHost = ArrayPool.Shared.Rent(_capacity); } - internal void Reset() - { - _disposed = false; + internal void Reset() + { + _disposed = false; } private void EnsureCapacity(int requiredCapacity) { - if (requiredCapacity <= _capacity) - { - return; - } + if (requiredCapacity <= _capacity) + { + return; + } - _exponentsBuffer.Dispose(); + _exponentsBuffer.Dispose(); _resultsBuffer.Dispose(); + _cycleMaskBuffer.Dispose(); ArrayPool.Shared.Return(_hostBuffer, clearArray: false); ArrayPool.Shared.Return(_primeBuffer, clearArray: false); ArrayPool.Shared.Return(_positionBuffer, clearArray: false); + ArrayPool.Shared.Return(_cycleMaskHost, clearArray: false); _capacity = requiredCapacity; _exponentsBuffer = _accelerator.Allocate1D(_capacity); _resultsBuffer = _accelerator.Allocate1D(_capacity); + _cycleMaskBuffer = _accelerator.Allocate1D(_capacity); _hostBuffer = ArrayPool.Shared.Rent(_capacity); _primeBuffer = ArrayPool.Shared.Rent(_capacity); _positionBuffer = ArrayPool.Shared.Rent(_capacity); + _cycleMaskHost = ArrayPool.Shared.Rent(_capacity); } public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan primes, Span hits) @@ -396,57 +409,69 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester throw new ObjectDisposedException(nameof(DivisorScanSession)); } - int primesLength = primes.Length; - if (primesLength == 0) - { - return; - } + int primesLength = primes.Length; + if (primesLength == 0) + { + return; + } - MersenneNumberDivisorByDivisorGpuTester owner = _owner; - int gpuBatchSize = owner._gpuBatchSize; - EnsureCapacity(gpuBatchSize); + MersenneNumberDivisorByDivisorGpuTester owner = _owner; + int gpuBatchSize = owner._gpuBatchSize; + EnsureCapacity(gpuBatchSize); - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return; - } + MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + ulong modulus = divisorData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return; + } bool cycleEnabled = owner._useDivisorCycles && divisorCycle != 0UL; - Accelerator accelerator = _accelerator; - Action, ArrayView> kernel = _kernel; - ArrayView1D exponentsView = _exponentsBuffer.View; - ArrayView1D resultsView = _resultsBuffer.View; + Accelerator accelerator = _accelerator; + Action, ArrayView> kernel = _kernel; + Action, ArrayView> cycleKernel = _cycleKernel; + ArrayView1D exponentsView = _exponentsBuffer.View; + ArrayView1D resultsView = _resultsBuffer.View; + ArrayView1D cycleMaskViewFull = _cycleMaskBuffer.View; - int offset = 0; - int batchSize; - ReadOnlySpan primesSlice; - Span hitsSlice; + int offset = 0; Span hostSpan = _hostBuffer.AsSpan(); Span primeSpan = _primeBuffer.AsSpan(); Span positionSpan = _positionBuffer.AsSpan(); bool hasPreviousResidue = false; ulong previousResidue = 0UL; ulong previousPrime = 0UL; + 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); int computeCount = 0; + Span cycleMaskSpan = default; + + if (cycleEnabled) + { + cycleMaskSpan = _cycleMaskHost.AsSpan(0, batchSize); + ArrayView1D primeInputView = _resultsBuffer.View.SubView(0, batchSize); + primeInputView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize); + ArrayView1D cycleMaskView = cycleMaskViewFull.SubView(0, batchSize); + cycleKernel(batchSize, divisorCycle, primeInputView, cycleMaskView); + accelerator.Synchronize(); + cycleMaskView.CopyToCPU(ref MemoryMarshal.GetReference(cycleMaskSpan), batchSize); + } for (int i = 0; i < batchSize; i++) { - ulong prime = primesSlice[i]; - if (cycleEnabled && prime % divisorCycle != 0UL) + if (cycleEnabled && cycleMaskSpan[i] == 0) { hitsSlice[i] = 0; continue; } + ulong prime = primesSlice[i]; primeSpan[computeCount] = prime; positionSpan[computeCount] = i; computeCount++; @@ -457,6 +482,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester Span exponentSlice = hostSpan.Slice(0, computeCount); ulong currentPrime = hasPreviousResidue ? previousPrime : 0UL; bool hasDeltaSource = hasPreviousResidue; + for (int i = 0; i < computeCount; i++) { ulong prime = primeSpan[i]; @@ -480,9 +506,10 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester accelerator.Synchronize(); resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); - ulong lastPrime = previousPrime; ulong residueAccumulator = previousResidue; + ulong lastPrime = previousPrime; bool hasAccumulator = hasPreviousResidue; + for (int i = 0; i < computeCount; i++) { int position = positionSpan[i]; @@ -537,6 +564,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester _sessionPool.Add(session); } + private static void CheckPrimeCycleKernel(Index1D index, ulong cycle, ArrayView primes, ArrayView mask) + { + if (cycle == 0UL) + { + mask[index] = 1; + return; + } + + mask[index] = primes[index] % cycle == 0UL ? (byte)1 : (byte)0; + } + private static void ComputePrimeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView exponents, ArrayView results) { ulong modulus = divisor.Modulus;