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

Merge pull request #75 from mlemanczyk/codex/implement-todos-in-mersennenumberdivisorbydivisorgputester

Browse Source 
Reuse divisor residues across prime candidates
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-09-24 16:54:52 +02:00
2 parents a39f7dda59 436d3856d6
commit 283913375b
1 file changed
+136 -109
Download Patch File Download Diff File Expand all files Collapse all files
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+136 -109
View File
@@ -27,14 +27,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
} }
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, byte, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>> _kernelCache = new(); 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>>> _kernelDeltaCache = 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 ConcurrentBag<DivisorScanSession> _sessionPool = new(); private readonly ConcurrentBag<DivisorScanSession> _sessionPool = new();
private Action<Index1D, ulong, byte, ArrayView<MontgomeryDivisorData>, ArrayView<byte>> GetKernel(Accelerator accelerator) => 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)); _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, byte, ArrayView<MontgomeryDivisorData>, ArrayView<byte>>(CheckKernel));
private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeDelta(Accelerator accelerator) => private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeExponent(Accelerator accelerator) =>
_kernelDeltaCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputePrimeDeltaKernel)); _kernelExponentCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputePrimeExponentKernel));
public int GpuBatchSize public int GpuBatchSize
{ {
@@ -341,10 +341,11 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
private readonly GpuContextPool.GpuContextLease _lease; private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator; private readonly Accelerator _accelerator;
private readonly Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel; private readonly Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel;
private MemoryBuffer1D<ulong, Stride1D.Dense> _exponentsBuffer; private MemoryBuffer1D<ulong, Stride1D.Dense> _exponentsBuffer;
private MemoryBuffer1D<ulong, Stride1D.Dense> _resultsBuffer; private MemoryBuffer1D<ulong, Stride1D.Dense> _resultsBuffer;
private ulong[] _hostBuffer; private ulong[] _hostBuffer;
private int[] _positionBuffer; private ulong[] _primeBuffer;
private int[] _positionBuffer;
private int _capacity; private int _capacity;
private bool _disposed; private bool _disposed;
@@ -353,13 +354,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
_owner = owner; _owner = owner;
_lease = GpuContextPool.RentPreferred(preferCpu: false); _lease = GpuContextPool.RentPreferred(preferCpu: false);
_accelerator = _lease.Accelerator; _accelerator = _lease.Accelerator;
_kernel = owner.GetKernelByPrimeDelta(_accelerator); _kernel = owner.GetKernelByPrimeExponent(_accelerator);
_capacity = Math.Max(1, owner._gpuBatchSize); _capacity = Math.Max(1, owner._gpuBatchSize);
_exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity); _exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity); _resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity); _hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_positionBuffer = ArrayPool<int>.Shared.Rent(_capacity); _primeBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
} _positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
}
internal void Reset() internal void Reset()
{ {
@@ -374,16 +376,18 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
} }
_exponentsBuffer.Dispose(); _exponentsBuffer.Dispose();
_resultsBuffer.Dispose(); _resultsBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false); ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
ArrayPool<int>.Shared.Return(_positionBuffer, clearArray: false); ArrayPool<ulong>.Shared.Return(_primeBuffer, clearArray: false);
ArrayPool<int>.Shared.Return(_positionBuffer, clearArray: false);
_capacity = requiredCapacity; _capacity = requiredCapacity;
_exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity); _exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity); _resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity); _hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_positionBuffer = ArrayPool<int>.Shared.Rent(_capacity); _primeBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
} _positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
}
public void CheckDivisor(ulong divisor, ReadOnlySpan<ulong> primes, Span<byte> hits) public void CheckDivisor(ulong divisor, ReadOnlySpan<ulong> primes, Span<byte> hits)
{ {
@@ -409,14 +413,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
return; return;
} }
bool useCycles = owner._useDivisorCycles; bool useCycles = owner._useDivisorCycles;
ulong cycle = 0UL; ulong cycle = 0UL;
bool cycleEnabled = false; bool cycleEnabled = false;
if (useCycles) if (useCycles)
{ {
cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus); cycle = MersenneDivisorCycles.CalculateCycleLengthGpu(modulus);
cycleEnabled = cycle != 0UL; cycleEnabled = cycle != 0UL;
} }
Accelerator accelerator = _accelerator; Accelerator accelerator = _accelerator;
Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> kernel = _kernel; Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> kernel = _kernel;
@@ -427,108 +431,131 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
int batchSize; int batchSize;
ReadOnlySpan<ulong> primesSlice; ReadOnlySpan<ulong> primesSlice;
Span<byte> hitsSlice; Span<byte> hitsSlice;
Span<ulong> hostSpan = _hostBuffer.AsSpan(); Span<ulong> hostSpan = _hostBuffer.AsSpan();
Span<int> positionSpan = _positionBuffer.AsSpan(); Span<ulong> primeSpan = _primeBuffer.AsSpan();
bool hasResidue = false; Span<int> positionSpan = _positionBuffer.AsSpan();
ulong lastPrimeWithResidue = 0UL; bool hasPreviousResidue = false;
ulong lastResidue = 0UL; 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);
while (offset < primesLength) int computeCount = 0;
{
batchSize = Math.Min(gpuBatchSize, primesLength - offset);
primesSlice = primes.Slice(offset, batchSize);
hitsSlice = hits.Slice(offset, batchSize);
int computeCount = 0; for (int i = 0; i < batchSize; i++)
ulong deltaBasePrime = lastPrimeWithResidue; {
bool deltaHasResidue = hasResidue; ulong prime = primesSlice[i];
if (cycleEnabled && prime % cycle != 0UL)
for (int i = 0; i < batchSize; i++)
{
ulong prime = primesSlice[i];
if (cycleEnabled && prime % cycle != 0UL)
{ {
hitsSlice[i] = 0; hitsSlice[i] = 0;
continue; continue;
} }
ulong delta; primeSpan[computeCount] = prime;
if (deltaHasResidue) positionSpan[computeCount] = i;
{ computeCount++;
delta = prime >= deltaBasePrime ? prime - deltaBasePrime : prime; }
}
else
{
delta = prime;
}
hostSpan[computeCount] = delta;
positionSpan[computeCount] = i;
computeCount++;
deltaBasePrime = prime;
deltaHasResidue = true;
}
if (computeCount > 0) if (computeCount > 0)
{ {
Span<ulong> exponentSlice = hostSpan.Slice(0, computeCount); Span<ulong> exponentSlice = hostSpan.Slice(0, computeCount);
ArrayView1D<ulong, Stride1D.Dense> exponentView = exponentsView.SubView(0, computeCount); ulong currentPrime = hasPreviousResidue ? previousPrime : 0UL;
exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); bool hasDeltaSource = hasPreviousResidue;
ArrayView1D<ulong, Stride1D.Dense> resultView = resultsView.SubView(0, computeCount); for (int i = 0; i < computeCount; i++)
kernel(computeCount, divisorData, exponentView, resultView); {
accelerator.Synchronize(); ulong prime = primeSpan[i];
resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount); if (hasDeltaSource)
{
exponentSlice[i] = prime - currentPrime;
}
else
{
exponentSlice[i] = prime;
hasDeltaSource = true;
}
// TODO: Can we move this somehow to GPU? MulMod64 takes half of the time in CheckDivisor currentPrime = prime;
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; ArrayView1D<ulong, Stride1D.Dense> exponentView = exponentsView.SubView(0, computeCount);
hasResidue = true; exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount);
lastResidue = residue; ArrayView1D<ulong, Stride1D.Dense> resultView = resultsView.SubView(0, computeCount);
lastPrimeWithResidue = prime; kernel(computeCount, divisorData, exponentView, resultView);
} accelerator.Synchronize();
} resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount);
offset += batchSize; ulong lastPrime = previousPrime;
} ulong residueAccumulator = previousResidue;
} bool hasAccumulator = hasPreviousResidue;
for (int i = 0; i < computeCount; i++)
{
int position = positionSpan[i];
ulong stepResidue = exponentSlice[i];
ulong prime = primeSpan[i];
ulong residue = hasAccumulator ? MultiplyMod(residueAccumulator, stepResidue, modulus) : stepResidue;
hitsSlice[position] = residue == 1UL ? (byte)1 : (byte)0;
exponentSlice[i] = residue;
residueAccumulator = residue;
lastPrime = prime;
hasAccumulator = true;
}
public void Dispose() if (hasAccumulator)
{ {
if (_disposed) previousResidue = residueAccumulator;
{ previousPrime = lastPrime;
return; hasPreviousResidue = true;
} }
}
_disposed = true; offset += batchSize;
_owner.ReturnSession(this); }
} }
} public void Dispose()
{
if (_disposed)
{
return;
}
_disposed = true;
_owner.ReturnSession(this);
}
private static ulong MultiplyMod(ulong left, ulong right, ulong modulus)
{
if (modulus == 0UL)
{
return 0UL;
}
UInt128 product = (UInt128)left * right;
return (ulong)(product % modulus);
}
}
internal void ReturnSession(DivisorScanSession session) internal void ReturnSession(DivisorScanSession session)
{ {
_sessionPool.Add(session); _sessionPool.Add(session);
} }
private static void ComputePrimeDeltaKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> deltas, ArrayView<ulong> results) private static void ComputePrimeExponentKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> exponents, ArrayView<ulong> results)
{ {
ulong modulus = divisor.Modulus; ulong modulus = divisor.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL) if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{ {
results[index] = 0UL; results[index] = 0UL;
return; return;
} }
ulong exponent = deltas[index]; ulong exponent = exponents[index];
results[index] = exponent.Pow2MontgomeryMod(divisor); results[index] = exponent.Pow2MontgomeryMod(divisor);
} }
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
{ {