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

Merge pull request #74 from mlemanczyk/codex/fix-evenperfectbitscanner-for-mersenne-mode

Browse Source 
Improve by-divisor residue reuse and fix Montgomery multiplication
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-09-24 15:44:49 +02:00
2 parents 82a0fb1dd7 df69dcbadb
commit 74b41aad85
3 files changed
+139 -46
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

EvenPerfectBitScanner/Program.cs
+10 -1
View File
@@ -813,7 +813,16 @@ internal static class Program
return; 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; int stateCount = states.Count;
ulong[] primeValues = new ulong[stateCount]; ulong[] primeValues = new ulong[stateCount];
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+126 -39
View File
@@ -1,3 +1,4 @@
using System;
using System.Buffers; using System.Buffers;
using System.Collections.Concurrent; using System.Collections.Concurrent;
using System.Runtime.CompilerServices; using System.Runtime.CompilerServices;
@@ -26,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, byte, ArrayView<ulong>, ArrayView<byte>>> _kernelByPrimeCache = new(); private readonly ConcurrentDictionary<Accelerator, Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>> _kernelDeltaCache = 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, byte, ArrayView<ulong>, ArrayView<byte>> GetKernelByPrime(Accelerator accelerator) => private Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> GetKernelByPrimeDelta(Accelerator accelerator) =>
_kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, byte, ArrayView<ulong>, ArrayView<byte>>(CheckKernelByPrime)); _kernelDeltaCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>>(ComputePrimeDeltaKernel));
public int GpuBatchSize public int GpuBatchSize
{ {
@@ -339,9 +340,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
private readonly MersenneNumberDivisorByDivisorGpuTester _owner; private readonly MersenneNumberDivisorByDivisorGpuTester _owner;
private readonly GpuContextPool.GpuContextLease _lease; private readonly GpuContextPool.GpuContextLease _lease;
private readonly Accelerator _accelerator; private readonly Accelerator _accelerator;
private readonly Action<Index1D, MontgomeryDivisorData, byte, ArrayView<ulong>, ArrayView<byte>> _kernel; private readonly Action<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> _kernel;
private readonly MemoryBuffer1D<ulong, Stride1D.Dense> _primesBuffer; private MemoryBuffer1D<ulong, Stride1D.Dense> _exponentsBuffer;
private readonly MemoryBuffer1D<byte, Stride1D.Dense> _hitsBuffer; private MemoryBuffer1D<ulong, Stride1D.Dense> _resultsBuffer;
private ulong[] _hostBuffer;
private int[] _positionBuffer;
private int _capacity;
private bool _disposed; private bool _disposed;
internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner)
@@ -349,9 +353,12 @@ 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.GetKernelByPrime(_accelerator); _kernel = owner.GetKernelByPrimeDelta(_accelerator);
_primesBuffer = _accelerator.Allocate1D<ulong>(owner._gpuBatchSize); _capacity = Math.Max(1, owner._gpuBatchSize);
_hitsBuffer = _accelerator.Allocate1D<byte>(owner._gpuBatchSize); _exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_hostBuffer = ArrayPool<ulong>.Shared.Rent(_capacity);
_positionBuffer = ArrayPool<int>.Shared.Rent(_capacity);
} }
internal void Reset() internal void Reset()
@@ -359,6 +366,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
_disposed = false; _disposed = false;
} }
private void EnsureCapacity(int requiredCapacity)
{
if (requiredCapacity <= _capacity)
{
return;
}
_exponentsBuffer.Dispose();
_resultsBuffer.Dispose();
ArrayPool<ulong>.Shared.Return(_hostBuffer, clearArray: false);
ArrayPool<int>.Shared.Return(_positionBuffer, clearArray: false);
_capacity = requiredCapacity;
_exponentsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_resultsBuffer = _accelerator.Allocate1D<ulong>(_capacity);
_hostBuffer = 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)
{ {
if (_disposed) if (_disposed)
@@ -366,23 +392,46 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
throw new ObjectDisposedException(nameof(DivisorScanSession)); throw new ObjectDisposedException(nameof(DivisorScanSession));
} }
int offset = 0; int primesLength = primes.Length;
if (primesLength == 0)
{
return;
}
MersenneNumberDivisorByDivisorGpuTester owner = _owner; MersenneNumberDivisorByDivisorGpuTester owner = _owner;
int gpuBatchSize = owner._gpuBatchSize; int gpuBatchSize = owner._gpuBatchSize;
bool useCycles = owner._useDivisorCycles; EnsureCapacity(gpuBatchSize);
Accelerator accelerator = _accelerator;
MemoryBuffer1D<byte, Stride1D.Dense> hitsBuffer = _hitsBuffer;
MemoryBuffer1D<ulong, Stride1D.Dense> primesBuffer = _primesBuffer;
ArrayView1D<ulong, Stride1D.Dense> primesBufferView = primesBuffer.View;
ArrayView1D<byte, Stride1D.Dense> hitsBufferView = hitsBuffer.View;
Action<Index1D, MontgomeryDivisorData, byte, ArrayView<ulong>, ArrayView<byte>> kernel = _kernel;
ArrayView1D<ulong, Stride1D.Dense> primesView;
ArrayView1D<byte, Stride1D.Dense> hitsView;
byte useCyclesFlag = useCycles ? (byte)1 : (byte)0;
MontgomeryDivisorData divisorData = CreateMontgomeryDivisorData(divisor); 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<Index1D, MontgomeryDivisorData, ArrayView<ulong>, ArrayView<ulong>> kernel = _kernel;
ArrayView1D<ulong, Stride1D.Dense> exponentsView = _exponentsBuffer.View;
ArrayView1D<ulong, Stride1D.Dense> resultsView = _resultsBuffer.View;
int offset = 0;
int batchSize;
ReadOnlySpan<ulong> primesSlice; ReadOnlySpan<ulong> primesSlice;
Span<byte> hitsSlice; Span<byte> hitsSlice;
int batchSize, primesLength = primes.Length; Span<ulong> hostSpan = _hostBuffer.AsSpan();
Span<int> positionSpan = _positionBuffer.AsSpan();
bool hasResidue = false;
ulong lastPrimeWithResidue = 0UL;
ulong lastResidue = 0UL;
while (offset < primesLength) while (offset < primesLength)
{ {
@@ -390,14 +439,60 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
primesSlice = primes.Slice(offset, batchSize); primesSlice = primes.Slice(offset, batchSize);
hitsSlice = hits.Slice(offset, batchSize); hitsSlice = hits.Slice(offset, batchSize);
primesView = primesBufferView.SubView(0, batchSize); int computeCount = 0;
hitsView = hitsBufferView.SubView(0, batchSize); ulong deltaBasePrime = lastPrimeWithResidue;
bool deltaHasResidue = hasResidue;
primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primesSlice), batchSize); for (int i = 0; i < batchSize; i++)
kernel(batchSize, divisorData, useCyclesFlag, primesView, hitsView); {
ulong prime = primesSlice[i];
if (cycleEnabled && prime % cycle != 0UL)
{
hitsSlice[i] = 0;
continue;
}
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<ulong> exponentSlice = hostSpan.Slice(0, computeCount);
ArrayView1D<ulong, Stride1D.Dense> exponentView = exponentsView.SubView(0, computeCount);
exponentView.CopyFromCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount);
ArrayView1D<ulong, Stride1D.Dense> resultView = resultsView.SubView(0, computeCount);
kernel(computeCount, divisorData, exponentView, resultView);
accelerator.Synchronize(); accelerator.Synchronize();
resultView.CopyToCPU(ref MemoryMarshal.GetReference(exponentSlice), computeCount);
hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSlice), batchSize); 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; offset += batchSize;
} }
@@ -413,6 +508,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
_disposed = true; _disposed = true;
_owner.ReturnSession(this); _owner.ReturnSession(this);
} }
} }
internal void ReturnSession(DivisorScanSession session) internal void ReturnSession(DivisorScanSession session)
@@ -420,26 +516,17 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester
_sessionPool.Add(session); _sessionPool.Add(session);
} }
private static void CheckKernelByPrime(Index1D index, MontgomeryDivisorData divisor, byte useCycleChecks, ArrayView<ulong> primes, ArrayView<byte> hits) private static void ComputePrimeDeltaKernel(Index1D index, MontgomeryDivisorData divisor, ArrayView<ulong> deltas, ArrayView<ulong> results)
{ {
ulong modulus = divisor.Modulus; ulong modulus = divisor.Modulus;
if (modulus <= 1UL || (modulus & 1UL) == 0UL) if (modulus <= 1UL || (modulus & 1UL) == 0UL)
{ {
hits[index] = 0; results[index] = 0UL;
return; return;
} }
if (useCycleChecks != 0) ulong exponent = deltas[index];
{ results[index] = exponent.Pow2MontgomeryMod(divisor);
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;
} }
private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus) private static MontgomeryDivisorData CreateMontgomeryDivisorData(ulong modulus)
PerfectNumbers.Core/ULongExtensions.cs
+3 -6
View File
@@ -169,12 +169,9 @@ public static class ULongExtensions
ulong w1 = xLow * yHigh; ulong w1 = xLow * yHigh;
ulong w2 = xHigh * yLow; ulong w2 = xHigh * yLow;
return (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32) + ulong result = (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32);
( result += (((xLow * yLow) >> 32) + (uint)w1 + (uint)w2) >> 32;
((xLow * yLow) >> 32) + return result;
(uint)w1 +
(uint)w2
) >> 32;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]