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

Tiny cleanup

Browse Source
This commit is contained in:
Marek Lemańczyk committed 2025-09-20 10:48:12 +02:00
1 parent d78503eade
commit fb7734528f
1 file changed
+229 -229
Download Patch File Download Diff File Expand all files Collapse all files
PerfectNumbers.Core/MersenneNumberTester.cs
+229 -229
View File
@@ -11,63 +11,63 @@ namespace PerfectNumbers.Core;
public enum GpuKernelType public enum GpuKernelType
{ {
Incremental, Incremental,
Pow2Mod, Pow2Mod,
} }
public sealed class MersenneNumberTester( public sealed class MersenneNumberTester(
bool useIncremental = true, bool useIncremental = true,
bool useOrderCache = false, bool useOrderCache = false,
GpuKernelType kernelType = GpuKernelType.Incremental, GpuKernelType kernelType = GpuKernelType.Incremental,
bool useModuloWorkaround = false, bool useModuloWorkaround = false,
bool useOrder = false, bool useOrder = false,
bool useGpuLucas = true, bool useGpuLucas = true,
bool useGpuScan = true, bool useGpuScan = true,
bool useGpuOrder = true, bool useGpuOrder = true,
bool useResidue = true, bool useResidue = true,
UInt128? maxK = null, UInt128? maxK = null,
UInt128? residueDivisorSets = null) UInt128? residueDivisorSets = null)
{ {
private readonly bool _useResidue = useResidue; private readonly bool _useResidue = useResidue;
private readonly bool _useIncremental = useIncremental && !useResidue; private readonly bool _useIncremental = useIncremental && !useResidue;
private readonly UInt128 _maxK = maxK ?? (UInt128)5_000_000UL; private readonly UInt128 _maxK = maxK ?? (UInt128)5_000_000UL;
private readonly UInt128 _residueDivisorSetCount = residueDivisorSets ?? (UInt128)PerfectNumberConstants.ExtraDivisorCycleSearchLimit; private readonly UInt128 _residueDivisorSetCount = residueDivisorSets ?? (UInt128)PerfectNumberConstants.ExtraDivisorCycleSearchLimit;
private readonly GpuKernelType _kernelType = kernelType; private readonly GpuKernelType _kernelType = kernelType;
private readonly bool _useModuloWorkaround = useModuloWorkaround; private readonly bool _useModuloWorkaround = useModuloWorkaround;
private readonly bool _useGpuLucas = useGpuLucas; private readonly bool _useGpuLucas = useGpuLucas;
private readonly bool _useGpuScan = useGpuScan; // device for pow2mod/incremental scanning private readonly bool _useGpuScan = useGpuScan; // device for pow2mod/incremental scanning
private readonly bool _useGpuOrder = useGpuOrder; // device for order computations private readonly bool _useGpuOrder = useGpuOrder; // device for order computations
private static readonly ConcurrentDictionary<UInt128, ulong> OrderCache = new(); private static readonly ConcurrentDictionary<UInt128, ulong> OrderCache = new();
private readonly MersenneNumberIncrementalGpuTester _incrementalGpuTester = new(kernelType, useGpuOrder); private readonly MersenneNumberIncrementalGpuTester _incrementalGpuTester = new(kernelType, useGpuOrder);
private readonly MersenneNumberIncrementalCpuTester _incrementalCpuTester = new(kernelType); private readonly MersenneNumberIncrementalCpuTester _incrementalCpuTester = new(kernelType);
private readonly MersenneNumberOrderGpuTester _orderGpuTester = new(kernelType, useGpuOrder); private readonly MersenneNumberOrderGpuTester _orderGpuTester = new(kernelType, useGpuOrder);
private readonly MersenneNumberOrderCpuTester _orderCpuTester = new(kernelType); private readonly MersenneNumberOrderCpuTester _orderCpuTester = new(kernelType);
private readonly MersenneNumberLucasLehmerGpuTester _lucasLehmerGpuTester = new(); private readonly MersenneNumberLucasLehmerGpuTester _lucasLehmerGpuTester = new();
private readonly MersenneNumberLucasLehmerCpuTester _lucasLehmerCpuTester = new(); private readonly MersenneNumberLucasLehmerCpuTester _lucasLehmerCpuTester = new();
private readonly MersenneNumberResidueGpuTester? _residueGpuTester = useResidue ? new(useGpuOrder) : null; private readonly MersenneNumberResidueGpuTester? _residueGpuTester = useResidue ? new(useGpuOrder) : null;
private readonly MersenneNumberResidueCpuTester? _residueCpuTester = useResidue ? new() : null; private readonly MersenneNumberResidueCpuTester? _residueCpuTester = useResidue ? new() : null;
// TODO: Ensure Program passes useResidue correctly and that residue-vs-incremental-vs-LL selection respects CLI flags. // TODO: Ensure Program passes useResidue correctly and that residue-vs-incremental-vs-LL selection respects CLI flags.
// Also consider injecting a shared cycles cache (MersenneDivisorCycles.Shared) to both CPU and GPU testers. // Also consider injecting a shared cycles cache (MersenneDivisorCycles.Shared) to both CPU and GPU testers.
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static UInt128 DivideRoundUp(UInt128 value, UInt128 divisor) private static UInt128 DivideRoundUp(UInt128 value, UInt128 divisor)
{ {
if (divisor == UInt128.Zero) if (divisor == UInt128.Zero)
{ {
return UInt128.Zero; return UInt128.Zero;
} }
UInt128 result = value / divisor; UInt128 result = value / divisor;
if (result * divisor == value) if (result * divisor == value)
{ {
return result; return result;
} }
return result + UInt128.One; return result + UInt128.One;
} }
public void WarmUpOrders(ulong exponent, ulong limit = 5_000_000UL) public void WarmUpOrders(ulong exponent, ulong limit = 5_000_000UL)
{ {
bool cacheEnabled = useOrderCache; bool cacheEnabled = useOrderCache;
if (!cacheEnabled) if (!cacheEnabled)
@@ -158,168 +158,168 @@ public sealed class MersenneNumberTester(
return; return;
} }
UInt128[] qs = qsBuffer; UInt128[] qs = qsBuffer;
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var accelerator = gpuLease.Accelerator; var accelerator = gpuLease.Accelerator;
var orderKernel = gpuLease.OrderKernel!; var orderKernel = gpuLease.OrderKernel!;
// Guard long-running kernels by chunking the warm-up across batches. // Guard long-running kernels by chunking the warm-up across batches.
// Reuse the same ScanBatchSize knob used for scanning. // Reuse the same ScanBatchSize knob used for scanning.
int batchSize = GpuConstants.ScanBatchSize; int batchSize = GpuConstants.ScanBatchSize;
ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL; ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL;
int offset = 0; int offset = 0;
while (offset < idx) while (offset < idx)
{ {
int count = Math.Min(batchSize, idx - offset); int count = Math.Min(batchSize, idx - offset);
var qBuffer = accelerator.Allocate1D<PerfectNumbers.Core.Gpu.GpuUInt128>(count); var qBuffer = accelerator.Allocate1D<PerfectNumbers.Core.Gpu.GpuUInt128>(count);
var orderBuffer = accelerator.Allocate1D<ulong>(count); var orderBuffer = accelerator.Allocate1D<ulong>(count);
// Convert to device-friendly GpuUInt128 on the fly // Convert to device-friendly GpuUInt128 on the fly
var tmp = System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Rent(count); var tmp = System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Rent(count);
for (int i = 0; i < count; i++) for (int i = 0; i < count; i++)
{ {
tmp[i] = (PerfectNumbers.Core.Gpu.GpuUInt128)qs[offset + i]; tmp[i] = (PerfectNumbers.Core.Gpu.GpuUInt128)qs[offset + i];
} }
// Copy only the populated elements to the device buffer to avoid overrun on pooled arrays. // Copy only the populated elements to the device buffer to avoid overrun on pooled arrays.
qBuffer.View.CopyFromCPU(ref tmp[0], count); qBuffer.View.CopyFromCPU(ref tmp[0], count);
System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Return(tmp); System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Return(tmp);
orderKernel(count, exponent, divMul, qBuffer.View, orderBuffer.View); orderKernel(count, exponent, divMul, qBuffer.View, orderBuffer.View);
accelerator.Synchronize(); accelerator.Synchronize();
ulong[] orders = orderBuffer.GetAsArray1D(); ulong[] orders = orderBuffer.GetAsArray1D();
for (int i = 0; i < count; i++) for (int i = 0; i < count; i++)
{ {
ulong order = orders[i]; ulong order = orders[i];
if (order == 0UL) if (order == 0UL)
{ {
order = qs[offset + i].CalculateOrder(); order = qs[offset + i].CalculateOrder();
} }
if (useOrderCache) if (useOrderCache)
{ {
OrderCache[qs[offset + i]] = order; OrderCache[qs[offset + i]] = order;
} }
} }
offset += count; offset += count;
orderBuffer.Dispose(); orderBuffer.Dispose();
qBuffer.Dispose(); qBuffer.Dispose();
} }
gpuLease.Dispose(); gpuLease.Dispose();
ArrayPool<UInt128>.Shared.Return(qs); ArrayPool<UInt128>.Shared.Return(qs);
} }
public bool IsMersennePrime(ulong exponent) => IsMersennePrime(exponent, out _); public bool IsMersennePrime(ulong exponent) => IsMersennePrime(exponent, out _);
public bool IsMersennePrime(ulong exponent, out bool divisorsExhausted) public bool IsMersennePrime(ulong exponent, out bool divisorsExhausted)
{ {
divisorsExhausted = false; divisorsExhausted = false;
// Safe early rejections using known orders for tiny primes: // Safe early rejections using known orders for tiny primes:
// 7 | M_p iff 3 | p. Avoid rejecting p == 3 where M_3 == 7 is prime. // 7 | M_p iff 3 | p. Avoid rejecting p == 3 where M_3 == 7 is prime.
if ((exponent % 3UL) == 0UL && exponent != 3UL) if ((exponent % 3UL) == 0UL && exponent != 3UL)
{ {
divisorsExhausted = true; divisorsExhausted = true;
return false; return false;
} }
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne()) if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
{ {
divisorsExhausted = true; divisorsExhausted = true;
return false; return false;
} }
// Skip even exponents: for even p > 2, M_p is composite; here p is large. // Skip even exponents: for even p > 2, M_p is composite; here p is large.
if ((exponent & 1UL) == 0UL) if ((exponent & 1UL) == 0UL)
{ {
divisorsExhausted = true; divisorsExhausted = true;
return false; return false;
} }
bool prePrime = true; bool prePrime = true;
UInt128 twoP = (UInt128)exponent << 1; // 2 * p UInt128 twoP = (UInt128)exponent << 1; // 2 * p
// UInt128 maxDivisor = new(ulong.MaxValue, ulong.MaxValue); // UInt128 maxDivisor = new(ulong.MaxValue, ulong.MaxValue);
UInt128 maxK = UInt128Numbers.Two.Pow(exponent) / 2 + 1; UInt128 maxK = UInt128Numbers.Two.Pow(exponent) / 2 + 1;
bool lastIsSeven = (exponent & 3UL) == 3UL; bool lastIsSeven = (exponent & 3UL) == 3UL;
// In residue mode the residue scan is the final primality check. // In residue mode the residue scan is the final primality check.
if (_useResidue) if (_useResidue)
{ {
if (_maxK < maxK) if (_maxK < maxK)
{ {
maxK = _maxK; maxK = _maxK;
} }
if (_maxK == UInt128.Zero || _residueDivisorSetCount == UInt128.Zero) if (_maxK == UInt128.Zero || _residueDivisorSetCount == UInt128.Zero)
{ {
return prePrime; return prePrime;
} }
UInt128 totalLimit = _maxK; UInt128 totalLimit = _maxK;
if (totalLimit > maxK) if (totalLimit > maxK)
{ {
totalLimit = maxK; totalLimit = maxK;
} }
if (totalLimit == UInt128.Zero) if (totalLimit == UInt128.Zero)
{ {
return true; return true;
} }
UInt128 requestedSets = _residueDivisorSetCount; UInt128 requestedSets = _residueDivisorSetCount;
if (_useGpuScan) if (_useGpuScan)
{ {
UInt128 maxUsefulSets = DivideRoundUp(totalLimit, (UInt128)GpuConstants.ScanBatchSize); UInt128 maxUsefulSets = DivideRoundUp(totalLimit, (UInt128)GpuConstants.ScanBatchSize);
if (maxUsefulSets == UInt128.Zero) if (maxUsefulSets == UInt128.Zero)
{ {
maxUsefulSets = UInt128.One; maxUsefulSets = UInt128.One;
} }
if (requestedSets == UInt128.Zero || requestedSets > maxUsefulSets) if (requestedSets == UInt128.Zero || requestedSets > maxUsefulSets)
{ {
requestedSets = maxUsefulSets; requestedSets = maxUsefulSets;
} }
} }
if (requestedSets == UInt128.Zero) if (requestedSets == UInt128.Zero)
{ {
requestedSets = UInt128.One; requestedSets = UInt128.One;
} }
UInt128 perSetLimit = DivideRoundUp(totalLimit, requestedSets); UInt128 perSetLimit = DivideRoundUp(totalLimit, requestedSets);
if (perSetLimit == UInt128.Zero) if (perSetLimit == UInt128.Zero)
{ {
perSetLimit = UInt128.One; perSetLimit = UInt128.One;
} }
UInt128 setsNeeded = DivideRoundUp(totalLimit, perSetLimit); UInt128 setsNeeded = DivideRoundUp(totalLimit, perSetLimit);
UInt128 effectiveSets = requestedSets; UInt128 effectiveSets = requestedSets;
if (setsNeeded < effectiveSets) if (setsNeeded < effectiveSets)
{ {
effectiveSets = setsNeeded; effectiveSets = setsNeeded;
} }
bool scanCompleted = false; bool scanCompleted = false;
if (_useGpuScan) if (_useGpuScan)
{ {
_residueGpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted); _residueGpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted);
} }
else else
{ {
_residueCpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted); _residueCpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted);
} }
if (!prePrime) if (!prePrime)
{ {
divisorsExhausted = true; divisorsExhausted = true;
return false; return false;
} }
divisorsExhausted = scanCompleted; divisorsExhausted = scanCompleted;
return true; return true;
} }
if (!_useIncremental) if (!_useIncremental)
{ {
@@ -337,64 +337,64 @@ public sealed class MersenneNumberTester(
_incrementalCpuTester.Scan(exponent, twoPPre, lastIsSevenPre, maxKPre, ref prePrime); _incrementalCpuTester.Scan(exponent, twoPPre, lastIsSevenPre, maxKPre, ref prePrime);
} }
if (!prePrime) if (!prePrime)
{ {
divisorsExhausted = true; divisorsExhausted = true;
return false; return false;
} }
bool lucasPrime = _useGpuLucas bool lucasPrime = _useGpuLucas
? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true) ? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true)
: _lucasLehmerCpuTester.IsPrime(exponent); : _lucasLehmerCpuTester.IsPrime(exponent);
divisorsExhausted = true; divisorsExhausted = true;
return lucasPrime; return lucasPrime;
} }
bool isPrime = true; bool isPrime = true;
if (useOrder) if (useOrder)
{ {
if (_useGpuOrder) if (_useGpuOrder)
{ {
_orderGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); _orderGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime);
} }
else else
{ {
_orderCpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); _orderCpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime);
} }
} }
else else
{ {
if (_useGpuScan) if (_useGpuScan)
{ {
_incrementalGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); _incrementalGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime);
} }
else else
{ {
_incrementalCpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); _incrementalCpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime);
} }
} }
divisorsExhausted = true; divisorsExhausted = true;
return isPrime; return isPrime;
} }
// Limit of k to scan in the pre-Lucas prescan // Limit of k to scan in the pre-Lucas prescan
public const int PreLucasPreScanK = 5_000_000; public const int PreLucasPreScanK = 5_000_000;
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static bool DividesMersenne(ulong exponent, ulong prime) private static bool DividesMersenne(ulong exponent, ulong prime)
{ {
return 2UL.ModPow64(exponent, prime) == 1UL; return 2UL.ModPow64(exponent, prime) == 1UL;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void LegacyOrderKernel(Index1D index, ulong exponent, ArrayView<UInt128> qs, ArrayView<ulong> orders) private static void LegacyOrderKernel(Index1D index, ulong exponent, ArrayView<UInt128> qs, ArrayView<ulong> orders)
{ {
UInt128 q = qs[index]; UInt128 q = qs[index];
UInt128 pow = exponent.PowMod(q); UInt128 pow = exponent.PowMod(q);
orders[index] = pow == 1UL ? exponent : 0UL; orders[index] = pow == 1UL ? exponent : 0UL;
} }
} }