mlemanczyk/even-perfect-numbers-scanner
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Implement GPU-safe fast division helper

Browse Source
This commit is contained in:
Marek Lemańczyk committed 2025-09-24 21:07:43 +02:00
1 parent cc7cb12f86
commit e316c0055f
2 files changed
+108 -54
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/Gpu/GpuKernelPool.cs
+19 -3
View File
@@ -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)
{
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+89 -51
View File
@@ -28,13 +28,16 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, byte, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>> _kernelCache = new();
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>> _kernelExponentCache = new();
private readonly ConcurrentDictionary<Accelerator, ConcurrentBag<BatchResources>> _resourcePools = new(AcceleratorReferenceComparer.Instance);
private readonly ConcurrentDictionary<Accelerator, ConcurrentBag<BatchResources>> _resourcePools = new(AcceleratorReferenceComparer.Instance);
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _cycleKernelCache = new();
private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
private Action<Index1D, ulong, byte, ArrayView<MontgomeryDivisorData>, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, byte, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>(CheckKernel));
private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeExponent(Accelerator accelerator) =>
_kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputePrimeExponentKernel));
private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> GetCycleKernel(Accelerator accelerator) =>
_cycleKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(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<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel;
private readonly MersenneNumberDivisorByDivisorGpuTester _owner;
private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator;
private readonly Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel;
private readonly Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> _cycleKernel;
private MemoryBuffer1D<ulong, Stride1D.Dense> _exponentsBuffer;
private MemoryBuffer1D<ulong, Stride1D.Dense> _resultsBuffer;
private MemoryBuffer1D<byte, Stride1D.Dense> _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<ulong>(_capacity);
_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_cycleKernel = owner.GetCycleKernel(_accelerator);
_capacity = Math.Max(1, owner._gpuBatchSize);
_exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_cycleMaskBuffer = _accelerator.Allocate1D<byte>(_capacity);
_hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_primeBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
_cycleMaskHost = ArrayPool<byte>.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<ulong>.Shared.Return(_hostBuffer, clearArray: false);
ArrayPool<ulong>.Shared.Return(_primeBuffer, clearArray: false);
ArrayPool<int>.Shared.Return(_positionBuffer, clearArray: false);
ArrayPool<byte>.Shared.Return(_cycleMaskHost, clearArray: false);
_capacity = requiredCapacity;
_exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_cycleMaskBuffer = _accelerator.Allocate1D<byte>(_capacity);
_hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_primeBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
_cycleMaskHost = ArrayPool<byte>.Shared.Rent(_capacity);
}
public void CheckDivisor(ulong divisor, ulong divisorCycle, ReadOnlySpan<ulong> primes, Span<byte> 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<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> kernel = _kernel;
ArrayView1D<ulong, Stride1D.Dense> exponentsView = _exponentsBuffer.View;
ArrayView1D<ulong, Stride1D.Dense> resultsView = _resultsBuffer.View;
Accelerator accelerator = _accelerator;
Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> kernel = _kernel;
Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> cycleKernel = _cycleKernel;
ArrayView1D<ulong, Stride1D.Dense> exponentsView = _exponentsBuffer.View;
ArrayView1D<ulong, Stride1D.Dense> resultsView = _resultsBuffer.View;
ArrayView1D<byte, Stride1D.Dense> cycleMaskViewFull = _cycleMaskBuffer.View;
int offset = 0;
int batchSize;
ReadOnlySpan<ulong> primesSlice;
Span<byte> hitsSlice;
int offset = 0;
Span<ulong> hostSpan = _hostBuffer.AsSpan();
Span<ulong> primeSpan = _primeBuffer.AsSpan();
Span<int> 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<ulong> primesSlice = primes.Slice(offset, batchSize);
Span<byte> hitsSlice = hits.Slice(offset, batchSize);
int computeCount = 0;
Span<byte> cycleMaskSpan = default;
if (cycleEnabled)
{
cycleMaskSpan = _cycleMaskHost.AsSpan(0, batchSize);
ArrayView1D<ulong, Stride1D.Dense> primeInputView = _resultsBuffer.View.SubView(0, batchSize);
primeInputView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize);
ArrayView1D<byte, Stride1D.Dense> 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<ulong> 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<ulong> primes, ArrayView<byte> 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<ulong> exponents, ArrayView<ulong> results)
{
ulong modulus = divisor.Modulus;