From a21ab3134e0bf1fe6a47a563dd9cde3d5fe3b850 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Fri, 19 Dec 2025 18:25:27 +0100 Subject: [PATCH] MersenneNumberTester converted to device dependent MersenneNumberTesterTemplate. --- EvenPerfectBitScanner.Tests/ProgramTests.cs | 15 +- EvenPerfectBitScanner/Program.cs | 81 ++- .../MersenneNumberTesterTests.cs | 12 +- .../MersenneResidueModeTests.cs | 6 +- PerfectNumbers.Core/MersenneNumberTester.cs | 642 ------------------ .../MersenneNumberTesterTemplate.cs | 352 ++++++++++ 6 files changed, 432 insertions(+), 676 deletions(-) delete mode 100644 PerfectNumbers.Core/MersenneNumberTester.cs create mode 100644 PerfectNumbers.Core/MersenneNumberTesterTemplate.cs diff --git a/EvenPerfectBitScanner.Tests/ProgramTests.cs b/EvenPerfectBitScanner.Tests/ProgramTests.cs index b0b2437..8745a19 100644 --- a/EvenPerfectBitScanner.Tests/ProgramTests.cs +++ b/EvenPerfectBitScanner.Tests/ProgramTests.cs @@ -227,7 +227,7 @@ public class ProgramTests var mersenneField = typeof(Program).GetField("MersenneTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; - mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester(useIncremental: true, useResidue: true, maxK: 1_000UL), trackAllValues: true)); + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTesterCpu(useIncremental: true, useResidue: true, maxK: 1_000UL), trackAllValues: true)); primeField.SetValue(null, new ThreadLocal(() => new PrimeTester(), trackAllValues: true)); CandidatesCalculator.OverrideResidueTrackers(new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true)); @@ -260,7 +260,7 @@ public class ProgramTests var mersenneField = typeof(Program).GetField("MersenneTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; - mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester( + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTesterCpu( useIncremental: true, useResidue: true, maxK: 1_024UL), trackAllValues: true)); @@ -297,7 +297,7 @@ public class ProgramTests var mersenneField = typeof(Program).GetField("MersenneTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; - mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester( + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTesterCpu( useIncremental: true, useResidue: true, maxK: 5_000_000UL), trackAllValues: true)); @@ -333,7 +333,7 @@ public class ProgramTests var mersenneField = typeof(Program).GetField("MersenneTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; - mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester( + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTesterCpu( useIncremental: true, useResidue: true, maxK: 5_000_000UL), trackAllValues: true)); @@ -370,7 +370,7 @@ public class ProgramTests var mersenneField = typeof(Program).GetField("MersenneTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; - mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester( + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTesterCpu( useIncremental: true, useGpuScan: false, useGpuOrder: false, @@ -406,7 +406,7 @@ public class ProgramTests var mersenneField = typeof(Program).GetField("MersenneTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; - mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester( + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTesterCpu( useIncremental: false, kernelType: configuration.KernelType, useGpuScan: true, @@ -493,7 +493,7 @@ public class ProgramTests var primeField = typeof(Program).GetField("PrimeTesters", BindingFlags.NonPublic | BindingFlags.Static)!; var compositeField = typeof(Program).GetField("_lastCompositeP", BindingFlags.NonPublic | BindingFlags.Static)!; - mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTester(useIncremental: true, useResidue: false, maxK: 1_000UL), trackAllValues: true)); + mersenneField.SetValue(null, new ThreadLocal(() => new MersenneNumberTesterCpu(useIncremental: true, useResidue: false, maxK: 1_000UL), trackAllValues: true)); primeField.SetValue(null, new ThreadLocal(() => new PrimeTester(), trackAllValues: true)); CandidatesCalculator.OverrideResidueTrackers(new ThreadLocal(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true)); @@ -547,3 +547,4 @@ public class ProgramTests } } + diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index d00d83a..de91414 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -14,7 +14,7 @@ namespace EvenPerfectBitScanner; internal static class Program { // private static ThreadLocal PrimeTesters = null!; - private static ThreadLocal MersenneTesters = null!; +private static ThreadLocal? MersenneTesters; private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP; private static MersenneNumberDivisorGpuTester? _divisorTester; private static UInt128 _divisor; @@ -27,6 +27,8 @@ internal static class Program private static CliArguments _cliArguments; private static int _testTargetPrimeCount; private static int _testProcessedPrimeCount; + private static ComputationDevice _mersenneDevice; + private static ComputationDevice _orderDevice; [ThreadStatic] private static bool _lastCompositeP; @@ -97,6 +99,8 @@ internal static class Program CandidatesCalculator.Configure(transformMode); ComputationDevice mersenneDevice = _cliArguments.MersenneDevice; ComputationDevice orderDevice = EnvironmentConfiguration.OrderDevice; + _mersenneDevice = mersenneDevice; + _orderDevice = orderDevice; bool mersenneOnGpu = mersenneDevice == ComputationDevice.Gpu; bool orderOnGpu = EnvironmentConfiguration.UseGpuOrder; bool orderOnHybrid = EnvironmentConfiguration.UseHybridOrder; @@ -214,34 +218,59 @@ internal static class Program // Initialize per-thread p residue tracker (Identity model) at currentP if (!useDivisor && !_cliArguments.UseByDivisor) { - MersenneTesters = new ThreadLocal(() => + MersenneTesters = new ThreadLocal(() => { // ProcessEightBitWindows windowed pow2 ladder is the default kernel. - var tester = new MersenneNumberTester( - useIncremental: !useLucas, - useOrderCache: false, - kernelType: kernelType, - useOrder: useOrder, - useGpuLucas: mersenneOnGpu, - useGpuScan: mersenneOnGpu, - useGpuOrder: orderOnGpu || orderOnHybrid, - useResidue: useResidue, - maxK: residueKMax, - orderDevice: orderDevice); + object tester = orderDevice switch + { + ComputationDevice.Gpu => new MersenneNumberTesterGpu( + useIncremental: !useLucas, + useOrderCache: false, + kernelType: kernelType, + useOrder: useOrder, + useGpuLucas: mersenneOnGpu, + useGpuScan: mersenneOnGpu, + useGpuOrder: orderOnGpu || orderOnHybrid, + useResidue: useResidue, + maxK: residueKMax, + orderDevice: orderDevice), + ComputationDevice.Hybrid => new MersenneNumberTesterHybrid( + useIncremental: !useLucas, + useOrderCache: false, + kernelType: kernelType, + useOrder: useOrder, + useGpuLucas: mersenneOnGpu, + useGpuScan: mersenneOnGpu, + useGpuOrder: orderOnGpu || orderOnHybrid, + useResidue: useResidue, + maxK: residueKMax, + orderDevice: orderDevice), + _ => new MersenneNumberTesterCpu( + useIncremental: !useLucas, + useOrderCache: false, + kernelType: kernelType, + useOrder: useOrder, + useGpuLucas: mersenneOnGpu, + useGpuScan: mersenneOnGpu, + useGpuOrder: orderOnGpu || orderOnHybrid, + useResidue: useResidue, + maxK: residueKMax, + orderDevice: orderDevice), + }; if (!useLucas) { Console.WriteLine("Warming up orders"); if (orderDevice == ComputationDevice.Hybrid) { - tester.WarmUpOrdersHybrid(currentP, orderWarmupLimitOverride ?? 5_000_000UL); + ((MersenneNumberTesterHybrid)tester).WarmUpOrders(currentP, orderWarmupLimitOverride ?? 5_000_000UL); } else if (orderDevice == ComputationDevice.Gpu) { - tester.WarmUpOrdersGpu(currentP, orderWarmupLimitOverride ?? 5_000_000UL); + ((MersenneNumberTesterGpu)tester).WarmUpOrders(currentP, orderWarmupLimitOverride ?? 5_000_000UL); } else { - tester.WarmUpOrdersCpu(currentP, orderWarmupLimitOverride ?? 5_000_000UL); + ((MersenneNumberTesterCpu)tester).WarmUpOrders(currentP, orderWarmupLimitOverride ?? 5_000_000UL); } } @@ -548,15 +577,15 @@ internal static class Program Console.WriteLine("Warming up orders..."); if (orderDevice == ComputationDevice.Hybrid) { - MersenneTesters.Value!.WarmUpOrdersHybrid(currentP, orderWarmupLimitOverride ?? 5_000_000UL); + ((MersenneNumberTesterHybrid)MersenneTesters.Value!).WarmUpOrders(currentP, orderWarmupLimitOverride ?? 5_000_000UL); } else if (orderDevice == ComputationDevice.Gpu) { - MersenneTesters.Value!.WarmUpOrdersGpu(currentP, orderWarmupLimitOverride ?? 5_000_000UL); + ((MersenneNumberTesterGpu)MersenneTesters.Value!).WarmUpOrders(currentP, orderWarmupLimitOverride ?? 5_000_000UL); } else { - MersenneTesters.Value!.WarmUpOrdersCpu(currentP, orderWarmupLimitOverride ?? 5_000_000UL); + ((MersenneNumberTesterCpu)MersenneTesters.Value!).WarmUpOrders(currentP, orderWarmupLimitOverride ?? 5_000_000UL); } } } @@ -664,6 +693,17 @@ internal static class Program } } + private static bool IsMersennePrime(PrimeOrderCalculatorAccelerator gpu, ulong exponent) + { + object tester = MersenneTesters!.Value!; + return _orderDevice switch + { + ComputationDevice.Gpu => ((MersenneNumberTesterGpu)tester).IsMersennePrime(gpu, exponent), + ComputationDevice.Hybrid => ((MersenneNumberTesterHybrid)tester).IsMersennePrime(gpu, exponent), + _ => ((MersenneNumberTesterCpu)tester).IsMersennePrime(gpu, exponent), + }; + } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void OnPrimeCandidateConfirmed() { @@ -781,7 +821,8 @@ internal static class Program return _byDivisorGpuOrderTester.IsPrime(p, out detailedCheck, out _); } - detailedCheck = MersenneTesters.Value!.IsMersennePrime(gpu, p); + detailedCheck = IsMersennePrime(gpu, p); return detailedCheck; } } + diff --git a/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneNumberTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneNumberTesterTests.cs index c81993b..68cbb75 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneNumberTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneNumberTesterTests.cs @@ -13,7 +13,7 @@ public class MersenneNumberTesterTests [InlineData(136_000_005UL)] public void IsMersennePrime_detects_small_prime_divisors_for_large_exponents(ulong p) { - var tester = new MersenneNumberTester(); + var tester = new MersenneNumberTesterCpu(); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { @@ -29,8 +29,8 @@ public class MersenneNumberTesterTests [InlineData(125UL)] public void WarmUpOrdersCpu_populates_cache_without_affecting_results(ulong p) { - var tester = new MersenneNumberTester(useOrderCache: true, useGpuOrder: true); - tester.WarmUpOrdersCpu(p, 1_000UL); + var tester = new MersenneNumberTesterCpu(useOrderCache: true, useGpuOrder: true); + tester.WarmUpOrders(p, 1_000UL); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { @@ -46,8 +46,8 @@ public class MersenneNumberTesterTests [InlineData(125UL)] public void WarmUpOrdersGpu_populates_cache_without_affecting_results(ulong p) { - var tester = new MersenneNumberTester(useOrderCache: true, useGpuOrder: true); - tester.WarmUpOrdersGpu(p, 1_000UL); + var tester = new MersenneNumberTesterGpu(useOrderCache: true, useGpuOrder: true); + tester.WarmUpOrders(p, 1_000UL); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); try { @@ -84,3 +84,5 @@ public class MersenneNumberTesterTests } } + + diff --git a/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneResidueModeTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneResidueModeTests.cs index e428761..2ff289a 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneResidueModeTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberTester/MersenneResidueModeTests.cs @@ -10,7 +10,7 @@ public class MersenneResidueModeTests [Fact] public void IsMersennePrime_residue_mode_rejects_composite_exponent() { - var tester = new MersenneNumberTester( + var tester = new MersenneNumberTesterCpu( useResidue: true, useIncremental: false, useGpuScan: false, @@ -29,7 +29,7 @@ public class MersenneResidueModeTests [Fact] public void IsMersennePrime_residue_mode_accepts_prime_exponent() { - var tester = new MersenneNumberTester( + var tester = new MersenneNumberTesterCpu( useResidue: true, useIncremental: true, useGpuScan: false, @@ -47,3 +47,5 @@ public class MersenneResidueModeTests } } + + diff --git a/PerfectNumbers.Core/MersenneNumberTester.cs b/PerfectNumbers.Core/MersenneNumberTester.cs deleted file mode 100644 index d76ee01..0000000 --- a/PerfectNumbers.Core/MersenneNumberTester.cs +++ /dev/null @@ -1,642 +0,0 @@ -using System.Runtime.CompilerServices; -using System.Collections.Concurrent; -using System.Buffers; -using ILGPU; -using ILGPU.Runtime; -using PerfectNumbers.Core.Cpu; -using PerfectNumbers.Core.Gpu; -using UInt128 = System.UInt128; -using PerfectNumbers.Core.Gpu.Accelerators; - -namespace PerfectNumbers.Core; - -public enum GpuKernelType -{ - Incremental, - Pow2Mod, -} - -public sealed class MersenneNumberTester( - bool useIncremental = true, - bool useOrderCache = false, - GpuKernelType kernelType = GpuKernelType.Incremental, - bool useOrder = false, - bool useGpuLucas = true, - bool useGpuScan = true, - bool useGpuOrder = true, - bool useResidue = true, - ulong maxK = 5_000_000UL, - ComputationDevice orderDevice = ComputationDevice.Gpu) -{ - private readonly bool _useResidue = useResidue; - private readonly bool _useIncremental = useIncremental && !useResidue; - private readonly ulong _maxK = maxK; - private readonly Func _sharesFactorWithExponentMinusOne = orderDevice switch - { - ComputationDevice.Hybrid => static (gpu, exponent) => exponent.SharesFactorWithExponentMinusOneHybrid(gpu), - _ when useGpuScan => static (gpu, exponent) => exponent.SharesFactorWithExponentMinusOneGpu(gpu), - _ => static (_, exponent) => exponent.SharesFactorWithExponentMinusOneCpu(), - }; - - private static readonly ConcurrentDictionary _orderCache = []; - // TODO: Swap this ConcurrentDictionary for the pooled dictionary variant highlighted in - // Pow2MontgomeryModBenchmarks once order warmups reuse deterministic divisor-cycle snapshots; the pooled approach removed - // the locking overhead when scanning p >= 138M in those benchmarks. - private readonly MersenneNumberIncrementalGpuTester _incrementalGpuTester = new(kernelType, useGpuOrder); - private readonly MersenneNumberIncrementalCpuTester _incrementalCpuTester = new(kernelType); - private readonly MersenneNumberOrderGpuTester _orderGpuTester = new(kernelType, useGpuOrder); - 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; - - // 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. - - private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; - - public void WarmUpOrdersHybrid(ulong exponent, ulong limit = 5_000_000UL) - { - bool cacheEnabled = useOrderCache; - if (!cacheEnabled) - { - return; - } - - UInt128 twoP = (UInt128)exponent << 1; // 2 * p - LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - if (!useGpuOrder) - { - UInt128 k = UInt128.One; - UInt128 q; - ulong remainder; - for (; k <= limit; k++) - { - q = twoP * k + UInt128.One; - remainder = q.Mod10(); - bool shouldCheck = lastDigit == LastDigit.Seven - ? remainder == 7UL || remainder == 9UL || remainder == 3UL - : remainder == 1UL || remainder == 3UL || remainder == 9UL; - - if (shouldCheck) - { - ulong mod8 = q.Mod8(); - if (mod8 != 1UL && mod8 != 7UL) - { - continue; - } - else if (q.Mod3() == 0UL || q.Mod5() == 0UL) - { - continue; - } - } - - if (!shouldCheck || !cacheEnabled || (cacheEnabled && _orderCache.ContainsKey(q))) - { - continue; - } - - ulong ord = q.CalculateOrderHybrid(gpu); - _orderCache[q] = ord; - } - - PrimeOrderCalculatorAccelerator.Return(gpu); - return; - } - - ArrayPool uInt128Pool = ThreadStaticPools.UInt128Pool; - UInt128[] qsBuffer = uInt128Pool.Rent((int)limit); - int idx = 0; - UInt128 k2 = 1UL; - UInt128 q2; - ulong remainder2; - for (; k2 <= limit; k2++) - { - q2 = twoP * k2 + 1UL; - remainder2 = q2.Mod10(); - bool shouldCheck2 = lastDigit == LastDigit.Seven - ? remainder2 == 7UL || remainder2 == 9UL || remainder2 == 3UL - : remainder2 == 1UL || remainder2 == 3UL || remainder2 == 9UL; - - if (shouldCheck2) - { - ulong mod8 = q2.Mod8(); - if (mod8 != 1UL && mod8 != 7UL) - { - shouldCheck2 = false; - } - else if (q2.Mod3() == 0UL || q2.Mod5() == 0UL) - { - shouldCheck2 = false; - } - } - - if (!shouldCheck2 || (cacheEnabled && _orderCache.ContainsKey(q2))) - { - continue; - } - - qsBuffer[idx++] = q2; - } - - if (idx == 0) - { - PrimeOrderCalculatorAccelerator.Return(gpu); - ThreadStaticPools.UInt128Pool.Return(qsBuffer); - return; - } - - UInt128[] qs = qsBuffer; - var acceleratorIndex = AcceleratorPool.Shared.Rent(); - var accelerator = _accelerators[acceleratorIndex]; - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - var orderKernel = GpuKernelPool.GetOrAddKernels(acceleratorIndex, stream, KernelType.OrderKernelScan).Order!; - - int batchSize = GpuConstants.ScanBatchSize; - ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL; - int offset = 0; - var gpuUInt128Pool = ThreadStaticPools.GpuUInt128Pool; - MemoryBuffer1D? qBuffer; - MemoryBuffer1D? orderBuffer; - - while (offset < idx) - { - int count = Math.Min(batchSize, idx - offset); - qBuffer = accelerator.Allocate1D(count); - orderBuffer = accelerator.Allocate1D(count); - - var tmp = gpuUInt128Pool.Rent(count); - for (int i = 0; i < count; i++) - { - tmp[i] = (GpuUInt128)qs[offset + i]; - } - - qBuffer.View.CopyFromCPU(stream, ref tmp[0], count); - gpuUInt128Pool.Return(tmp); - orderKernel(stream, count, exponent, divMul, qBuffer.View, orderBuffer.View); - stream.Synchronize(); - - ulong[] orders = orderBuffer.GetAsArray1D(); - for (int i = 0; i < count; i++) - { - ulong order = orders[i]; - if (order == 0UL) - { - order = qs[offset + i].CalculateOrderHybrid(gpu); - } - - if (useOrderCache) - { - _orderCache[qs[offset + i]] = order; - } - } - - offset += count; - orderBuffer.Dispose(); - qBuffer.Dispose(); - } - - AcceleratorStreamPool.Return(acceleratorIndex,stream); - PrimeOrderCalculatorAccelerator.Return(gpu); - ThreadStaticPools.UInt128Pool.Return(qs); - } - - public void WarmUpOrdersCpu(ulong exponent, ulong limit = 5_000_000UL) - { - bool cacheEnabled = useOrderCache; - if (!cacheEnabled) - { - return; - } - - UInt128 twoP = (UInt128)exponent << 1; // 2 * p - LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - // CPU path: compute orders directly and cache - if (!useGpuOrder) - { - UInt128 k = UInt128.One; - UInt128 q; - ulong remainder; - for (; k <= limit; k++) - { - q = twoP * k + UInt128.One; - remainder = q.Mod10(); - // TODO: Replace these repeated Mod10/Mod8/Mod3/Mod5 calls with a residue automaton walk so - // order warmups reuse the optimized cycle stepping validated in the residue benchmarks. - bool shouldCheck = lastDigit == LastDigit.Seven - ? remainder == 7UL || remainder == 9UL || remainder == 3UL - : remainder == 1UL || remainder == 3UL || remainder == 9UL; - - // Early modular filters for odd order and small-prime rejections - if (shouldCheck) - { - ulong mod8 = q.Mod8(); - if (mod8 != 1UL && mod8 != 7UL) - { - continue; - } - else if (q.Mod3() == 0UL || q.Mod5() == 0UL) - { - continue; - } - } - - if (!shouldCheck || !cacheEnabled || (cacheEnabled && _orderCache.ContainsKey(q))) - { - continue; - } - - ulong ord = q.CalculateOrderCpu(); - _orderCache[q] = ord; - } - - return; - } - - ArrayPool uInt128Pool = ThreadStaticPools.UInt128Pool; - UInt128[] qsBuffer = uInt128Pool.Rent((int)limit); - int idx = 0; - UInt128 k2 = 1UL; - UInt128 q2; - ulong remainder2; - for (; k2 <= limit; k2++) - { - q2 = twoP * k2 + 1UL; - remainder2 = q2.Mod10(); - // TODO: Reuse the same residue-automaton fast path here so the GPU warmup staging avoids `%` and - // branches the benchmarks showed slower than the tracked stepping helper. - bool shouldCheck2 = lastDigit == LastDigit.Seven - ? remainder2 == 7UL || remainder2 == 9UL || remainder2 == 3UL - : remainder2 == 1UL || remainder2 == 3UL || remainder2 == 9UL; - - // Early modular filters for odd order and small-prime rejections - if (shouldCheck2) - { - ulong mod8 = q2.Mod8(); - if (mod8 != 1UL && mod8 != 7UL) - { - shouldCheck2 = false; - } - else if (q2.Mod3() == 0UL || q2.Mod5() == 0UL) - { - shouldCheck2 = false; - } - } - - if (!shouldCheck2 || (cacheEnabled && _orderCache.ContainsKey(q2))) - { - continue; - } - - qsBuffer[idx++] = q2; - } - - if (idx == 0) - { - ThreadStaticPools.UInt128Pool.Return(qsBuffer); - return; - } - - UInt128[] qs = qsBuffer; - - var acceleratorIndex = AcceleratorPool.Shared.Rent(); - var accelerator = _accelerators[acceleratorIndex]; - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - var orderKernel = GpuKernelPool.GetOrAddKernels(acceleratorIndex, stream, KernelType.OrderKernelScan).Order!; - - // Guard long-running kernels by chunking the warm-up across batches. - // Reuse the same ScanBatchSize knob used for scanning. - int batchSize = GpuConstants.ScanBatchSize; - ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL; - int offset = 0; - var gpuUInt128Pool = ThreadStaticPools.GpuUInt128Pool; - MemoryBuffer1D? qBuffer; - MemoryBuffer1D? orderBuffer; - - while (offset < idx) - { - int count = Math.Min(batchSize, idx - offset); - qBuffer = accelerator.Allocate1D(count); - orderBuffer = accelerator.Allocate1D(count); - - // Convert to device-friendly GpuUInt128 on the fly - var tmp = gpuUInt128Pool.Rent(count); - for (int i = 0; i < count; i++) - { - tmp[i] = (GpuUInt128)qs[offset + i]; - } - - // Copy only the populated elements to the device buffer to avoid overrun on pooled arrays. - qBuffer.View.CopyFromCPU(stream, ref tmp[0], count); - gpuUInt128Pool.Return(tmp); - orderKernel(stream, count, exponent, divMul, qBuffer.View, orderBuffer.View); - stream.Synchronize(); - - // TODO: Change this to .CopyToCPU and use pre-allocated buffer with fixed, max capacity, using stream. Synchronize on the stream after copy. - ulong[] orders = orderBuffer.GetAsArray1D(); - for (int i = 0; i < count; i++) - { - ulong order = orders[i]; - if (order == 0UL) - { - order = qs[offset + i].CalculateOrderCpu(); - } - - if (useOrderCache) - { - _orderCache[qs[offset + i]] = order; - } - } - - offset += count; - orderBuffer.Dispose(); - qBuffer.Dispose(); - } - - AcceleratorStreamPool.Return(acceleratorIndex,stream); - uInt128Pool.Return(qs); - } - - public void WarmUpOrdersGpu(ulong exponent, ulong limit = 5_000_000UL) - { - bool cacheEnabled = useOrderCache; - if (!cacheEnabled) - { - return; - } - - UInt128 twoP = (UInt128)exponent << 1; // 2 * p - LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - // CPU path: compute orders directly and cache - var gpu = PrimeOrderCalculatorAccelerator.Rent(1); - if (!useGpuOrder) - { - UInt128 k = UInt128.One; - UInt128 q; - ulong remainder; - for (; k <= limit; k++) - { - q = twoP * k + UInt128.One; - remainder = q.Mod10(); - // TODO: Replace these repeated Mod10/Mod8/Mod3/Mod5 calls with a residue automaton walk so - // order warmups reuse the optimized cycle stepping validated in the residue benchmarks. - bool shouldCheck = lastDigit == LastDigit.Seven - ? remainder == 7UL || remainder == 9UL || remainder == 3UL - : remainder == 1UL || remainder == 3UL || remainder == 9UL; - - // Early modular filters for odd order and small-prime rejections - if (shouldCheck) - { - ulong mod8 = q.Mod8(); - if (mod8 != 1UL && mod8 != 7UL) - { - continue; - } - else if (q.Mod3() == 0UL || q.Mod5() == 0UL) - { - continue; - } - } - - if (!shouldCheck || !cacheEnabled || (cacheEnabled && _orderCache.ContainsKey(q))) - { - continue; - } - - ulong ord = q.CalculateOrderGpu(gpu); - _orderCache[q] = ord; - } - - PrimeOrderCalculatorAccelerator.Return(gpu); - return; - } - - ArrayPool uInt128Pool = ThreadStaticPools.UInt128Pool; - UInt128[] qsBuffer = uInt128Pool.Rent((int)limit); - int idx = 0; - UInt128 k2 = 1UL; - UInt128 q2; - ulong remainder2; - for (; k2 <= limit; k2++) - { - q2 = twoP * k2 + 1UL; - remainder2 = q2.Mod10(); - // TODO: Reuse the same residue-automaton fast path here so the GPU warmup staging avoids `%` and - // branches the benchmarks showed slower than the tracked stepping helper. - bool shouldCheck2 = lastDigit == LastDigit.Seven - ? remainder2 == 7UL || remainder2 == 9UL || remainder2 == 3UL - : remainder2 == 1UL || remainder2 == 3UL || remainder2 == 9UL; - - // Early modular filters for odd order and small-prime rejections - if (shouldCheck2) - { - ulong mod8 = q2.Mod8(); - if (mod8 != 1UL && mod8 != 7UL) - { - shouldCheck2 = false; - } - else if (q2.Mod3() == 0UL || q2.Mod5() == 0UL) - { - shouldCheck2 = false; - } - } - - if (!shouldCheck2 || (cacheEnabled && _orderCache.ContainsKey(q2))) - { - continue; - } - - qsBuffer[idx++] = q2; - } - - if (idx == 0) - { - PrimeOrderCalculatorAccelerator.Return(gpu); - ThreadStaticPools.UInt128Pool.Return(qsBuffer); - return; - } - - UInt128[] qs = qsBuffer; - - // GpuPrimeWorkLimiter.Acquire(); - var acceleratorIndex = AcceleratorPool.Shared.Rent(); - var accelerator = _accelerators[acceleratorIndex]; - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - var orderKernel = GpuKernelPool.GetOrAddKernels(acceleratorIndex, stream, KernelType.OrderKernelScan).Order!; - - // Guard long-running kernels by chunking the warm-up across batches. - // Reuse the same ScanBatchSize knob used for scanning. - int batchSize = GpuConstants.ScanBatchSize; - ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL; - int offset = 0; - var gpuUInt128Pool = ThreadStaticPools.GpuUInt128Pool; - MemoryBuffer1D? qBuffer; - MemoryBuffer1D? orderBuffer; - - while (offset < idx) - { - int count = Math.Min(batchSize, idx - offset); - qBuffer = accelerator.Allocate1D(count); - orderBuffer = accelerator.Allocate1D(count); - - // Convert to device-friendly GpuUInt128 on the fly - var tmp = gpuUInt128Pool.Rent(count); - for (int i = 0; i < count; i++) - { - tmp[i] = (GpuUInt128)qs[offset + i]; - } - - // Copy only the populated elements to the device buffer to avoid overrun on pooled arrays. - qBuffer.View.CopyFromCPU(stream, ref tmp[0], count); - gpuUInt128Pool.Return(tmp); - orderKernel(stream, count, exponent, divMul, qBuffer.View, orderBuffer.View); - stream.Synchronize(); - - // TODO: Change this to .CopyToCPU and use pre-allocated buffer with fixed, max capacity, using stream. Synchronize on the stream after copy. - ulong[] orders = orderBuffer.GetAsArray1D(); - for (int i = 0; i < count; i++) - { - ulong order = orders[i]; - if (order == 0UL) - { - order = qs[offset + i].CalculateOrderGpu(gpu); - } - - if (useOrderCache) - { - _orderCache[qs[offset + i]] = order; - } - } - - offset += count; - orderBuffer.Dispose(); - qBuffer.Dispose(); - } - - AcceleratorStreamPool.Return(acceleratorIndex,stream); - PrimeOrderCalculatorAccelerator.Return(gpu); - // GpuPrimeWorkLimiter.Release(); - uInt128Pool.Return(qs); - } - - public bool IsMersennePrime(PrimeOrderCalculatorAccelerator gpu, 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) - { - // TODO: Replace this `% 3` check with ULongExtensions.Mod3 to align the early rejection with - // the benchmarked bitmask helper instead of generic modulo for CPU workloads. - return false; - } - - if ((exponent & 3UL) == 1UL && _sharesFactorWithExponentMinusOne(gpu, exponent)) - { - 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; - LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - - // In residue mode the residue scan is the final primality check. - if (_useResidue) - { - if (_maxK < maxK) - { - maxK = _maxK; - } - - if (useGpuScan) - { - _residueGpuTester!.Scan(gpu, exponent, twoP, lastDigit, maxK, ref prePrime); - } - else - { - _residueCpuTester!.Scan(exponent, twoP, lastDigit, maxK, ref prePrime); - } - - return prePrime; - } - - if (!_useIncremental) - { - // Optional prefilter before Lucas–Lehmer using the same GPU/CPU scan path. - // Scan a limited range of k to quickly reject obvious cases. - UInt128 twoPPre = (UInt128)exponent << 1; - UInt128 maxKPre = (UInt128)PreLucasPreScanK; - LastDigit lastDigitPre = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; - if (useGpuScan) - { - _incrementalGpuTester.Scan(exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime); - } - else - { - _incrementalCpuTester.Scan(gpu, exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime); - } - - if (!prePrime) - { - return false; - } - - return useGpuLucas - ? _lucasLehmerGpuTester.IsPrime(gpu, exponent, runOnGpu: true) - : _lucasLehmerCpuTester.IsPrime(gpu, exponent); - } - - bool isPrime = true; - - if (useOrder) - { - if (useGpuOrder) - { - _orderGpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); - } - else - { - _orderCpuTester.Scan(gpu, exponent, twoP, lastDigit, maxK, ref isPrime); - } - } - else - { - if (useGpuScan) - { - _incrementalGpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); - } - else - { - _incrementalCpuTester.Scan(gpu, exponent, twoP, lastDigit, maxK, ref isPrime); - } - } - - return isPrime; - } - - // Limit of k to scan in the pre-Lucas prescan - public const int PreLucasPreScanK = 5_000_000; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool DividesMersenne(ulong exponent, ulong prime) - { - return 2UL.ModPow64(exponent, prime) == 1UL; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static void LegacyOrderKernel(Index1D index, ulong exponent, ArrayView qs, ArrayView orders) - { - UInt128 q = qs[index]; - UInt128 pow = exponent.PowMod(q); - orders[index] = pow == 1UL ? exponent : 0UL; - } -} diff --git a/PerfectNumbers.Core/MersenneNumberTesterTemplate.cs b/PerfectNumbers.Core/MersenneNumberTesterTemplate.cs new file mode 100644 index 0000000..e424403 --- /dev/null +++ b/PerfectNumbers.Core/MersenneNumberTesterTemplate.cs @@ -0,0 +1,352 @@ +using System.Runtime.CompilerServices; +using System.Collections.Concurrent; +using System.Buffers; +using ILGPU; +using ILGPU.Runtime; +using PerfectNumbers.Core.Cpu; +using PerfectNumbers.Core.Gpu; +using UInt128 = System.UInt128; +using PerfectNumbers.Core.Gpu.Accelerators; + +namespace PerfectNumbers.Core; + +public enum GpuKernelType +{ + Incremental, + Pow2Mod, +} + +public enum CalculationMethod +{ + ByDivisor, + Divisor, + Incremental, + LucasLehmer, + Pow2Mod, + Residue, +} + +[DeviceDependentTemplate(typeof(ComputationDevice))] +public sealed class MersenneNumberTesterTemplate( + bool useIncremental = true, + bool useOrderCache = false, + GpuKernelType kernelType = GpuKernelType.Incremental, + bool useOrder = false, + bool useGpuLucas = true, + bool useGpuScan = true, + bool useGpuOrder = true, + bool useResidue = true, + ulong maxK = 5_000_000UL, + ComputationDevice orderDevice = ComputationDevice.Gpu) +{ + private readonly bool _useResidue = useResidue; + private readonly bool _useIncremental = useIncremental && !useResidue; + private readonly ulong _maxK = maxK; + private readonly Func _sharesFactorWithExponentMinusOne = orderDevice switch + { + ComputationDevice.Hybrid => static (gpu, exponent) => exponent.SharesFactorWithExponentMinusOneHybrid(gpu), + _ when useGpuScan => static (gpu, exponent) => exponent.SharesFactorWithExponentMinusOneGpu(gpu), + _ => static (_, exponent) => exponent.SharesFactorWithExponentMinusOneCpu(), + }; + + private static readonly ConcurrentDictionary _orderCache = []; + // TODO: Swap this ConcurrentDictionary for the pooled dictionary variant highlighted in + // Pow2MontgomeryModBenchmarks once order warmups reuse deterministic divisor-cycle snapshots; the pooled approach removed + // the locking overhead when scanning p >= 138M in those benchmarks. + private readonly MersenneNumberIncrementalGpuTester _incrementalGpuTester = new(kernelType, useGpuOrder); + private readonly MersenneNumberIncrementalCpuTester _incrementalCpuTester = new(kernelType); + private readonly MersenneNumberOrderGpuTester _orderGpuTester = new(kernelType, useGpuOrder); + 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; + + // 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. + + private static readonly Accelerator[] _accelerators = AcceleratorPool.Shared.Accelerators; + + public void WarmUpOrders(ulong exponent, ulong limit = 5_000_000UL) + { + if (!useOrderCache) + { + return; + } + + UInt128 twoP = (UInt128)exponent << 1; // 2 * p + bool lastIsSeven = (exponent & 3UL) == 3UL; +#if DEVICE_GPU || DEVICE_HYBRID + var gpu = PrimeOrderCalculatorAccelerator.Rent(1); +#endif + if (!useGpuOrder) + { + UInt128 k = UInt128.One; + for (; k <= limit; k++) + { + UInt128 q = twoP * k + UInt128.One; + ulong remainder = q.Mod10(); + bool shouldCheck = lastIsSeven + ? remainder == 7UL || remainder == 9UL || remainder == 3UL + : remainder == 1UL || remainder == 3UL || remainder == 9UL; + + if (shouldCheck) + { + ulong mod8 = q.Mod8(); + if (mod8 != 1UL && mod8 != 7UL) + { + continue; + } + + if (q.Mod3() == 0UL || q.Mod5() == 0UL) + { + continue; + } + } + + if (!shouldCheck || _orderCache.ContainsKey(q)) + { + continue; + } + + _orderCache[q] = +#if DEVICE_GPU + q.CalculateOrderGpu(gpu); +#elif DEVICE_HYBRID + q.CalculateOrderHybrid(gpu); +#else + q.CalculateOrderCpu(); +#endif + } +#if DEVICE_GPU || DEVICE_HYBRID + PrimeOrderCalculatorAccelerator.Return(gpu); +#endif + return; + } + + ArrayPool uInt128Pool = ThreadStaticPools.UInt128Pool; + UInt128[] qsBuffer = uInt128Pool.Rent((int)limit); + int idx = 0; + UInt128 k2 = 1UL; + for (; k2 <= limit; k2++) + { + UInt128 q2 = twoP * k2 + 1UL; + ulong remainder2 = q2.Mod10(); + bool shouldCheck2 = lastIsSeven + ? remainder2 == 7UL || remainder2 == 9UL || remainder2 == 3UL + : remainder2 == 1UL || remainder2 == 3UL || remainder2 == 9UL; + + if (shouldCheck2) + { + ulong mod8 = q2.Mod8(); + if (mod8 != 1UL && mod8 != 7UL) + { + shouldCheck2 = false; + } + else if (q2.Mod3() == 0UL || q2.Mod5() == 0UL) + { + shouldCheck2 = false; + } + } + + if (!shouldCheck2 || _orderCache.ContainsKey(q2)) + { + continue; + } + + qsBuffer[idx++] = q2; + } + + if (idx == 0) + { +#if DEVICE_GPU || DEVICE_HYBRID + PrimeOrderCalculatorAccelerator.Return(gpu); +#endif + ThreadStaticPools.UInt128Pool.Return(qsBuffer); + return; + } + + UInt128[] qs = qsBuffer; + + var acceleratorIndex = AcceleratorPool.Shared.Rent(); + var accelerator = _accelerators[acceleratorIndex]; + var stream = AcceleratorStreamPool.Rent(acceleratorIndex); + var orderKernel = GpuKernelPool.GetOrAddKernels(acceleratorIndex, stream, KernelType.OrderKernelScan).Order!; + + int batchSize = GpuConstants.ScanBatchSize; + ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL; + int offset = 0; + var gpuUInt128Pool = ThreadStaticPools.GpuUInt128Pool; + MemoryBuffer1D? qBuffer; + MemoryBuffer1D? orderBuffer; + + while (offset < idx) + { + int count = Math.Min(batchSize, idx - offset); + qBuffer = accelerator.Allocate1D(count); + orderBuffer = accelerator.Allocate1D(count); + + var tmp = gpuUInt128Pool.Rent(count); + for (int i = 0; i < count; i++) + { + tmp[i] = (GpuUInt128)qs[offset + i]; + } + + qBuffer.View.CopyFromCPU(stream, ref tmp[0], count); + gpuUInt128Pool.Return(tmp); + orderKernel(stream, count, exponent, divMul, qBuffer.View, orderBuffer.View); + stream.Synchronize(); + + ulong[] orders = orderBuffer.GetAsArray1D(); + for (int i = 0; i < count; i++) + { + ulong order = orders[i]; + if (order == 0UL) + { + order = +#if DEVICE_GPU + qs[offset + i].CalculateOrderGpu(gpu); +#elif DEVICE_HYBRID + qs[offset + i].CalculateOrderHybrid(gpu); +#else + qs[offset + i].CalculateOrderCpu(); +#endif + } + + if (useOrderCache) + { + _orderCache[qs[offset + i]] = order; + } + } + + offset += count; + orderBuffer.Dispose(); + qBuffer.Dispose(); + } + + AcceleratorStreamPool.Return(acceleratorIndex,stream); +#if DEVICE_GPU || DEVICE_HYBRID + PrimeOrderCalculatorAccelerator.Return(gpu); +#endif + uInt128Pool.Return(qs); + } + + public bool IsMersennePrime(PrimeOrderCalculatorAccelerator gpu, 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) + { + // TODO: Replace this `% 3` check with ULongExtensions.Mod3 to align the early rejection with + // the benchmarked bitmask helper instead of generic modulo for CPU workloads. + return false; + } + + if ((exponent & 3UL) == 1UL && _sharesFactorWithExponentMinusOne(gpu, exponent)) + { + 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; + LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + + // In residue mode the residue scan is the final primality check. + if (_useResidue) + { + if (_maxK < maxK) + { + maxK = _maxK; + } + + if (useGpuScan) + { + _residueGpuTester!.Scan(gpu, exponent, twoP, lastDigit, maxK, ref prePrime); + } + else + { + _residueCpuTester!.Scan(exponent, twoP, lastDigit, maxK, ref prePrime); + } + + return prePrime; + } + + if (!_useIncremental) + { + // Optional prefilter before Lucas–Lehmer using the same GPU/CPU scan path. + // Scan a limited range of k to quickly reject obvious cases. + UInt128 twoPPre = (UInt128)exponent << 1; + UInt128 maxKPre = (UInt128)PreLucasPreScanK; + LastDigit lastDigitPre = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + if (useGpuScan) + { + _incrementalGpuTester.Scan(exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime); + } + else + { + _incrementalCpuTester.Scan(gpu, exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime); + } + + if (!prePrime) + { + return false; + } + + return useGpuLucas + ? _lucasLehmerGpuTester.IsPrime(gpu, exponent, runOnGpu: true) + : _lucasLehmerCpuTester.IsPrime(gpu, exponent); + } + + bool isPrime = true; + + if (useOrder) + { + if (useGpuOrder) + { + _orderGpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); + } + else + { + _orderCpuTester.Scan(gpu, exponent, twoP, lastDigit, maxK, ref isPrime); + } + } + else + { + if (useGpuScan) + { + _incrementalGpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); + } + else + { + _incrementalCpuTester.Scan(gpu, exponent, twoP, lastDigit, maxK, ref isPrime); + } + } + + return isPrime; + } + + // Limit of k to scan in the pre-Lucas prescan + public const int PreLucasPreScanK = 5_000_000; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool DividesMersenne(ulong exponent, ulong prime) + { + return 2UL.ModPow64(exponent, prime) == 1UL; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void LegacyOrderKernel(Index1D index, ulong exponent, ArrayView qs, ArrayView orders) + { + UInt128 q = qs[index]; + UInt128 pow = exponent.PowMod(q); + orders[index] = pow == 1UL ? exponent : 0UL; + } +}