From 49fcd937479c5f6945fab4b68eb84e9ee8ccb5db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Thu, 25 Sep 2025 05:28:23 +0200 Subject: [PATCH] GPU logic updated to use prime % divisor result precomputed on CPU instead of re-calculating it on GPU side. Cycle lease removed. --- PerfectNumbers.Core/DivisorCycleCache.cs | 70 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 1362 ++++++++--------- 2 files changed, 684 insertions(+), 748 deletions(-) diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 460db80..dd14fd5 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -12,61 +12,16 @@ namespace PerfectNumbers.Core; public sealed class DivisorCycleCache { - internal sealed class CycleBlock - { - internal CycleBlock(int index, ulong start, ulong[] cycles) - { - Index = index; - Start = start; - Cycles = cycles; - End = start + (ulong)cycles.Length - 1UL; - } + 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; - internal int Index { get; } - - internal ulong Start { get; } - - internal ulong End { get; } - - internal ulong[] Cycles { get; } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal ulong GetCycle(ulong divisor) - { - if (divisor < Start || divisor > End) - { - return 0UL; - } - - return Cycles[(int)(divisor - Start)]; - } - } - - public struct Lease : IDisposable - { - private CycleBlock? _block; - - internal Lease(CycleBlock block) - { - _block = block; - } - - public readonly bool IsValid => _block is not null; - - public readonly ulong Start => _block?.Start ?? 0UL; - - public readonly ulong End => _block?.End ?? 0UL; - - public readonly ulong[]? Values => _block?.Cycles; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public readonly ulong GetCycle(ulong divisor) => _block?.GetCycle(divisor) ?? 0UL; - - public void Dispose() - { - _block = null; - } - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ulong GetCycle(ulong divisor) => Cycles[(int)(divisor - Start)]; + } private const int CycleGenerationBatchSize = 262_144; private const int GpuCycleStepsPerInvocation = 524_288 / 16; @@ -115,15 +70,14 @@ public sealed class DivisorCycleCache } } - public Lease Acquire(ulong divisor) + public CycleBlock Acquire(ulong divisor) { if (divisor <= _baseBlock.End) { - return new Lease(_baseBlock); + return _baseBlock; } - CycleBlock block = AcquireDynamicBlock(divisor); - return new Lease(block); + return AcquireDynamicBlock(divisor); } private CycleBlock AcquireDynamicBlock(ulong divisor) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 6915f25..a5931ae 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -11,727 +11,709 @@ 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; - } + 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>> _kernelExponentCache = new(); - private readonly ConcurrentDictionary> _resourcePools = new(AcceleratorReferenceComparer.Instance); - private readonly ConcurrentDictionary, ArrayView>> _cycleKernelCache = new(); - private readonly ConcurrentBag _sessionPool = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelCache = new(); + private readonly ConcurrentDictionary, ArrayView>> _kernelExponentCache = new(); + 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)); + 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 - { - 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; - } + int IMersenneNumberDivisorByDivisorTester.BatchSize + { + get => GpuBatchSize; + set => GpuBatchSize = value; + } - public bool UseDivisorCycles - { - get => _useDivisorCycles; - set => _useDivisorCycles = value; - } + public bool UseDivisorCycles + { + get => _useDivisorCycles; + set => _useDivisorCycles = value; + } - public void ConfigureFromMaxPrime(ulong maxPrime) - { - lock (_sync) - { - _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); - _isConfigured = true; - } - } + public void ConfigureFromMaxPrime(ulong maxPrime) + { + lock (_sync) + { + _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); + _isConfigured = true; + } + } - public bool IsPrime(ulong prime, out bool divisorsExhausted) - { - ulong allowedMax; - bool useCycles; - int batchCapacity; + 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."); - } + lock (_sync) + { + if (!_isConfigured) + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } - allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); - useCycles = _useDivisorCycles; - batchCapacity = _gpuBatchSize; - } + allowedMax = ComputeAllowedMaxDivisor(prime, _divisorLimit); + useCycles = _useDivisorCycles; + batchCapacity = _gpuBatchSize; + } - if (allowedMax < 3UL) - { - divisorsExhausted = true; - return true; - } + if (allowedMax < 3UL) + { + divisorsExhausted = true; + return true; + } - bool composite; - bool coveredRange; - ulong processedCount; - ulong lastProcessed; + bool composite; + bool coveredRange; + ulong processedCount; + ulong lastProcessed; - var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); - var accelerator = gpuLease.Accelerator; - var kernel = GetKernel(accelerator); - BatchResources resources = RentBatchResources(accelerator, batchCapacity); + var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); + var accelerator = gpuLease.Accelerator; + var kernel = GetKernel(accelerator); + BatchResources resources = RentBatchResources(accelerator, batchCapacity); - try - { - composite = CheckDivisors( - prime, - allowedMax, - useCycles, - accelerator, - kernel, - resources.DivisorsBuffer, - resources.HitsBuffer, - resources.Divisors, - resources.Hits, - resources.DivisorData, - out lastProcessed, - out coveredRange, - out processedCount - ); - } - finally - { - ReturnBatchResources(accelerator, resources); - gpuLease.Dispose(); - } + try + { + composite = CheckDivisors( + prime, + allowedMax, + useCycles, + accelerator, + kernel, + resources.DivisorsBuffer, + resources.HitsBuffer, + resources.Divisors, + resources.Hits, + resources.DivisorData, + out lastProcessed, + out coveredRange, + out processedCount + ); + } + finally + { + ReturnBatchResources(accelerator, resources); + gpuLease.Dispose(); + } - if (composite) - { - divisorsExhausted = true; - return false; - } + if (composite) + { + divisorsExhausted = true; + return false; + } - divisorsExhausted = coveredRange; - return true; - } + divisorsExhausted = coveredRange; + return true; + } - private static bool CheckDivisors( - ulong prime, - ulong allowedMax, - bool useCycles, - Accelerator accelerator, - Action, ArrayView> kernel, - MemoryBuffer1D divisorsBuffer, - MemoryBuffer1D hitsBuffer, - ulong[] divisors, - byte[] hits, - MontgomeryDivisorData[] divisorData, - out ulong lastProcessed, - out bool coveredRange, - out ulong processedCount) - { - lastProcessed = 0UL; - processedCount = 0UL; + private static bool CheckDivisors( + ulong prime, + ulong allowedMax, + bool useCycles, + Accelerator accelerator, + Action, ArrayView> kernel, + MemoryBuffer1D divisorsBuffer, + MemoryBuffer1D hitsBuffer, + ulong[] divisors, + byte[] hits, + MontgomeryDivisorData[] divisorData, + out ulong lastProcessed, + out bool coveredRange, + out ulong processedCount) + { + lastProcessed = 0UL; + processedCount = 0UL; - if (allowedMax < 3UL) - { - coveredRange = true; - return false; - } + if (allowedMax < 3UL) + { + coveredRange = true; + return false; + } - int batchCapacity = (int)divisorsBuffer.Length; + int batchCapacity = (int)divisorsBuffer.Length; - bool composite = false; - bool processedAll = false; - ulong currentDivisor, divisor = 3UL, next, nextDivisor; - int batchSize, i; - bool reachedEndInBatch; - bool requiresCycleKernelCheck; + bool composite = false; + bool processedAll = false; + ulong currentDivisor = 3UL, cycle = 0UL, next, nextDivisor; + int batchSize, i; + bool reachedEndInBatch; - Span divisorDataSpan; - Span divisorSpan; - Span hitsSpan; - ArrayView1D divisorView; - ArrayView1D hitsView; + Span divisorDataSpan; + Span divisorSpan; + Span hitsSpan; + ArrayView1D divisorView; + ArrayView1D hitsView; DivisorCycleCache divisorCyclesCache = DivisorCycleCache.Shared; + DivisorCycleCache.CycleBlock? cycleLease = useCycles ? divisorCyclesCache.Acquire(currentDivisor) : null; - while (divisor <= allowedMax) - { - batchSize = 0; - reachedEndInBatch = false; - requiresCycleKernelCheck = false; + while (currentDivisor <= allowedMax) + { + batchSize = 0; + reachedEndInBatch = false; - while (batchSize < batchCapacity && divisor <= allowedMax) - { - currentDivisor = divisor; - nextDivisor = currentDivisor + 2UL; - processedCount++; + while (batchSize < batchCapacity && currentDivisor <= allowedMax) + { + nextDivisor = currentDivisor + 2UL; + processedCount++; - bool includeDivisor = true; + bool includeDivisor = true; if (useCycles) { - DivisorCycleCache.Lease cycleLease = divisorCyclesCache.Acquire(currentDivisor); - ulong cycle = cycleLease.GetCycle(currentDivisor); - if (cycle == 0UL) + if (currentDivisor > cycleLease!.End) { - requiresCycleKernelCheck = true; + cycleLease = divisorCyclesCache.Acquire(currentDivisor); } - else if (prime % cycle != 0UL) + + cycle = prime % cycleLease.GetCycle(currentDivisor); + if (cycle != 0UL) { includeDivisor = false; } - - cycleLease.Dispose(); - } - - if (includeDivisor) - { - divisors[batchSize++] = currentDivisor; - } - - lastProcessed = currentDivisor; - next = nextDivisor; - - if (next <= divisor) - { - reachedEndInBatch = true; - break; - } - - if (next > allowedMax) - { - divisor = next; - reachedEndInBatch = true; - break; - } - - divisor = next; - } - - if (batchSize == 0) - { - if (reachedEndInBatch) - { - processedAll = true; - break; - } - - continue; - } - - divisorDataSpan = divisorData.AsSpan(0, batchSize); - divisorView = divisorsBuffer.View.SubView(0, batchSize); - hitsView = hitsBuffer.View.SubView(0, batchSize); - - divisorSpan = divisors.AsSpan(0, batchSize); - hitsSpan = hits.AsSpan(0, batchSize); - - for (i = 0; i < batchSize; i++) - { - divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]); - } - - divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); - byte useCycleChecks = useCycles && requiresCycleKernelCheck ? (byte)1 : (byte)0; - kernel(batchSize, prime, useCycleChecks, divisorView, hitsView); - 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; - } - } - - if (!composite) - { - lastProcessed = divisorSpan[batchSize - 1]; - } - if (composite) - { - break; - } - - if (reachedEndInBatch) - { - processedAll = true; - break; - } - } - - coveredRange = composite || processedAll || divisor > allowedMax; - - return composite; - } - - 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 void CheckKernel(Index1D index, ulong prime, byte useCycleChecks, ArrayView divisors, ArrayView hits) - { - MontgomeryDivisorData divisor = divisors[index]; - ulong modulus = divisor.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - hits[index] = 0; - return; - } - - ulong exponent = prime; - if (useCycleChecks != 0) - { - ulong cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); - if (cycle == 0UL) - { - hits[index] = 0; - return; - } - ulong remainder = prime % cycle; - if (remainder != 0UL) - { - hits[index] = 0; - return; - } - - exponent = remainder; - } - - hits[index] = exponent.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; - } - - public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession - { - 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 byte[] _cycleMaskHost; - private int _capacity; - private bool _disposed; - - internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) - { - _owner = owner; - _lease = GpuContextPool.RentPreferred(preferCpu: false); - _accelerator = _lease.Accelerator; - _kernel = owner.GetKernelByPrimeExponent(_accelerator); - _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; - } - - private void EnsureCapacity(int requiredCapacity) - { - if (requiredCapacity <= _capacity) - { - return; - } - - _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) - { - if (_disposed) - { - throw new ObjectDisposedException(nameof(DivisorScanSession)); - } - - int primesLength = primes.Length; - if (primesLength == 0) - { - return; - } - - MersenneNumberDivisorByDivisorGpuTester owner = _owner; - int gpuBatchSize = owner._gpuBatchSize; - EnsureCapacity(gpuBatchSize); - - MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); - ulong modulus = divisorData.Modulus; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - return; - } - - ulong cycle = divisorCycle; - bool cycleEnabled = owner._useDivisorCycles && cycle != 0UL; - - 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; - Span hostSpan = _hostBuffer.AsSpan(); - Span primeSpan = _primeBuffer.AsSpan(); - Span positionSpan = _positionBuffer.AsSpan(); - - while (offset < primesLength) - { - int batchSize = Math.Min(gpuBatchSize, primesLength - offset); - ReadOnlySpan primesSlice = primes.Slice(offset, batchSize); - Span hitsSlice = hits.Slice(offset, batchSize); - - int computeCount = 0; - - if (cycleEnabled) - { - Span 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, cycle, primeInputView, cycleMaskView); - accelerator.Synchronize(); - cycleMaskView.CopyToCPU(ref MemoryMarshal.GetReference(cycleMaskSpan), batchSize); - - for (int i = 0; i < batchSize; i++) - { - if (cycleMaskSpan[i] == 0) - { - hitsSlice[i] = 0; - continue; - } - - primeSpan[computeCount] = primesSlice[i]; - positionSpan[computeCount] = i; - computeCount++; - } - } - else - { - for (int i = 0; i < batchSize; i++) - { - primeSpan[computeCount] = primesSlice[i]; - positionSpan[computeCount] = i; - computeCount++; - } - } - - if (computeCount > 0) - { - Span exponentSlice = hostSpan[..computeCount]; - for (int i = 0; i < computeCount; i++) - { - ulong prime = primeSpan[i]; - exponentSlice[i] = cycleEnabled ? prime % cycle : prime; - } - - 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 residue = exponentSlice[i]; - hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0; - } - } - - offset += batchSize; - } - } - - public void Dispose() - { - if (_disposed) - { - return; - } - - _disposed = true; - _owner.ReturnSession(this); - } - - - } - - internal void ReturnSession(DivisorScanSession session) - { - _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; - if (modulus <= 1UL || (modulus & 1UL) == 0UL) - { - results[index] = 0UL; - return; - } - - ulong exponent = exponents[index]; - results[index] = exponent.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 - { - lock (_sync) - { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - return _divisorLimit; - } - } - } - - private BatchResources RentBatchResources(Accelerator accelerator, int capacity) - { - var bag = _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()); - if (bag.TryTake(out BatchResources? resources)) - { - resources.EnsureCapacity(capacity); - return resources; - } - - return new BatchResources(accelerator, capacity); - } - - private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) - { - _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()).Add(resources); - } - - private sealed class BatchResources : IDisposable - { - private readonly Accelerator _accelerator; - - internal BatchResources(Accelerator accelerator, int capacity) - { - _accelerator = accelerator; - 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 MemoryBuffer1D DivisorsBuffer { get; private set; } - - internal MemoryBuffer1D HitsBuffer { get; private set; } - - internal ulong[] Divisors { get; private set; } - - internal byte[] Hits { get; private set; } - - internal MontgomeryDivisorData[] DivisorData { get; private set; } - - internal int Capacity { get; private set; } - - public void EnsureCapacity(int requiredCapacity) - { - if (requiredCapacity <= Capacity) - { - return; - } - - Resize(requiredCapacity); - } - - private void Resize(int newCapacity) - { - DivisorsBuffer.Dispose(); - HitsBuffer.Dispose(); - ArrayPool.Shared.Return(Divisors, clearArray: false); - ArrayPool.Shared.Return(Hits, clearArray: false); - ArrayPool.Shared.Return(DivisorData, clearArray: false); - - DivisorsBuffer = _accelerator.Allocate1D(newCapacity); - HitsBuffer = _accelerator.Allocate1D(newCapacity); - Divisors = ArrayPool.Shared.Rent(newCapacity); - Hits = ArrayPool.Shared.Rent(newCapacity); - DivisorData = ArrayPool.Shared.Rent(newCapacity); - Capacity = newCapacity; - } - - 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); - } - } - - private sealed class AcceleratorReferenceComparer : IEqualityComparer - { - public static AcceleratorReferenceComparer Instance { get; } = new(); - - public bool Equals(Accelerator? x, Accelerator? y) => ReferenceEquals(x, y); - - public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj); - } - - 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 DivisorScanSession 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); - } - } - - IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() - { - return CreateDivisorSession(); - } + } + + if (includeDivisor) + { + divisors[batchSize++] = currentDivisor; + } + + lastProcessed = currentDivisor; + next = nextDivisor; + + if (next <= currentDivisor) + { + reachedEndInBatch = true; + break; + } + + if (next > allowedMax) + { + currentDivisor = next; + reachedEndInBatch = true; + break; + } + + currentDivisor = next; + } + + if (batchSize == 0) + { + if (reachedEndInBatch) + { + processedAll = true; + break; + } + + continue; + } + + divisorDataSpan = divisorData.AsSpan(0, batchSize); + divisorView = divisorsBuffer.View.SubView(0, batchSize); + hitsView = hitsBuffer.View.SubView(0, batchSize); + + divisorSpan = divisors.AsSpan(0, batchSize); + hitsSpan = hits.AsSpan(0, batchSize); + + for (i = 0; i < batchSize; i++) + { + divisorDataSpan[i] = CreateMontgomeryDivisorData(divisorSpan[i]); + } + + divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorDataSpan), batchSize); + byte useCycleChecks = useCycles ? (byte)1 : (byte)0; + kernel(batchSize, prime, useCycleChecks, cycle, divisorView, hitsView); + 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; + } + } + + if (!composite) + { + lastProcessed = divisorSpan[batchSize - 1]; + } + if (composite) + { + break; + } + + if (reachedEndInBatch) + { + processedAll = true; + break; + } + } + + coveredRange = composite || processedAll || currentDivisor > allowedMax; + + return composite; + } + + 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 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; + } + + if (useCycleChecks != 0) + { + prime = cycle; + } + + hits[index] = prime.Pow2MontgomeryMod(divisor) == 1UL ? (byte)1 : (byte)0; + } + + public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession + { + 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 byte[] _cycleMaskHost; + private int _capacity; + private bool _disposed; + + internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) + { + _owner = owner; + _lease = GpuContextPool.RentPreferred(preferCpu: false); + _accelerator = _lease.Accelerator; + _kernel = owner.GetKernelByPrimeExponent(_accelerator); + _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; + } + + private void EnsureCapacity(int requiredCapacity) + { + if (requiredCapacity <= _capacity) + { + return; + } + + _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) + { + if (_disposed) + { + throw new ObjectDisposedException(nameof(DivisorScanSession)); + } + + int primesLength = primes.Length; + if (primesLength == 0) + { + return; + } + + MersenneNumberDivisorByDivisorGpuTester owner = _owner; + int gpuBatchSize = owner._gpuBatchSize; + EnsureCapacity(gpuBatchSize); + + MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); + ulong modulus = divisorData.Modulus; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + return; + } + + ulong cycle = divisorCycle; + bool cycleEnabled = owner._useDivisorCycles && cycle != 0UL; + + 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; + Span hostSpan = _hostBuffer.AsSpan(); + Span primeSpan = _primeBuffer.AsSpan(); + Span positionSpan = _positionBuffer.AsSpan(); + + while (offset < primesLength) + { + int batchSize = Math.Min(gpuBatchSize, primesLength - offset); + ReadOnlySpan primesSlice = primes.Slice(offset, batchSize); + Span hitsSlice = hits.Slice(offset, batchSize); + + int computeCount = 0; + + if (cycleEnabled) + { + Span 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, cycle, primeInputView, cycleMaskView); + accelerator.Synchronize(); + cycleMaskView.CopyToCPU(ref MemoryMarshal.GetReference(cycleMaskSpan), batchSize); + + for (int i = 0; i < batchSize; i++) + { + if (cycleMaskSpan[i] == 0) + { + hitsSlice[i] = 0; + continue; + } + + primeSpan[computeCount] = primesSlice[i]; + positionSpan[computeCount] = i; + computeCount++; + } + } + else + { + for (int i = 0; i < batchSize; i++) + { + primeSpan[computeCount] = primesSlice[i]; + positionSpan[computeCount] = i; + computeCount++; + } + } + + if (computeCount > 0) + { + Span exponentSlice = hostSpan[..computeCount]; + for (int i = 0; i < computeCount; i++) + { + ulong prime = primeSpan[i]; + exponentSlice[i] = cycleEnabled ? prime % cycle : prime; + } + + 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 residue = exponentSlice[i]; + hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0; + } + } + + offset += batchSize; + } + } + + public void Dispose() + { + if (_disposed) + { + return; + } + + _disposed = true; + _owner.ReturnSession(this); + } + + + } + + internal void ReturnSession(DivisorScanSession session) + { + _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; + if (modulus <= 1UL || (modulus & 1UL) == 0UL) + { + results[index] = 0UL; + return; + } + + ulong exponent = exponents[index]; + results[index] = exponent.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 + { + lock (_sync) + { + if (!_isConfigured) + { + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); + } + + return _divisorLimit; + } + } + } + + private BatchResources RentBatchResources(Accelerator accelerator, int capacity) + { + var bag = _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()); + if (bag.TryTake(out BatchResources? resources)) + { + resources.EnsureCapacity(capacity); + return resources; + } + + return new BatchResources(accelerator, capacity); + } + + private void ReturnBatchResources(Accelerator accelerator, BatchResources resources) + { + _resourcePools.GetOrAdd(accelerator, static _ => new ConcurrentBag()).Add(resources); + } + + private sealed class BatchResources : IDisposable + { + private readonly Accelerator _accelerator; + + internal BatchResources(Accelerator accelerator, int capacity) + { + _accelerator = accelerator; + 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 MemoryBuffer1D DivisorsBuffer { get; private set; } + + internal MemoryBuffer1D HitsBuffer { get; private set; } + + internal ulong[] Divisors { get; private set; } + + internal byte[] Hits { get; private set; } + + internal MontgomeryDivisorData[] DivisorData { get; private set; } + + internal int Capacity { get; private set; } + + public void EnsureCapacity(int requiredCapacity) + { + if (requiredCapacity <= Capacity) + { + return; + } + + Resize(requiredCapacity); + } + + private void Resize(int newCapacity) + { + DivisorsBuffer.Dispose(); + HitsBuffer.Dispose(); + ArrayPool.Shared.Return(Divisors, clearArray: false); + ArrayPool.Shared.Return(Hits, clearArray: false); + ArrayPool.Shared.Return(DivisorData, clearArray: false); + + DivisorsBuffer = _accelerator.Allocate1D(newCapacity); + HitsBuffer = _accelerator.Allocate1D(newCapacity); + Divisors = ArrayPool.Shared.Rent(newCapacity); + Hits = ArrayPool.Shared.Rent(newCapacity); + DivisorData = ArrayPool.Shared.Rent(newCapacity); + Capacity = newCapacity; + } + + 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); + } + } + + private sealed class AcceleratorReferenceComparer : IEqualityComparer + { + public static AcceleratorReferenceComparer Instance { get; } = new(); + + public bool Equals(Accelerator? x, Accelerator? y) => ReferenceEquals(x, y); + + public int GetHashCode(Accelerator obj) => RuntimeHelpers.GetHashCode(obj); + } + + 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 DivisorScanSession 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); + } + } + + IMersenneNumberDivisorByDivisorTester.IDivisorScanSession IMersenneNumberDivisorByDivisorTester.CreateDivisorSession() + { + return CreateDivisorSession(); + } }