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

Merge branch 'main' into codex/analyze-evenperfectbitscanner-execution-with-mersenne

Browse Source
This commit is contained in:
Marek Lemańczyk authored and GitHub committed 2025-09-21 13:05:22 +02:00
2 parents f91f526121 458cb77227
commit e7a3713db8
3 files changed
+287 -422
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

EvenPerfectBitScanner/Program.cs
+1 -3
View File
@@ -16,8 +16,6 @@ internal static class Program
private static ThreadLocal<MersenneNumberTester> MersenneTesters = null!; private static ThreadLocal<MersenneNumberTester> MersenneTesters = null!;
private static ThreadLocal<ModResidueTracker> PResidue = null!; // p mod d tracker (per-thread) private static ThreadLocal<ModResidueTracker> PResidue = null!; // p mod d tracker (per-thread)
private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP;
private const int ConsoleInterval = 100_000;
private const int WriteBatchSize = 100; private const int WriteBatchSize = 100;
private const string DefaultCyclesPath = "divisor_cycles.bin"; private const string DefaultCyclesPath = "divisor_cycles.bin";
private static string ResultsFileName = "even_perfect_bit_scan_results.csv"; private static string ResultsFileName = "even_perfect_bit_scan_results.csv";
@@ -1082,7 +1080,7 @@ internal static class Program
lock (Sync) lock (Sync)
{ {
if (_consoleCounter >= ConsoleInterval) if (_consoleCounter >= PerfectNumberConstants.ConsoleInterval)
{ {
printToConsole = true; printToConsole = true;
primeFlag = passedAllTests && !lastWasComposite && _primeFoundAfterInit; primeFlag = passedAllTests && !lastWasComposite && _primeFoundAfterInit;
PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs
+283 -417
View File
@@ -10,422 +10,288 @@ namespace PerfectNumbers.Core.Gpu;
public sealed class MersenneNumberDivisorByDivisorGpuTester public sealed class MersenneNumberDivisorByDivisorGpuTester
{ {
private int _gpuBatchSize = GpuConstants.ScanBatchSize; private int _gpuBatchSize = GpuConstants.ScanBatchSize;
private readonly object _sync = new(); private readonly object _sync = new();
private ulong _divisorLimit; private ulong _divisorLimit;
private bool _isConfigured; private bool _isConfigured;
private ulong _lastStatusDivisor; private ulong _lastStatusDivisor;
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _kernelCache = new(); private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _kernelCache = new();
private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>> _kernelByPrimeCache = new();
private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> GetKernel(Accelerator accelerator) =>
private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> GetKernel(Accelerator accelerator) => _kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
_kernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(CheckKernel));
public int GpuBatchSize
private Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> GetKernelByPrime(Accelerator accelerator) => {
_kernelByPrimeCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>>(CheckKernelByPrime)); get => _gpuBatchSize;
set => _gpuBatchSize = Math.Max(1, value);
public int GpuBatchSize }
{
get => _gpuBatchSize; public void ConfigureFromMaxPrime(ulong maxPrime)
set => _gpuBatchSize = Math.Max(1, value); {
} lock (_sync)
{
public void ConfigureFromMaxPrime(ulong maxPrime) _divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime);
{ _lastStatusDivisor = 0UL;
lock (_sync) _isConfigured = true;
{ }
_divisorLimit = ComputeDivisorLimitFromMaxPrime(maxPrime); }
_lastStatusDivisor = 0UL;
_isConfigured = true; public bool IsPrime(ulong prime, out bool divisorsExhausted)
} {
} lock (_sync)
{
public ulong DivisorLimit if (!_isConfigured)
{ {
get throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
{ }
lock (_sync)
{ ulong allowedMax = ComputeAllowedMaxDivisor(prime);
if (!_isConfigured) if (allowedMax < 3UL)
{ {
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); divisorsExhausted = true;
} return true;
}
return _divisorLimit;
} bool composite = CheckDivisors(prime, allowedMax, out ulong lastProcessed, out bool coveredRange);
}
} if (lastProcessed != 0UL)
{
public ulong GetAllowedMaxDivisor(ulong prime) ReportStatus(lastProcessed);
{ }
lock (_sync)
{ if (composite)
if (!_isConfigured) {
{ divisorsExhausted = true;
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); return false;
} }
return ComputeAllowedMaxDivisor(prime); divisorsExhausted = coveredRange;
} return true;
} }
}
public DivisorScanSession CreateDivisorSession()
{ private bool CheckDivisors(ulong prime, ulong allowedMax, out ulong lastProcessed, out bool coveredRange)
lock (_sync) {
{ allowedMax = Math.Min(allowedMax, _divisorLimit);
if (!_isConfigured) lastProcessed = 0UL;
{
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); if (allowedMax < 3UL)
} {
coveredRange = true;
return new DivisorScanSession(this); return false;
} }
}
var gpuLease = GpuContextPool.RentPreferred(preferCpu: false);
public bool IsPrime(ulong prime, out bool divisorsExhausted) var accelerator = gpuLease.Accelerator;
{ var kernel = GetKernel(accelerator);
lock (_sync)
{ int batchCapacity = _gpuBatchSize;
if (!_isConfigured) var divisorsBuffer = accelerator.Allocate1D<ulong>(batchCapacity);
{ var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity);
throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester.");
} ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity);
ulong allowedMax = ComputeAllowedMaxDivisor(prime);
if (allowedMax < 3UL) bool composite = false;
{ bool processedAll = false;
divisorsExhausted = true; ulong divisor = 3UL;
return true;
} try
{
bool composite = CheckDivisors(prime, allowedMax, out ulong lastProcessed, out bool coveredRange); while (divisor <= allowedMax)
{
if (!composite && lastProcessed != 0UL) int batchSize = 0;
{ bool reachedEndInBatch = false;
ReportStatus(lastProcessed);
} while (batchSize < batchCapacity && divisor <= allowedMax)
{
if (composite) ulong current = divisor;
{ divisors[batchSize++] = current;
divisorsExhausted = true;
return false; ulong next = current + 2UL;
} if (next <= current)
{
divisorsExhausted = coveredRange; reachedEndInBatch = true;
return true; break;
} }
}
if (next > allowedMax)
private bool CheckDivisors(ulong prime, ulong allowedMax, out ulong lastProcessed, out bool coveredRange) {
{ divisor = next;
allowedMax = Math.Min(allowedMax, _divisorLimit); reachedEndInBatch = true;
lastProcessed = 0UL; break;
}
if (allowedMax < 3UL)
{ divisor = next;
coveredRange = true; }
return false;
} if (batchSize == 0)
{
var gpuLease = GpuContextPool.RentPreferred(preferCpu: false); break;
var accelerator = gpuLease.Accelerator; }
var kernel = GetKernel(accelerator);
var divisorView = divisorsBuffer.View.SubView(0, batchSize);
int batchCapacity = _gpuBatchSize; var hitsView = hitsBuffer.View.SubView(0, batchSize);
var divisorsBuffer = accelerator.Allocate1D<ulong>(batchCapacity);
var hitsBuffer = accelerator.Allocate1D<byte>(batchCapacity); Span<ulong> divisorSpan = divisors.AsSpan(0, batchSize);
Span<byte> hitsSpan = hits.AsSpan(0, batchSize);
ulong[] divisors = ArrayPool<ulong>.Shared.Rent(batchCapacity);
byte[] hits = ArrayPool<byte>.Shared.Rent(batchCapacity); divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorSpan), batchSize);
kernel(batchSize, prime, divisorView, hitsView);
bool composite = false; accelerator.Synchronize();
bool processedAll = false;
ulong divisor = 3UL; hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize);
try for (int i = 0; i < batchSize; i++)
{ {
while (divisor <= allowedMax) lastProcessed = divisorSpan[i];
{ if (hitsSpan[i] != 0)
int batchSize = 0; {
bool reachedEndInBatch = false; composite = true;
break;
while (batchSize < batchCapacity && divisor <= allowedMax) }
{ }
ulong current = divisor;
divisors[batchSize++] = current; if (composite)
{
ulong next = current + 2UL; break;
if (next <= current) }
{
reachedEndInBatch = true; if (reachedEndInBatch)
break; {
} processedAll = true;
break;
if (next > allowedMax) }
{ }
divisor = next;
reachedEndInBatch = true; coveredRange = composite || processedAll || divisor > allowedMax;
break; }
} finally
{
divisor = next; ArrayPool<byte>.Shared.Return(hits, clearArray: true);
} ArrayPool<ulong>.Shared.Return(divisors, clearArray: true);
hitsBuffer.Dispose();
if (batchSize == 0) divisorsBuffer.Dispose();
{ gpuLease.Dispose();
break; }
}
return composite;
var divisorView = divisorsBuffer.View.SubView(0, batchSize); }
var hitsView = hitsBuffer.View.SubView(0, batchSize);
private void ReportStatus(ulong divisor)
Span<ulong> divisorSpan = divisors.AsSpan(0, batchSize); {
Span<byte> hitsSpan = hits.AsSpan(0, batchSize); if (++_lastStatusDivisor == PerfectNumberConstants.ConsoleInterval)
{
divisorView.CopyFromCPU(ref MemoryMarshal.GetReference(divisorSpan), batchSize); Console.WriteLine($"...processed by-divisor candidate = {divisor}");
kernel(batchSize, prime, divisorView, hitsView); _lastStatusDivisor = 0;
accelerator.Synchronize(); }
}
hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitsSpan), batchSize);
private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime)
for (int i = 0; i < batchSize; i++) {
{ if (maxPrime <= 1UL)
lastProcessed = divisorSpan[i]; {
if (hitsSpan[i] != 0) return 0UL;
{ }
composite = true;
break; if (maxPrime - 1UL >= 64UL)
} {
} return ulong.MaxValue;
}
if (composite)
{ return (1UL << (int)(maxPrime - 1UL)) - 1UL;
break; }
}
private ulong ComputeAllowedMaxDivisor(ulong prime)
if (reachedEndInBatch) {
{ if (prime <= 1UL)
processedAll = true; {
break; return 0UL;
} }
}
if (prime - 1UL >= 64UL)
coveredRange = composite || processedAll || divisor > allowedMax; {
} return _divisorLimit;
finally }
{
ArrayPool<byte>.Shared.Return(hits, clearArray: true); ulong maxByPrime = (1UL << (int)(prime - 1UL)) - 1UL;
ArrayPool<ulong>.Shared.Return(divisors, clearArray: true); return Math.Min(maxByPrime, _divisorLimit);
hitsBuffer.Dispose(); }
divisorsBuffer.Dispose();
gpuLease.Dispose(); private static void CheckKernel(Index1D index, ulong prime, ArrayView<ulong> divisors, ArrayView<byte> hits)
} {
ulong divisor = divisors[index];
return composite; if (divisor <= 1UL || (divisor & 1UL) == 0UL)
} {
hits[index] = 0;
public sealed class DivisorScanSession : IDisposable return;
{ }
private readonly MersenneNumberDivisorByDivisorGpuTester _owner;
private readonly GpuContextPool.GpuContextLease _lease; ulong pow = Pow2Mod(prime, divisor);
private readonly Accelerator _accelerator; hits[index] = pow == 1UL ? (byte)1 : (byte)0;
private readonly Action<Index1D, ulong, ArrayView<ulong>, ArrayView<byte>> _kernel; }
private readonly MemoryBuffer1D<ulong, Stride1D.Dense> _primesBuffer;
private readonly MemoryBuffer1D<byte, Stride1D.Dense> _hitsBuffer; private static ulong Pow2Mod(ulong exponent, ulong modulus)
private readonly ulong[] _primesHost; {
private readonly byte[] _hitsHost; if (modulus <= 1UL)
private bool _disposed; {
return 0UL;
internal DivisorScanSession(MersenneNumberDivisorByDivisorGpuTester owner) }
{
_owner = owner; ulong result = 1UL % modulus;
_lease = GpuContextPool.RentPreferred(preferCpu: false); ulong baseVal = 2UL % modulus;
_accelerator = _lease.Accelerator; ulong exp = exponent;
_kernel = owner.GetKernelByPrime(_accelerator); while (exp > 0UL)
_primesBuffer = _accelerator.Allocate1D<ulong>(owner._gpuBatchSize); {
_hitsBuffer = _accelerator.Allocate1D<byte>(owner._gpuBatchSize); if ((exp & 1UL) != 0UL)
_primesHost = ArrayPool<ulong>.Shared.Rent(owner._gpuBatchSize); {
_hitsHost = ArrayPool<byte>.Shared.Rent(owner._gpuBatchSize); result = MulMod(result, baseVal, modulus);
} }
public void CheckDivisor(ulong divisor, ReadOnlySpan<ulong> primes, Span<byte> hits) exp >>= 1;
{ if (exp == 0UL)
if (_disposed) {
{ break;
throw new ObjectDisposedException(nameof(DivisorScanSession)); }
}
baseVal = MulMod(baseVal, baseVal, modulus);
if (primes.Length == 0) }
{
return; return result;
} }
int offset = 0; private static ulong MulMod(ulong a, ulong b, ulong modulus)
while (offset < primes.Length) {
{ if (modulus == 0UL)
int batchSize = Math.Min(_owner._gpuBatchSize, primes.Length - offset); {
Span<ulong> primeSpan = _primesHost.AsSpan(0, batchSize); return 0UL;
primes.Slice(offset, batchSize).CopyTo(primeSpan); }
var primesView = _primesBuffer.View.SubView(0, batchSize); ulong result = 0UL;
var hitsView = _hitsBuffer.View.SubView(0, batchSize); ulong x = a % modulus;
ulong y = b;
primesView.CopyFromCPU(ref MemoryMarshal.GetReference(primeSpan), batchSize); while (y > 0UL)
_kernel(batchSize, divisor, primesView, hitsView); {
_accelerator.Synchronize(); if ((y & 1UL) != 0UL)
{
Span<byte> hitSpan = _hitsHost.AsSpan(0, batchSize); result += x;
hitsView.CopyToCPU(ref MemoryMarshal.GetReference(hitSpan), batchSize); if (result >= modulus)
hitSpan.CopyTo(hits.Slice(offset, batchSize)); {
offset += batchSize; result -= modulus;
} }
} }
public void Dispose() x <<= 1;
{ if (x >= modulus)
if (_disposed) {
{ x -= modulus;
return; }
}
y >>= 1;
_disposed = true; }
ArrayPool<ulong>.Shared.Return(_primesHost, clearArray: true);
ArrayPool<byte>.Shared.Return(_hitsHost, clearArray: true); return result;
_primesBuffer.Dispose(); }
_hitsBuffer.Dispose();
_lease.Dispose();
}
}
private void ReportStatus(ulong divisor)
{
if (divisor <= _lastStatusDivisor)
{
return;
}
_lastStatusDivisor = divisor;
Console.WriteLine($"...processed by-divisor candidate = {divisor}");
}
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 ulong ComputeAllowedMaxDivisor(ulong prime)
{
if (prime <= 1UL)
{
return 0UL;
}
if (prime - 1UL >= 64UL)
{
return _divisorLimit;
}
ulong maxByPrime = (1UL << (int)(prime - 1UL)) - 1UL;
return Math.Min(maxByPrime, _divisorLimit);
}
private static void CheckKernel(Index1D index, ulong prime, ArrayView<ulong> divisors, ArrayView<byte> hits)
{
ulong divisor = divisors[index];
if (divisor <= 1UL || (divisor & 1UL) == 0UL)
{
hits[index] = 0;
return;
}
ulong pow = Pow2Mod(prime, divisor);
hits[index] = pow == 1UL ? (byte)1 : (byte)0;
}
private static void CheckKernelByPrime(Index1D index, ulong divisor, ArrayView<ulong> primes, ArrayView<byte> hits)
{
ulong prime = primes[index];
if (divisor <= 1UL || (divisor & 1UL) == 0UL)
{
hits[index] = 0;
return;
}
ulong pow = Pow2Mod(prime, divisor);
hits[index] = pow == 1UL ? (byte)1 : (byte)0;
}
private static ulong Pow2Mod(ulong exponent, ulong modulus)
{
if (modulus <= 1UL)
{
return 0UL;
}
ulong result = 1UL % modulus;
ulong baseVal = 2UL % modulus;
ulong exp = exponent;
while (exp > 0UL)
{
if ((exp & 1UL) != 0UL)
{
result = MulMod(result, baseVal, modulus);
}
exp >>= 1;
if (exp == 0UL)
{
break;
}
baseVal = MulMod(baseVal, baseVal, modulus);
}
return result;
}
private static ulong MulMod(ulong a, ulong b, ulong modulus)
{
if (modulus == 0UL)
{
return 0UL;
}
ulong result = 0UL;
ulong x = a % modulus;
ulong y = b;
while (y > 0UL)
{
if ((y & 1UL) != 0UL)
{
result += x;
if (result >= modulus)
{
result -= modulus;
}
}
x <<= 1;
if (x >= modulus)
{
x -= modulus;
}
y >>= 1;
}
return result;
}
} }
PerfectNumbers.Core/PerfectNumberConstants.cs
+3 -2
View File
@@ -3,8 +3,9 @@ namespace PerfectNumbers.Core;
public static class PerfectNumberConstants public static class PerfectNumberConstants
{ {
public const ulong BiggestKnownEvenPerfectP = 136279841UL; public const ulong BiggestKnownEvenPerfectP = 136279841UL;
public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int ConsoleInterval = 100_000;
public const int MaxQForDivisorCycles = 4_000_000;
public const ulong ExtraDivisorCycleSearchLimit = 64UL; public const ulong ExtraDivisorCycleSearchLimit = 64UL;
public const int MaxQForDivisorCycles = 4_000_000;
public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000;
} }