MersenneNumberTester converted to device dependent MersenneNumberTesterTemplate.

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<MersenneNumberTester>(() => new MersenneNumberTester(useIncremental: true, useResidue: true, maxK: 1_000UL), trackAllValues: true));
+        mersenneField.SetValue(null, new ThreadLocal<object>(() => new MersenneNumberTesterCpu(useIncremental: true, useResidue: true, maxK: 1_000UL), trackAllValues: true));
         primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
         CandidatesCalculator.OverrideResidueTrackers(new ThreadLocal<ModResidueTracker>(() => 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<MersenneNumberTester>(() => new MersenneNumberTester(
+        mersenneField.SetValue(null, new ThreadLocal<object>(() => 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<MersenneNumberTester>(() => new MersenneNumberTester(
+        mersenneField.SetValue(null, new ThreadLocal<object>(() => 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<MersenneNumberTester>(() => new MersenneNumberTester(
+        mersenneField.SetValue(null, new ThreadLocal<object>(() => 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<MersenneNumberTester>(() => new MersenneNumberTester(
+        mersenneField.SetValue(null, new ThreadLocal<object>(() => 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<MersenneNumberTester>(() => new MersenneNumberTester(
+        mersenneField.SetValue(null, new ThreadLocal<object>(() => 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<MersenneNumberTester>(() => new MersenneNumberTester(useIncremental: true, useResidue: false, maxK: 1_000UL), trackAllValues: true));
+        mersenneField.SetValue(null, new ThreadLocal<object>(() => new MersenneNumberTesterCpu(useIncremental: true, useResidue: false, maxK: 1_000UL), trackAllValues: true));
         primeField.SetValue(null, new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true));
         CandidatesCalculator.OverrideResidueTrackers(new ThreadLocal<ModResidueTracker>(() => new ModResidueTracker(ResidueModel.Identity, 2UL, true), trackAllValues: true));
 
@@ -547,3 +547,4 @@ public class ProgramTests
     }
 }
 
+

EvenPerfectBitScanner/Program.cs

@@ -14,7 +14,7 @@ namespace EvenPerfectBitScanner;
 internal static class Program
 {
 	// private static ThreadLocal<PrimeTester> PrimeTesters = null!;
-	private static ThreadLocal<MersenneNumberTester> MersenneTesters = null!;
+private static ThreadLocal<object>? 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,10 +218,12 @@ internal static class Program
 			// Initialize per-thread p residue tracker (Identity model) at currentP
 			if (!useDivisor && !_cliArguments.UseByDivisor)
 			{
-				MersenneTesters = new ThreadLocal<MersenneNumberTester>(() =>
+				MersenneTesters = new ThreadLocal<object>(() =>
 				{
 					// ProcessEightBitWindows windowed pow2 ladder is the default kernel.
-					var tester = new MersenneNumberTester(
+					object tester = orderDevice switch
+					{
+						ComputationDevice.Gpu => new MersenneNumberTesterGpu(
 							useIncremental: !useLucas,
 							useOrderCache: false,
 							kernelType: kernelType,
@@ -227,21 +233,44 @@ internal static class Program
 							useGpuOrder: orderOnGpu || orderOnHybrid,
 							useResidue: useResidue,
 							maxK: residueKMax,
-												orderDevice: orderDevice);
+							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;
 	}
 }
+

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);
+		var tester = new MersenneNumberTesterCpu(useOrderCache: true, useGpuOrder: true);
-		tester.WarmUpOrdersCpu(p, 1_000UL);
+		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);
+		var tester = new MersenneNumberTesterGpu(useOrderCache: true, useGpuOrder: true);
-		tester.WarmUpOrdersGpu(p, 1_000UL);
+		tester.WarmUpOrders(p, 1_000UL);
 		var gpu = PrimeOrderCalculatorAccelerator.Rent(1);
 		try
 		{
@@ -84,3 +84,5 @@ public class MersenneNumberTesterTests
 	}
 }
 
+
+

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
 	}
 }
 
+
+

PerfectNumbers.Core/MersenneNumberTester.cs

@@ -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<PrimeOrderCalculatorAccelerator, ulong, bool> _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<UInt128, ulong> _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<UInt128> 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<GpuUInt128, Stride1D.Dense>? qBuffer;
-		MemoryBuffer1D<ulong, Stride1D.Dense>? orderBuffer;
-		
-		while (offset < idx)
-		{
-			int count = Math.Min(batchSize, idx - offset);
-			qBuffer = accelerator.Allocate1D<GpuUInt128>(count);
-			orderBuffer = accelerator.Allocate1D<ulong>(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<UInt128> 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<GpuUInt128, Stride1D.Dense>? qBuffer;
-		MemoryBuffer1D<ulong, Stride1D.Dense>? orderBuffer;
-		
-		while (offset < idx)
-		{
-			int count = Math.Min(batchSize, idx - offset);
-			qBuffer = accelerator.Allocate1D<GpuUInt128>(count);
-			orderBuffer = accelerator.Allocate1D<ulong>(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<UInt128> 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<GpuUInt128, Stride1D.Dense>? qBuffer;
-		MemoryBuffer1D<ulong, Stride1D.Dense>? orderBuffer;
-		
-		while (offset < idx)
-		{
-			int count = Math.Min(batchSize, idx - offset);
-			qBuffer = accelerator.Allocate1D<GpuUInt128>(count);
-			orderBuffer = accelerator.Allocate1D<ulong>(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<UInt128> qs, ArrayView<ulong> orders)
-	{
-		UInt128 q = qs[index];
-		UInt128 pow = exponent.PowMod(q);
-		orders[index] = pow == 1UL ? exponent : 0UL;
-	}
-}

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<PrimeOrderCalculatorAccelerator, ulong, bool> _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<UInt128, ulong> _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<UInt128> 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<GpuUInt128, Stride1D.Dense>? qBuffer;
+		MemoryBuffer1D<ulong, Stride1D.Dense>? orderBuffer;
+		
+		while (offset < idx)
+		{
+			int count = Math.Min(batchSize, idx - offset);
+			qBuffer = accelerator.Allocate1D<GpuUInt128>(count);
+			orderBuffer = accelerator.Allocate1D<ulong>(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<UInt128> qs, ArrayView<ulong> orders)
+	{
+		UInt128 q = qs[index];
+		UInt128 pow = exponent.PowMod(q);
+		orders[index] = pow == 1UL ? exponent : 0UL;
+	}
+}