Merge pull request #24 from mlemanczyk/codex/modify-evenperfectbitscanner-for-residue-mode
Report residue scan exhaustion through detailed checks
This commit is contained in:
7 files changed
+465
-224
No files matched your search
@@ -142,6 +142,95 @@ public class ProgramTests
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Divisor_mode_sets_detailedCheck_when_divisor_is_found()
|
||||
{
|
||||
var useDivisorField = typeof(Program).GetField("_useDivisor", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var divisorField = typeof(Program).GetField("_divisor", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var testerField = typeof(Program).GetField("_divisorTester", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var candidatesField = typeof(MersenneNumberDivisorGpuTester).GetField("_divisorCandidates", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var residueField = typeof(Program).GetField("PResidue", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu");
|
||||
|
||||
useDivisorField.SetValue(null, true);
|
||||
divisorField.SetValue(null, UInt128.Zero);
|
||||
testerField.SetValue(null, new MersenneNumberDivisorGpuTester());
|
||||
primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
|
||||
residueField.SetValue(null, new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true));
|
||||
candidatesField.SetValue(null, new (ulong, uint)[] { (23UL, 11U) });
|
||||
forceCpuProp!.SetValue(null, true);
|
||||
|
||||
try
|
||||
{
|
||||
Program.IsEvenPerfectCandidate(11UL, 32UL, out bool searched, out bool detailedCheck).Should().BeFalse();
|
||||
searched.Should().BeTrue();
|
||||
detailedCheck.Should().BeTrue();
|
||||
}
|
||||
finally
|
||||
{
|
||||
useDivisorField.SetValue(null, false);
|
||||
testerField.SetValue(null, null);
|
||||
candidatesField.SetValue(null, null);
|
||||
primeField.SetValue(null, null);
|
||||
residueField.SetValue(null, null);
|
||||
forceCpuProp!.SetValue(null, false);
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Residue_mode_sets_detailedCheck_based_on_exhaustion()
|
||||
{
|
||||
var useDivisorField = typeof(Program).GetField("_useDivisor", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var residueModeField = typeof(Program).GetField("_useResidueMode", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var mersenneField = typeof(Program).GetField("MersenneTesters", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var residueField = typeof(Program).GetField("PResidue", BindingFlags.NonPublic | BindingFlags.Static)!;
|
||||
var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu");
|
||||
|
||||
const ulong exponent = 7UL;
|
||||
|
||||
useDivisorField.SetValue(null, false);
|
||||
residueModeField.SetValue(null, true);
|
||||
primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
|
||||
residueField.SetValue(null, new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, exponent, true), trackAllValues: true));
|
||||
forceCpuProp!.SetValue(null, true);
|
||||
|
||||
try
|
||||
{
|
||||
mersenneField.SetValue(null, new ThreadLocal<MersenneNumberTester>(() => new MersenneNumberTester(
|
||||
useGpuScan: false,
|
||||
useGpuOrder: false,
|
||||
useResidue: true,
|
||||
maxK: 0UL,
|
||||
residueDivisorSets: 0UL), trackAllValues: true));
|
||||
|
||||
Program.IsEvenPerfectCandidate(exponent, 0UL, out bool searched, out bool detailedCheck).Should().BeTrue();
|
||||
searched.Should().BeTrue();
|
||||
detailedCheck.Should().BeFalse();
|
||||
|
||||
mersenneField.SetValue(null, new ThreadLocal<MersenneNumberTester>(() => new MersenneNumberTester(
|
||||
useGpuScan: false,
|
||||
useGpuOrder: false,
|
||||
useResidue: true,
|
||||
maxK: 4UL,
|
||||
residueDivisorSets: 2UL), trackAllValues: true));
|
||||
|
||||
Program.IsEvenPerfectCandidate(exponent, 0UL, out searched, out detailedCheck).Should().BeTrue();
|
||||
searched.Should().BeTrue();
|
||||
detailedCheck.Should().BeTrue();
|
||||
}
|
||||
finally
|
||||
{
|
||||
mersenneField.SetValue(null, null);
|
||||
primeField.SetValue(null, null);
|
||||
residueField.SetValue(null, null);
|
||||
residueModeField.SetValue(null, false);
|
||||
forceCpuProp!.SetValue(null, false);
|
||||
useDivisorField.SetValue(null, false);
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void Divisor_mode_accepts_large_search_limit()
|
||||
{
|
||||
@@ -223,7 +312,7 @@ public class ProgramTests
|
||||
var forceCpuProp = typeof(GpuContextPool).GetProperty("ForceCpu");
|
||||
|
||||
useDivisorField.SetValue(null, false);
|
||||
mersenneField.SetValue(null, new ThreadLocal<MersenneNumberTester>(() => new MersenneNumberTester(useIncremental: true, useResidue: true, maxK: 1_000UL), trackAllValues: true));
|
||||
mersenneField.SetValue(null, new ThreadLocal<MersenneNumberTester>(() => new MersenneNumberTester(useIncremental: true, useResidue: true, maxK: 1_000UL, residueDivisorSets: 1UL), trackAllValues: true));
|
||||
primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
|
||||
residueField.SetValue(null, new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true));
|
||||
forceCpuProp!.SetValue(null, false);
|
||||
@@ -232,7 +321,7 @@ public class ProgramTests
|
||||
{
|
||||
Program.IsEvenPerfectCandidate(29UL, 0UL, out bool searched, out bool detailed).Should().BeFalse();
|
||||
searched.Should().BeTrue();
|
||||
detailed.Should().BeFalse();
|
||||
detailed.Should().BeTrue();
|
||||
|
||||
Program.IsEvenPerfectCandidate(31UL, 0UL, out searched, out detailed).Should().BeTrue();
|
||||
searched.Should().BeTrue();
|
||||
|
||||
@@ -435,19 +435,20 @@ internal static class Program
|
||||
// Initialize per-thread p residue tracker (Identity model) at currentP
|
||||
if (!useDivisor)
|
||||
{
|
||||
MersenneTesters = new ThreadLocal<MersenneNumberTester>(() =>
|
||||
{
|
||||
var tester = new MersenneNumberTester(
|
||||
useIncremental: !useLucas,
|
||||
useOrderCache: false,
|
||||
kernelType: kernelType,
|
||||
useModuloWorkaround: useModuloWorkaround,
|
||||
useOrder: useOrder,
|
||||
useGpuLucas: mersenneOnGpu,
|
||||
useGpuScan: mersenneOnGpu,
|
||||
useGpuOrder: orderOnGpu,
|
||||
useResidue: useResidue,
|
||||
maxK: residueKMax);
|
||||
MersenneTesters = new ThreadLocal<MersenneNumberTester>(() =>
|
||||
{
|
||||
var tester = new MersenneNumberTester(
|
||||
useIncremental: !useLucas,
|
||||
useOrderCache: false,
|
||||
kernelType: kernelType,
|
||||
useModuloWorkaround: useModuloWorkaround,
|
||||
useOrder: useOrder,
|
||||
useGpuLucas: mersenneOnGpu,
|
||||
useGpuScan: mersenneOnGpu,
|
||||
useGpuOrder: orderOnGpu,
|
||||
useResidue: useResidue,
|
||||
maxK: residueKMax,
|
||||
residueDivisorSets: divisorCyclesSearchLimit);
|
||||
if (!useLucas)
|
||||
{
|
||||
tester.WarmUpOrders(currentP, _orderWarmupLimitOverride ?? 5_000_000UL);
|
||||
@@ -731,7 +732,7 @@ internal static class Program
|
||||
Console.WriteLine(" --divisor-cycles-device=cpu|gpu device for cycles generation (default gpu)");
|
||||
Console.WriteLine(" --divisor-cycles-batch-size=<value> batch size for cycles generation (default 512)");
|
||||
Console.WriteLine(" --divisor-cycles-continue continue divisor cycles generation");
|
||||
Console.WriteLine(" --divisor-cycles-limit=<value> cycle search iterations when --mersenne=divisor");
|
||||
Console.WriteLine(" --divisor-cycles-limit=<value> cycle search iterations when --mersenne=divisor or residue");
|
||||
|
||||
Console.WriteLine(" --use-order test primality via q order");
|
||||
Console.WriteLine(" --workaround-mod avoid '%' operator on the GPU");
|
||||
@@ -1077,15 +1078,16 @@ internal static class Program
|
||||
return false;
|
||||
}
|
||||
|
||||
searchedMersenne = true;
|
||||
if (_useDivisor)
|
||||
{
|
||||
return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck);
|
||||
}
|
||||
searchedMersenne = true;
|
||||
if (_useDivisor)
|
||||
{
|
||||
return _divisorTester!.IsPrime(p, _divisor, divisorCyclesSearchLimit, out detailedCheck);
|
||||
}
|
||||
|
||||
detailedCheck = MersenneTesters.Value!.IsMersennePrime(p);
|
||||
return detailedCheck;
|
||||
}
|
||||
bool mersennePrime = MersenneTesters.Value!.IsMersennePrime(p, out bool residueExhausted);
|
||||
detailedCheck = residueExhausted;
|
||||
return mersennePrime;
|
||||
}
|
||||
|
||||
// Use ModResidueTracker with a small set of primes to pre-filter composite p.
|
||||
private static bool IsCompositeByResidues(ulong p)
|
||||
|
||||
@@ -24,8 +24,18 @@ public class MersenneNumberResidueCpuTesterTests
|
||||
private static void RunCase(MersenneNumberResidueCpuTester tester, ulong exponent, ulong maxK, bool expectedPrime)
|
||||
{
|
||||
bool isPrime = true;
|
||||
tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime);
|
||||
bool exhausted = false;
|
||||
tester.Scan(
|
||||
exponent,
|
||||
(UInt128)exponent << 1,
|
||||
LastDigitIsSeven(exponent),
|
||||
(UInt128)maxK,
|
||||
1UL,
|
||||
(UInt128)maxK,
|
||||
ref isPrime,
|
||||
ref exhausted);
|
||||
isPrime.Should().Be(expectedPrime);
|
||||
exhausted.Should().BeTrue();
|
||||
}
|
||||
}
|
||||
|
||||
@@ -27,8 +27,18 @@ public class MersenneNumberResidueGpuTesterTests
|
||||
private static void RunCase(MersenneNumberResidueGpuTester tester, ulong exponent, ulong maxK, bool expectedPrime)
|
||||
{
|
||||
bool isPrime = true;
|
||||
tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime);
|
||||
bool exhausted = false;
|
||||
tester.Scan(
|
||||
exponent,
|
||||
(UInt128)exponent << 1,
|
||||
LastDigitIsSeven(exponent),
|
||||
(UInt128)maxK,
|
||||
1UL,
|
||||
(UInt128)maxK,
|
||||
ref isPrime,
|
||||
ref exhausted);
|
||||
isPrime.Should().Be(expectedPrime);
|
||||
exhausted.Should().BeTrue();
|
||||
}
|
||||
}
|
||||
|
||||
@@ -5,25 +5,32 @@ 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 maxK, ref bool isPrime)
|
||||
{
|
||||
// Initialize/update Mersenne residue tracker and start a merge walk for ascending q divisors
|
||||
public void Scan(
|
||||
ulong exponent,
|
||||
UInt128 twoP,
|
||||
bool lastIsSeven,
|
||||
UInt128 perSetLimit,
|
||||
ulong setCount,
|
||||
UInt128 overallLimit,
|
||||
ref bool isPrime,
|
||||
ref bool divisorsExhausted)
|
||||
{
|
||||
if (setCount == 0UL || perSetLimit == UInt128.Zero || overallLimit == UInt128.Zero)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
_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);
|
||||
// Initialize/update Mersenne residue tracker and start a merge walk for ascending q divisors
|
||||
|
||||
// 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);
|
||||
_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);
|
||||
|
||||
// Allowed last-digit sets for q depending on last digit of M_p
|
||||
// lastIsSeven == true => allow {7,9}
|
||||
@@ -31,139 +38,162 @@ 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;
|
||||
|
||||
maxK++;
|
||||
ModResidueTracker tracker = _mersenneResidueTracker!;
|
||||
UInt128 limitInclusive = overallLimit + UInt128.One;
|
||||
for (ulong setIndex = 0; setIndex < setCount && Volatile.Read(ref isPrime); setIndex++)
|
||||
{
|
||||
UInt128 setOffset = perSetLimit * (UInt128)setIndex;
|
||||
UInt128 k = setOffset + UInt128.One;
|
||||
if (k >= limitInclusive)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
ModResidueTracker tracker = _mersenneResidueTracker!;
|
||||
while (k < maxK && (localIsPrime = Volatile.Read(ref isPrime)))
|
||||
{
|
||||
remaining = maxK - k;
|
||||
if (remaining == UInt128.Zero)
|
||||
{
|
||||
break;
|
||||
}
|
||||
UInt128 setLimitExclusive = setOffset + perSetLimit + UInt128.One;
|
||||
if (setLimitExclusive > limitInclusive)
|
||||
{
|
||||
setLimitExclusive = limitInclusive;
|
||||
}
|
||||
|
||||
lanes = remaining > 4UL ? 4 : (int)remaining;
|
||||
UInt128 q = twoP * k + UInt128.One;
|
||||
q.Mod10_8_5_3(out ulong r10, out ulong r8, out ulong r5, out ulong r3);
|
||||
|
||||
// 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;
|
||||
while (k < setLimitExclusive && (localIsPrime = Volatile.Read(ref isPrime)))
|
||||
{
|
||||
remaining = setLimitExclusive - k;
|
||||
if (remaining == UInt128.Zero)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
// r8
|
||||
r8_0 = r8 & 7UL;
|
||||
r8_1 = (r8_0 + step8) & 7UL;
|
||||
r8_2 = (r8_1 + step8) & 7UL;
|
||||
r8_3 = (r8_2 + step8) & 7UL;
|
||||
lanes = remaining > 4UL ? 4 : (int)remaining;
|
||||
if (lanes <= 0)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
// 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;
|
||||
// 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;
|
||||
|
||||
// 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;
|
||||
// r8
|
||||
r8_0 = r8 & 7UL;
|
||||
r8_1 = (r8_0 + step8) & 7UL;
|
||||
r8_2 = (r8_1 + step8) & 7UL;
|
||||
r8_3 = (r8_2 + step8) & 7UL;
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
// 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 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;
|
||||
}
|
||||
}
|
||||
// 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 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 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 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 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;
|
||||
}
|
||||
}
|
||||
|
||||
if (!localIsPrime)
|
||||
{
|
||||
Volatile.Write(ref isPrime, false);
|
||||
return;
|
||||
}
|
||||
else if (!Volatile.Read(ref isPrime))
|
||||
{
|
||||
return;
|
||||
}
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
// 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))
|
||||
};
|
||||
// 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 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;
|
||||
}
|
||||
}
|
||||
if (!localIsPrime)
|
||||
{
|
||||
Volatile.Write(ref isPrime, false);
|
||||
divisorsExhausted = true;
|
||||
return;
|
||||
}
|
||||
else if (!Volatile.Read(ref isPrime))
|
||||
{
|
||||
divisorsExhausted = true;
|
||||
return;
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -11,16 +11,28 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
|
||||
// candidates without launching heavy kernels. Consider a device-visible constant buffer per-accelerator
|
||||
// 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 maxK, ref bool isPrime)
|
||||
{
|
||||
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
|
||||
// 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,
|
||||
ulong setCount,
|
||||
UInt128 overallLimit,
|
||||
ref bool isPrime,
|
||||
ref bool divisorsExhausted)
|
||||
{
|
||||
if (setCount == 0UL || perSetLimit == UInt128.Zero || overallLimit == UInt128.Zero)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
|
||||
var accelerator = gpuLease.Accelerator;
|
||||
// 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;
|
||||
@@ -31,57 +43,71 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
|
||||
|
||||
var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
|
||||
ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
|
||||
UInt128 remaining;
|
||||
int currentSize;
|
||||
int i;
|
||||
UInt128 batchSize128 = (UInt128)batchSize, q;
|
||||
UInt128 batchSize128 = (UInt128)batchSize;
|
||||
UInt128 limitInclusive = overallLimit + UInt128.One;
|
||||
Span<ulong> orders = orderArray.AsSpan(0, batchSize);
|
||||
try
|
||||
{
|
||||
while (kStart < maxK && Volatile.Read(ref isPrime))
|
||||
for (ulong setIndex = 0; setIndex < setCount && Volatile.Read(ref isPrime); setIndex++)
|
||||
{
|
||||
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))
|
||||
UInt128 setOffset = perSetLimit * (UInt128)setIndex;
|
||||
UInt128 setStart = setOffset + UInt128.One;
|
||||
if (setStart >= limitInclusive)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
for (i = 0; i < currentSize; i++)
|
||||
UInt128 setLimitExclusive = setOffset + perSetLimit + UInt128.One;
|
||||
if (setLimitExclusive > limitInclusive)
|
||||
{
|
||||
if (orders[i] != 0UL)
|
||||
{
|
||||
Volatile.Write(ref isPrime, false);
|
||||
break;
|
||||
}
|
||||
setLimitExclusive = limitInclusive;
|
||||
}
|
||||
|
||||
kStart += batchSize128;
|
||||
UInt128 kStart = setStart;
|
||||
while (kStart < setLimitExclusive && Volatile.Read(ref isPrime))
|
||||
{
|
||||
UInt128 remaining = setLimitExclusive - kStart;
|
||||
int currentSize = remaining > batchSize128 ? batchSize : (int)remaining;
|
||||
if (currentSize <= 0)
|
||||
{
|
||||
break;
|
||||
}
|
||||
orders = orderArray.AsSpan(0, currentSize);
|
||||
|
||||
UInt128 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;
|
||||
}
|
||||
|
||||
for (int i = 0; i < currentSize; i++)
|
||||
{
|
||||
if (orders[i] != 0UL)
|
||||
{
|
||||
Volatile.Write(ref isPrime, false);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
kStart += batchSize128;
|
||||
}
|
||||
}
|
||||
}
|
||||
finally
|
||||
{
|
||||
ArrayPool<ulong>.Shared.Return(orderArray);
|
||||
orderBuffer.Dispose();
|
||||
gpuLease.Dispose();
|
||||
}
|
||||
}
|
||||
orderBuffer.Dispose();
|
||||
gpuLease.Dispose();
|
||||
}
|
||||
|
||||
divisorsExhausted = true;
|
||||
}
|
||||
}
|
||||
@@ -24,13 +24,15 @@ public sealed class MersenneNumberTester(
|
||||
bool useGpuLucas = true,
|
||||
bool useGpuScan = true,
|
||||
bool useGpuOrder = true,
|
||||
bool useResidue = true,
|
||||
ulong maxK = 5_000_000UL)
|
||||
bool useResidue = true,
|
||||
ulong maxK = 5_000_000UL,
|
||||
ulong residueDivisorSets = 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 ulong _residueDivisorSetCount = residueDivisorSets;
|
||||
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
|
||||
@@ -42,11 +44,28 @@ public sealed class MersenneNumberTester(
|
||||
private readonly MersenneNumberOrderCpuTester _orderCpuTester = new(kernelType);
|
||||
private readonly MersenneNumberLucasLehmerGpuTester _lucasLehmerGpuTester = new();
|
||||
private readonly MersenneNumberLucasLehmerCpuTester _lucasLehmerCpuTester = new();
|
||||
private readonly MersenneNumberResidueGpuTester? _residueGpuTester = useResidue ? new(useGpuOrder) : null;
|
||||
private readonly MersenneNumberResidueCpuTester? _residueCpuTester = useResidue ? new() : null;
|
||||
private readonly MersenneNumberResidueGpuTester? _residueGpuTester = useResidue ? new(useGpuOrder) : 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.
|
||||
// Also consider injecting a shared cycles cache (MersenneDivisorCycles.Shared) to both CPU and GPU testers.
|
||||
// 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)
|
||||
{
|
||||
@@ -191,27 +210,34 @@ public sealed class MersenneNumberTester(
|
||||
ArrayPool<UInt128>.Shared.Return(qs);
|
||||
}
|
||||
|
||||
public bool IsMersennePrime(ulong exponent)
|
||||
public bool IsMersennePrime(ulong exponent) => IsMersennePrime(exponent, out _);
|
||||
|
||||
public bool IsMersennePrime(ulong exponent, out bool divisorsExhausted)
|
||||
{
|
||||
divisorsExhausted = 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;
|
||||
}
|
||||
|
||||
bool prePrime = true;
|
||||
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;
|
||||
@@ -225,16 +251,60 @@ public sealed class MersenneNumberTester(
|
||||
maxK = _maxK;
|
||||
}
|
||||
|
||||
if (_maxK == 0UL || _residueDivisorSetCount == 0UL)
|
||||
{
|
||||
return prePrime;
|
||||
}
|
||||
|
||||
UInt128 totalLimit = (UInt128)_maxK;
|
||||
if (totalLimit > maxK)
|
||||
{
|
||||
totalLimit = maxK;
|
||||
}
|
||||
|
||||
if (totalLimit == UInt128.Zero)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
ulong requestedSets = _residueDivisorSetCount;
|
||||
if (requestedSets == 0UL)
|
||||
{
|
||||
requestedSets = 1UL;
|
||||
}
|
||||
|
||||
UInt128 perSetLimit = DivideRoundUp(totalLimit, (UInt128)requestedSets);
|
||||
if (perSetLimit == UInt128.Zero)
|
||||
{
|
||||
perSetLimit = UInt128.One;
|
||||
}
|
||||
|
||||
UInt128 setsNeeded = DivideRoundUp(totalLimit, perSetLimit);
|
||||
ulong effectiveSets = requestedSets;
|
||||
if (setsNeeded < (UInt128)effectiveSets)
|
||||
{
|
||||
effectiveSets = (ulong)setsNeeded;
|
||||
}
|
||||
|
||||
bool scanCompleted = false;
|
||||
|
||||
if (_useGpuScan)
|
||||
{
|
||||
_residueGpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime);
|
||||
_residueGpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted);
|
||||
}
|
||||
else
|
||||
{
|
||||
_residueCpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime);
|
||||
_residueCpuTester!.Scan(exponent, twoP, lastIsSeven, perSetLimit, effectiveSets, totalLimit, ref prePrime, ref scanCompleted);
|
||||
}
|
||||
|
||||
return prePrime;
|
||||
if (!prePrime)
|
||||
{
|
||||
divisorsExhausted = true;
|
||||
return false;
|
||||
}
|
||||
|
||||
divisorsExhausted = scanCompleted;
|
||||
return true;
|
||||
}
|
||||
|
||||
if (!_useIncremental)
|
||||
@@ -253,15 +323,18 @@ public sealed class MersenneNumberTester(
|
||||
_incrementalCpuTester.Scan(exponent, twoPPre, lastIsSevenPre, maxKPre, ref prePrime);
|
||||
}
|
||||
|
||||
if (!prePrime)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
if (!prePrime)
|
||||
{
|
||||
divisorsExhausted = true;
|
||||
return false;
|
||||
}
|
||||
|
||||
return _useGpuLucas
|
||||
? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true)
|
||||
: _lucasLehmerCpuTester.IsPrime(exponent);
|
||||
}
|
||||
bool lucasPrime = _useGpuLucas
|
||||
? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true)
|
||||
: _lucasLehmerCpuTester.IsPrime(exponent);
|
||||
divisorsExhausted = true;
|
||||
return lucasPrime;
|
||||
}
|
||||
|
||||
bool isPrime = true;
|
||||
|
||||
@@ -288,6 +361,7 @@ public sealed class MersenneNumberTester(
|
||||
}
|
||||
}
|
||||
|
||||
divisorsExhausted = true;
|
||||
return isPrime;
|
||||
}
|
||||
|
||||
|
||||
Reference in new issue
Block a user