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

Revert "Merge pull request #24 from mlemanczyk/codex/modify-evenperfectbitscanner-for-residue-mode"

Browse Source 
This reverts commit 2d024cc708, reversing
changes made to b1c29ad16e.
This commit is contained in:
Marek Lemańczyk committed 2025-09-21 02:25:18 +02:00
1 parent deb25796a8
commit 8a6fc0c6a5
7 files changed
+265 -551
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs
+136 -166
View File
@@ -5,32 +5,25 @@ public class MersenneNumberResidueCpuTester
private ModResidueTracker? _mersenneResidueTracker;
// CPU residue variant using tracker + unrolled residue updates (no method calls in the loop).
public void Scan(
ulong exponent,
UInt128 twoP,
bool lastIsSeven,
UInt128 perSetLimit,
UInt128 setCount,
UInt128 overallLimit,
ref bool isPrime,
ref bool divisorsExhausted)
{
if (setCount == UInt128.Zero || perSetLimit == UInt128.Zero || overallLimit == UInt128.Zero)
{
return;
}
public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
{
// Initialize/update Mersenne residue tracker and start a merge walk for ascending q divisors
// Initialize/update Mersenne residue tracker and start a merge walk for ascending q divisors
_mersenneResidueTracker ??= new(ResidueModel.Mersenne);
// For small divisors q <= 4,000,000 use cached cycle length to fast-path checks
// and avoid per-q powmods when possible when scanning lanes.
_mersenneResidueTracker.BeginMerge(exponent);
_mersenneResidueTracker ??= new(ResidueModel.Mersenne);
// For small divisors q <= 4,000,000 use cached cycle length to fast-path checks
// and avoid per-q powmods when possible when scanning lanes.
_mersenneResidueTracker.BeginMerge(exponent);
ulong step10 = (ulong)(twoP % 10UL) % 10UL;
ulong step8 = (ulong)(twoP % 8UL);
ulong step3 = (ulong)(twoP % 3UL);
ulong step5 = (ulong)(twoP % 5UL);
// Start at k = 1: q0 = 2*p*1 + 1
UInt128 k = UInt128.One;
UInt128 q = twoP + UInt128.One;
// Current residues
q.Mod10_8_5_3(out ulong r10, out ulong r8, out ulong r5, out ulong r3);
// Step residues for q += 2*p
ulong step10 = (ulong)(twoP % 10UL) % 10UL;
ulong step8 = (ulong)(twoP % 8UL);
ulong step3 = (ulong)(twoP % 3UL);
ulong step5 = (ulong)(twoP % 5UL);
// Allowed last-digit sets for q depending on last digit of M_p
// lastIsSeven == true => allow {7,9}
@@ -38,162 +31,139 @@ public class MersenneNumberResidueCpuTester
int allowMask = lastIsSeven ? CpuConstants.LastSevenMask10 : CpuConstants.LastOneMask10;
// Predeclare temps to avoid redeclaration overhead in the loop
UInt128 remaining;
int i, lanes;
// per-lane residues
ulong r10_0, r10_1, r10_2, r10_3, r8_0, r8_1, r8_2, r8_3, r3_0, r3_1, r3_2, r3_3, r5_0, r5_1, r5_2, r5_3;
// q lanes
UInt128 twoP2 = twoP + twoP,
twoP3 = twoP2 + twoP,
twoP4 = twoP3 + twoP,
qIncremental;
UInt128 remaining;
int i, lanes;
// per-lane residues
ulong r10_0, r10_1, r10_2, r10_3, r8_0, r8_1, r8_2, r8_3, r3_0, r3_1, r3_2, r3_3, r5_0, r5_1, r5_2, r5_3;
// q lanes
UInt128 twoP2 = twoP + twoP,
twoP3 = twoP2 + twoP,
twoP4 = twoP3 + twoP,
qIncremental;
bool localIsPrime, isDivider;
bool localIsPrime, isDivider;
ModResidueTracker tracker = _mersenneResidueTracker!;
UInt128 limitInclusive = overallLimit + UInt128.One;
for (UInt128 setIndex = UInt128.Zero; setIndex < setCount && Volatile.Read(ref isPrime); setIndex++)
{
UInt128 setOffset = perSetLimit * setIndex;
UInt128 k = setOffset + UInt128.One;
if (k >= limitInclusive)
{
break;
}
maxK++;
UInt128 setLimitExclusive = setOffset + perSetLimit + UInt128.One;
if (setLimitExclusive > limitInclusive)
{
setLimitExclusive = limitInclusive;
}
ModResidueTracker tracker = _mersenneResidueTracker!;
while (k < maxK && (localIsPrime = Volatile.Read(ref isPrime)))
{
remaining = maxK - k;
if (remaining == UInt128.Zero)
{
break;
}
UInt128 q = twoP * k + UInt128.One;
q.Mod10_8_5_3(out ulong r10, out ulong r8, out ulong r5, out ulong r3);
lanes = remaining > 4UL ? 4 : (int)remaining;
while (k < setLimitExclusive && (localIsPrime = Volatile.Read(ref isPrime)))
{
remaining = setLimitExclusive - k;
if (remaining == UInt128.Zero)
{
break;
}
// Compute lane residues incrementally (avoid multiplications, minimize reductions)
// r10
r10_0 = r10;
r10_1 = r10_0 + step10; if (r10_1 >= 20UL) r10_1 -= 20UL; if (r10_1 >= 10UL) r10_1 -= 10UL;
r10_2 = r10_1 + step10; if (r10_2 >= 20UL) r10_2 -= 20UL; if (r10_2 >= 10UL) r10_2 -= 10UL;
r10_3 = r10_2 + step10; if (r10_3 >= 20UL) r10_3 -= 20UL; if (r10_3 >= 10UL) r10_3 -= 10UL;
lanes = remaining > 4UL ? 4 : (int)remaining;
if (lanes <= 0)
{
break;
}
// r8
r8_0 = r8 & 7UL;
r8_1 = (r8_0 + step8) & 7UL;
r8_2 = (r8_1 + step8) & 7UL;
r8_3 = (r8_2 + step8) & 7UL;
// Compute lane residues incrementally (avoid multiplications, minimize reductions)
// r10
r10_0 = r10;
r10_1 = r10_0 + step10; if (r10_1 >= 20UL) r10_1 -= 20UL; if (r10_1 >= 10UL) r10_1 -= 10UL;
r10_2 = r10_1 + step10; if (r10_2 >= 20UL) r10_2 -= 20UL; if (r10_2 >= 10UL) r10_2 -= 10UL;
r10_3 = r10_2 + step10; if (r10_3 >= 20UL) r10_3 -= 20UL; if (r10_3 >= 10UL) r10_3 -= 10UL;
// r3 (pattern 0,1,2,0)
r3_0 = r3; if (r3_0 >= 3UL) r3_0 -= 3UL;
r3_1 = r3_0 + step3; if (r3_1 >= 3UL) r3_1 -= 3UL;
r3_2 = r3_1 + step3; if (r3_2 >= 3UL) r3_2 -= 3UL;
r3_3 = r3; if (r3_3 >= 3UL) r3_3 -= 3UL;
// r8
r8_0 = r8 & 7UL;
r8_1 = (r8_0 + step8) & 7UL;
r8_2 = (r8_1 + step8) & 7UL;
r8_3 = (r8_2 + step8) & 7UL;
// r5 (pattern 0,1,2,3)
r5_0 = r5; if (r5_0 >= 5UL) r5_0 -= 5UL;
r5_1 = r5_0 + step5; if (r5_1 >= 10UL) r5_1 -= 10UL; if (r5_1 >= 5UL) r5_1 -= 5UL;
r5_2 = r5_1 + step5; if (r5_2 >= 10UL) r5_2 -= 10UL; if (r5_2 >= 5UL) r5_2 -= 5UL;
r5_3 = r5_2 + step5; if (r5_3 >= 10UL) r5_3 -= 10UL; if (r5_3 >= 5UL) r5_3 -= 5UL;
// r3 (pattern 0,1,2,0)
r3_0 = r3; if (r3_0 >= 3UL) r3_0 -= 3UL;
r3_1 = r3_0 + step3; if (r3_1 >= 3UL) r3_1 -= 3UL;
r3_2 = r3_1 + step3; if (r3_2 >= 3UL) r3_2 -= 3UL;
r3_3 = r3; if (r3_3 >= 3UL) r3_3 -= 3UL;
// lane 0
if (lanes >= 1 &&
(((allowMask >> (int)r10_0) & 1) != 0 && (r8_0 == 1UL || r8_0 == 7UL) && r3_0 != 0UL && r5_0 != 0UL))
{
// cycle-based quick check for small q
if (localIsPrime && tracker.MergeOrAppend(exponent, q, out isDivider) && isDivider &&
(q <= ulong.MaxValue ? ((ulong)q).IsPrimeCandidate() : q.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
// r5 (pattern 0,1,2,3)
r5_0 = r5; if (r5_0 >= 5UL) r5_0 -= 5UL;
r5_1 = r5_0 + step5; if (r5_1 >= 10UL) r5_1 -= 10UL; if (r5_1 >= 5UL) r5_1 -= 5UL;
r5_2 = r5_1 + step5; if (r5_2 >= 10UL) r5_2 -= 10UL; if (r5_2 >= 5UL) r5_2 -= 5UL;
r5_3 = r5_2 + step5; if (r5_3 >= 10UL) r5_3 -= 10UL; if (r5_3 >= 5UL) r5_3 -= 5UL;
// lane 1
qIncremental = q + twoP;
if (lanes >= 2 &&
((allowMask >> (int)r10_1) & 1) != 0 && (r8_1 == 1UL || r8_1 == 7UL) && r3_1 != 0UL && r5_1 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
// lane 0
if (lanes >= 1 &&
(((allowMask >> (int)r10_0) & 1) != 0 && (r8_0 == 1UL || r8_0 == 7UL) && r3_0 != 0UL && r5_0 != 0UL))
{
if (localIsPrime && tracker.MergeOrAppend(exponent, q, out isDivider) && isDivider &&
(q <= ulong.MaxValue ? ((ulong)q).IsPrimeCandidate() : q.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
// lane 2
qIncremental += twoP;
if (lanes >= 3 &&
((allowMask >> (int)r10_2) & 1) != 0 && (r8_2 == 1UL || r8_2 == 7UL) && r3_2 != 0UL && r5_2 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
// lane 1
qIncremental = q + twoP;
if (lanes >= 2 &&
((allowMask >> (int)r10_1) & 1) != 0 && (r8_1 == 1UL || r8_1 == 7UL) && r3_1 != 0UL && r5_1 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
// lane 3
qIncremental += twoP;
if (lanes >= 4 &&
((allowMask >> (int)r10_3) & 1) != 0 && (r8_3 == 1UL || r8_3 == 7UL) && r3_3 != 0UL && r5_3 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
// lane 2
qIncremental += twoP;
if (lanes >= 3 &&
((allowMask >> (int)r10_2) & 1) != 0 && (r8_2 == 1UL || r8_2 == 7UL) && r3_2 != 0UL && r5_2 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
if (!localIsPrime)
{
Volatile.Write(ref isPrime, false);
return;
}
else if (!Volatile.Read(ref isPrime))
{
return;
}
// lane 3
qIncremental += twoP;
if (lanes >= 4 &&
((allowMask >> (int)r10_3) & 1) != 0 && (r8_3 == 1UL || r8_3 == 7UL) && r3_3 != 0UL && r5_3 != 0UL &&
localIsPrime)
{
if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
{
localIsPrime = false;
}
}
// advance base by processed lanes
k += (UInt128)lanes;
q += lanes switch
{
0 => 0UL,
1 => twoP,
2 => twoP2,
3 => twoP3,
4 => twoP4,
_ => throw new ArgumentException("Unsupproted value", nameof(twoP))
};
if (!localIsPrime)
{
Volatile.Write(ref isPrime, false);
divisorsExhausted = true;
return;
}
else if (!Volatile.Read(ref isPrime))
{
divisorsExhausted = true;
return;
}
// advance residues by 'lanes' steps (consistent with base advance)
for (i = 0; i < lanes; i++)
{
r10 += step10; if (r10 >= 20UL) r10 -= 20UL; if (r10 >= 10UL) r10 -= 10UL;
r8 = (r8 + step8) & 7UL;
r3 += step3; if (r3 >= 3UL) r3 -= 3UL;
r5 += step5; if (r5 >= 10UL) r5 -= 10UL; if (r5 >= 5UL) r5 -= 5UL;
}
}
// advance base by processed lanes
k += (UInt128)lanes;
q += lanes switch
{
0 => 0UL,
1 => twoP,
2 => twoP2,
3 => twoP3,
4 => twoP4,
_ => throw new ArgumentException("Unsupproted value", nameof(twoP))
};
// advance residues by 'lanes' steps (consistent with base advance)
for (i = 0; i < lanes; i++)
{
r10 += step10; if (r10 >= 20UL) r10 -= 20UL; if (r10 >= 10UL) r10 -= 10UL;
r8 = (r8 + step8) & 7UL;
r3 += step3; if (r3 >= 3UL) r3 -= 3UL;
r5 += step5; if (r5 >= 10UL) r5 -= 10UL; if (r5 >= 5UL) r5 -= 5UL;
}
divisorsExhausted = true;
}
}
}
return;
}
}
PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs
+50 -86
View File
@@ -12,36 +12,22 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
// via GpuKernelPool to avoid host round-trips.
// GPU residue variant: check 2^p % q == 1 for q = 2*p*k + 1.
public void Scan(
ulong exponent,
UInt128 twoP,
bool lastIsSeven,
UInt128 perSetLimit,
UInt128 setCount,
UInt128 overallLimit,
ref bool isPrime,
ref bool divisorsExhausted)
{
if (setCount == UInt128.Zero || perSetLimit == UInt128.Zero || overallLimit == UInt128.Zero)
{
return;
}
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var execution = gpuLease.EnterExecutionScope();
public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
{
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
var accelerator = gpuLease.Accelerator;
var stream = gpuLease.Stream;
// Ensure device has small cycles and primes tables for in-kernel lookup
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
int batchSize = GpuConstants.ScanBatchSize;
byte last = lastIsSeven ? (byte)1 : (byte)0;
var kernel = gpuLease.Pow2ModKernel;
ulong step10 = ((exponent % 10UL) << 1) % 10UL;
ulong step8 = ((exponent & 7UL) << 1) & 7UL;
ulong step3 = ((exponent % 3UL) << 1) % 3UL;
ulong step5 = ((exponent % 5UL) << 1) % 5UL;
GpuUInt128 twoPGpu = (GpuUInt128)twoP;
// Ensure device has small cycles and primes tables for in-kernel lookup
var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
int batchSize = GpuConstants.ScanBatchSize;
UInt128 kStart = 1UL;
byte last = lastIsSeven ? (byte)1 : (byte)0;
var kernel = gpuLease.Pow2ModKernel;
ulong step10 = ((exponent % 10UL) << 1) % 10UL;
ulong step8 = ((exponent & 7UL) << 1) & 7UL;
ulong step3 = ((exponent % 3UL) << 1) % 3UL;
ulong step5 = ((exponent % 5UL) << 1) % 5UL;
GpuUInt128 twoPGpu = (GpuUInt128)twoP;
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
// Device memory returned by Allocate1D is not guaranteed to be zeroed.
@@ -52,79 +38,57 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
// testers such as the order kernels.
orderBuffer.MemSetToZero();
ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
UInt128 batchSize128 = (UInt128)batchSize;
UInt128 limitInclusive = overallLimit == UInt128.MaxValue
? overallLimit
: checked(overallLimit + UInt128.One);
UInt128 remaining;
int currentSize;
int i;
UInt128 batchSize128 = (UInt128)batchSize, q;
Span<ulong> orders = orderArray.AsSpan(0, batchSize);
try
{
for (UInt128 setIndex = UInt128.Zero; setIndex < setCount && Volatile.Read(ref isPrime); setIndex++)
while (kStart < maxK && Volatile.Read(ref isPrime))
{
UInt128 setOffset = checked(perSetLimit * setIndex);
UInt128 setStart = checked(setOffset + UInt128.One);
if (setStart >= limitInclusive)
remaining = maxK - kStart;
if (remaining > batchSize128)
{
currentSize = batchSize;
}
else
{
currentSize = (int)remaining;
orders = orderArray.AsSpan(0, currentSize);
}
q = twoP * kStart + UInt128.One;
q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
kernel(currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL,
ra, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2);
accelerator.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
if (!Volatile.Read(ref isPrime))
{
break;
}
UInt128 setLimitExclusive = checked(setOffset + perSetLimit);
setLimitExclusive = setLimitExclusive >= limitInclusive
? limitInclusive
: checked(setLimitExclusive + UInt128.One);
if (setLimitExclusive > limitInclusive)
for (i = 0; i < currentSize; i++)
{
setLimitExclusive = limitInclusive;
}
UInt128 kStart = setStart;
while (kStart < setLimitExclusive && Volatile.Read(ref isPrime))
{
UInt128 remaining = checked(setLimitExclusive - kStart);
int currentSize = remaining > batchSize128 ? batchSize : (int)remaining;
if (currentSize <= 0)
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
break;
}
orders = orderArray.AsSpan(0, currentSize);
UInt128 q = checked(twoP * kStart);
q = checked(q + UInt128.One);
q.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
kernel(stream, currentSize, exponent, twoPGpu, (GpuUInt128)kStart, last, 0UL,
ra, orderBuffer.View, smallCyclesView, primeViews.LastOne, primeViews.LastSeven, primeViews.LastOnePow2, primeViews.LastSevenPow2);
stream.Synchronize();
orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
stream.Synchronize();
if (!Volatile.Read(ref isPrime))
{
break;
}
for (int i = 0; i < currentSize; i++)
{
if (orders[i] != 0UL)
{
Volatile.Write(ref isPrime, false);
break;
}
}
kStart = checked(kStart + (UInt128)currentSize);
}
kStart += batchSize128;
}
}
finally
{
ArrayPool<ulong>.Shared.Return(orderArray);
orderBuffer.Dispose();
execution.Dispose();
gpuLease.Dispose();
}
divisorsExhausted = true;
}
orderBuffer.Dispose();
gpuLease.Dispose();
}
}
}
PerfectNumbers.Core/MersenneNumberTester.cs
+62 -150
View File
@@ -16,28 +16,26 @@ public enum GpuKernelType
}
public sealed class MersenneNumberTester(
bool useIncremental = true,
bool useOrderCache = false,
GpuKernelType kernelType = GpuKernelType.Incremental,
bool useModuloWorkaround = false,
bool useOrder = false,
bool useGpuLucas = true,
bool useGpuScan = true,
bool useGpuOrder = true,
bool useResidue = true,
UInt128? maxK = null,
UInt128? residueDivisorSets = null)
bool useIncremental = true,
bool useOrderCache = false,
GpuKernelType kernelType = GpuKernelType.Incremental,
bool useModuloWorkaround = false,
bool useOrder = false,
bool useGpuLucas = true,
bool useGpuScan = true,
bool useGpuOrder = true,
bool useResidue = true,
ulong maxK = 5_000_000UL)
{
private readonly bool _useResidue = useResidue;
private readonly bool _useIncremental = useIncremental && !useResidue;
private readonly UInt128 _maxK = maxK ?? (UInt128)5_000_000UL;
private readonly UInt128 _residueDivisorSetCount = residueDivisorSets ?? (UInt128)PerfectNumberConstants.ExtraDivisorCycleSearchLimit;
private readonly bool _useResidue = useResidue;
private readonly bool _useIncremental = useIncremental && !useResidue;
private readonly ulong _maxK = maxK;
private readonly GpuKernelType _kernelType = kernelType;
private readonly bool _useModuloWorkaround = useModuloWorkaround;
private readonly bool _useGpuLucas = useGpuLucas;
private readonly bool _useGpuScan = useGpuScan; // device for pow2mod/incremental scanning
private readonly bool _useGpuOrder = useGpuOrder; // device for order computations
private static readonly ConcurrentDictionary<UInt128, ulong> OrderCache = new();
private readonly bool _useModuloWorkaround = useModuloWorkaround;
private readonly bool _useGpuLucas = useGpuLucas;
private readonly bool _useGpuScan = useGpuScan; // device for pow2mod/incremental scanning
private readonly bool _useGpuOrder = useGpuOrder; // device for order computations
private static readonly ConcurrentDictionary<UInt128, ulong> OrderCache = new();
private readonly MersenneNumberIncrementalGpuTester _incrementalGpuTester = new(kernelType, useGpuOrder);
private readonly MersenneNumberIncrementalCpuTester _incrementalCpuTester = new(kernelType);
private readonly MersenneNumberOrderGpuTester _orderGpuTester = new(kernelType, useGpuOrder);
@@ -50,24 +48,7 @@ public sealed class MersenneNumberTester(
// 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.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static UInt128 DivideRoundUp(UInt128 value, UInt128 divisor)
{
if (divisor == UInt128.Zero)
{
return UInt128.Zero;
}
UInt128 result = value / divisor;
if (result * divisor == value)
{
return result;
}
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;
if (!cacheEnabled)
@@ -218,116 +199,51 @@ public sealed class MersenneNumberTester(
}
}
public bool IsMersennePrime(ulong exponent) => IsMersennePrime(exponent, out _);
public bool IsMersennePrime(ulong exponent)
{
// 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.
if ((exponent % 3UL) == 0UL && exponent != 3UL)
{
return false;
}
public bool IsMersennePrime(ulong exponent, out bool divisorsExhausted)
{
divisorsExhausted = false;
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
{
return false;
}
// 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.
if ((exponent % 3UL) == 0UL && exponent != 3UL)
{
divisorsExhausted = true;
return false;
}
if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
{
divisorsExhausted = true;
return false;
}
// Skip even exponents: for even p > 2, M_p is composite; here p is large.
if ((exponent & 1UL) == 0UL)
{
divisorsExhausted = true;
return false;
}
// Skip even exponents: for even p > 2, M_p is composite; here p is large.
if ((exponent & 1UL) == 0UL)
{
return false;
}
bool prePrime = true;
UInt128 twoP = (UInt128)exponent << 1; // 2 * p
// UInt128 maxDivisor = new(ulong.MaxValue, ulong.MaxValue);
UInt128 maxK = UInt128Numbers.Two.Pow(exponent) / 2 + 1;
bool lastIsSeven = (exponent & 3UL) == 3UL;
UInt128 twoP = (UInt128)exponent << 1; // 2 * p
// UInt128 maxDivisor = new(ulong.MaxValue, ulong.MaxValue);
UInt128 maxK = UInt128Numbers.Two.Pow(exponent) / 2 + 1;
bool lastIsSeven = (exponent & 3UL) == 3UL;
// In residue mode the residue scan is the final primality check.
if (_useResidue)
{
if (_maxK < maxK)
{
maxK = _maxK;
}
// In residue mode the residue scan is the final primality check.
if (_useResidue)
{
if ((UInt128)_maxK < maxK)
{
maxK = _maxK;
}
if (_maxK == UInt128.Zero || _residueDivisorSetCount == UInt128.Zero)
{
return prePrime;
}
if (_useGpuScan)
{
_residueGpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime);
}
else
{
_residueCpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime);
}
UInt128 totalLimit = _maxK;
if (totalLimit > maxK)
{
totalLimit = maxK;
}
if (totalLimit == UInt128.Zero)
{
return true;
}
UInt128 requestedSets = _residueDivisorSetCount;
if (_useGpuScan)
{
UInt128 maxUsefulSets = DivideRoundUp(totalLimit, (UInt128)GpuConstants.ScanBatchSize);
if (maxUsefulSets == UInt128.Zero)
{
maxUsefulSets = UInt128.One;
}
if (requestedSets == UInt128.Zero || requestedSets > maxUsefulSets)
{
requestedSets = maxUsefulSets;
}
}
if (requestedSets == UInt128.Zero)
{
requestedSets = UInt128.One;
}
UInt128 perSetLimit = DivideRoundUp(totalLimit, requestedSets);
if (perSetLimit == UInt128.Zero)
{
perSetLimit = UInt128.One;
}
UInt128 setsNeeded = DivideRoundUp(totalLimit, perSetLimit);
UInt128 effectiveSets = requestedSets;
if (setsNeeded < effectiveSets)
{
effectiveSets = setsNeeded;
}
bool scanCompleted = false;
if (_useGpuScan)
{
_residueGpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted);
}
else
{
_residueCpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted);
}
if (!prePrime)
{
divisorsExhausted = true;
return false;
}
divisorsExhausted = scanCompleted;
return true;
}
return prePrime;
}
if (!_useIncremental)
{
@@ -347,15 +263,12 @@ public sealed class MersenneNumberTester(
if (!prePrime)
{
divisorsExhausted = true;
return false;
}
bool lucasPrime = _useGpuLucas
? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true)
: _lucasLehmerCpuTester.IsPrime(exponent);
divisorsExhausted = true;
return lucasPrime;
return _useGpuLucas
? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true)
: _lucasLehmerCpuTester.IsPrime(exponent);
}
bool isPrime = true;
@@ -383,9 +296,8 @@ public sealed class MersenneNumberTester(
}
}
divisorsExhausted = true;
return isPrime;
}
return isPrime;
}
// Limit of k to scan in the pre-Lucas prescan
public const int PreLucasPreScanK = 5_000_000;