From 1b285841682cefb1bf536436c27192ec2f3a56eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 17 Sep 2025 19:02:21 +0200 Subject: [PATCH] Add divisor-mode detailed check coverage for factors --- EvenPerfectBitScanner.Tests/ProgramTests.cs | 93 +++++- EvenPerfectBitScanner/Program.cs | 46 +-- .../MersenneNumberResidueCpuTesterTests.cs | 12 +- .../MersenneNumberResidueGpuTesterTests.cs | 12 +- .../Cpu/MersenneNumberResidueCpuTester.cs | 302 ++++++++++-------- .../Gpu/MersenneNumberResidueGpuTester.cs | 110 ++++--- PerfectNumbers.Core/MersenneNumberTester.cs | 114 +++++-- 7 files changed, 465 insertions(+), 224 deletions(-) diff --git a/EvenPerfectBitScanner.Tests/ProgramTests.cs b/EvenPerfectBitScanner.Tests/ProgramTests.cs index 1d72484..c61b3f0 100644 --- a/EvenPerfectBitScanner.Tests/ProgramTests.cs +++ b/EvenPerfectBitScanner.Tests/ProgramTests.cs @@ -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(() => new PrimeTester(), trackAllValues: true)); + residueField.SetValue(null, new ThreadLocal(() => 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(() => new PrimeTester(), trackAllValues: true)); + residueField.SetValue(null, new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, exponent, true), trackAllValues: true)); + forceCpuProp!.SetValue(null, true); + + try + { + mersenneField.SetValue(null, new ThreadLocal(() => 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(() => 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(() => new MersenneNumberTester(useIncremental: true, useResidue: true, maxK: 1_000UL), trackAllValues: true)); + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester(useIncremental: true, useResidue: true, maxK: 1_000UL, residueDivisorSets: 1UL), trackAllValues: true)); primeField.SetValue(null, new ThreadLocal(() => new PrimeTester(), trackAllValues: true)); residueField.SetValue(null, new ThreadLocal(() => 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(); diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 7502d7e..f9ddf0d 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -435,19 +435,20 @@ internal static class Program // Initialize per-thread p residue tracker (Identity model) at currentP if (!useDivisor) { - MersenneTesters = new ThreadLocal(() => - { - 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(() => + { + 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= batch size for cycles generation (default 512)"); Console.WriteLine(" --divisor-cycles-continue continue divisor cycles generation"); - Console.WriteLine(" --divisor-cycles-limit= cycle search iterations when --mersenne=divisor"); + Console.WriteLine(" --divisor-cycles-limit= 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) diff --git a/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs index e083d11..8885e0c 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs @@ -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(); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs index 3aa2447..5e44cb6 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs @@ -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(); } } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs index 9d912c5..6e37b3b 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs @@ -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; + } +} + } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index b8cfbd2..598de60 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -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(batchSize); ulong[] orderArray = ArrayPool.Shared.Rent(batchSize); - UInt128 remaining; - int currentSize; - int i; - UInt128 batchSize128 = (UInt128)batchSize, q; + UInt128 batchSize128 = (UInt128)batchSize; + UInt128 limitInclusive = overallLimit + UInt128.One; Span 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.Shared.Return(orderArray); - orderBuffer.Dispose(); - gpuLease.Dispose(); - } - } + orderBuffer.Dispose(); + gpuLease.Dispose(); + } + + divisorsExhausted = true; + } } diff --git a/PerfectNumbers.Core/MersenneNumberTester.cs b/PerfectNumbers.Core/MersenneNumberTester.cs index 31a5d2e..fefc38f 100644 --- a/PerfectNumbers.Core/MersenneNumberTester.cs +++ b/PerfectNumbers.Core/MersenneNumberTester.cs @@ -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.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; }