Initial commit.

PerfectNumbers.Core/AlphaCache.cs

@@ -0,0 +1,41 @@
+using PeterO.Numbers;
+
+namespace PerfectNumbers.Core;
+
+public sealed class AlphaCache
+{
+    private readonly Dictionary<(EInteger P, int K), AlphaValues> _cache = new();
+
+    public AlphaValues Get(EInteger p, int k)
+    {
+        if (!_cache.TryGetValue((p, k), out var values))
+        {
+            ERational alphaPValue = AlphaCalculations.ComputeAlphaP(p, k);
+            ERational alphaMValue = RationalNumbers.Two.Divide(alphaPValue);
+            values = new AlphaValues(alphaMValue, alphaPValue);
+            _cache[(p, k)] = values;
+        }
+
+        return values;
+    }
+
+    public void Clear()
+    {
+        _cache.Clear();
+    }
+
+    public int Count => _cache.Count;
+
+    public readonly struct AlphaValues
+    {
+        public readonly ERational AlphaM;
+        public readonly ERational AlphaP;
+
+        public AlphaValues(ERational alphaM, ERational alphaP)
+        {
+            AlphaM = alphaM;
+            AlphaP = alphaP;
+        }
+    }
+}
+

PerfectNumbers.Core/AlphaCalculations.cs

@@ -0,0 +1,29 @@
+using PeterO.Numbers;
+
+namespace PerfectNumbers.Core;
+
+public static class AlphaCalculations
+{
+    public static ERational ComputeAlphaP(EInteger p, int k)
+    {
+        int exponent = checked(4 * k + 1);
+        EInteger y = PowerCache.Get(p, exponent);
+        return ERational.Create(y * p - 1, y * (p - 1));
+    }
+
+    public static ERational ComputeAlphaFactor(EInteger q, int a2)
+    {
+        EInteger num = PowerCache.Get(q, a2 + 1) - EInteger.One;
+        EInteger den = PowerCache.Get(q, a2) * (q - EInteger.One);
+        return ERational.Create(num, den);
+    }
+
+    public static ERational ComputeAlphaM(ReadOnlySpan<EInteger> primes, int maxR, int a2)
+    {
+        ERational result = ERational.One;
+        for (int i = 0; i < maxR; i++)
+            result = result.Multiply(ComputeAlphaFactor(primes[i], a2));
+        return result;
+    }
+}
+

PerfectNumbers.Core/AlphaMScannerCandidateStatus.cs

@@ -0,0 +1,13 @@
+namespace PerfectNumbers.Core;
+
+public enum AlphaMScannerCandidateStatus
+{
+    UNKNOWN,
+    PRODUCT_ABOVE_TARGET,
+    PRODUCT_BELOW_TARGET,
+    ALPHA_P_AND_ALPHA_M_OUT_OF_RANGE,
+    ALPHA_P_OUT_OF_RANGE,
+    ALPHA_M_OUT_OF_RANGE,
+    OK,
+}
+

PerfectNumbers.Core/Cpu/CpuConstants.cs

@@ -0,0 +1,9 @@
+namespace PerfectNumbers.Core.Cpu;
+
+public static class CpuConstants
+{
+	// Allowed q mod 10 sets depending on last digit of M_p
+	public const int LastSevenMask10 = (1 << 7) | (1 << 9);
+	// For M_p ending with 1, restrict to {1,3,9}; 7 considered only for q == 7 (not possible for large p)
+	public const int LastOneMask10 = (1 << 1) | (1 << 3) | (1 << 9);
+}

PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs

@@ -0,0 +1,83 @@
+namespace PerfectNumbers.Core.Cpu;
+
+public class MersenneNumberIncrementalCpuTester(GpuKernelType kernelType)
+{
+	private static readonly MersenneDivisorCycles _cycles = MersenneDivisorCycles.Shared;
+	private readonly GpuKernelType _kernelType = kernelType;
+
+	public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
+	{
+		UInt128 k, one = k = UInt128.One;
+		ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
+		// Use residue automaton always on CPU
+		var auto = new MersenneResidueAutomaton(exponent);
+
+		bool reject, shouldCheck;
+		UInt128 q = auto.CurrentQ(),
+				divPow, halfPow, phi, qCycle;
+
+		ulong div, mod8, phi64;
+
+		while (k <= maxK && Volatile.Read(ref isPrime))
+		{
+			qCycle = MersenneDivisorCycles.GetCycle(q);
+			shouldCheck = lastIsSeven
+				? (auto.Mod10R != 5UL)
+				: (auto.Mod10R != 5UL);
+			if (shouldCheck)
+			{
+				mod8 = auto.Mod8R;
+				if (mod8 != 1UL && mod8 != 7UL)
+				{
+					shouldCheck = false;
+				}
+				else if (auto.Mod3R == 0UL || auto.Mod5R == 0UL)
+				{
+					shouldCheck = false;
+				}
+			}
+
+			if (shouldCheck)
+			{
+				if (_kernelType == GpuKernelType.Pow2Mod)
+				{
+					reject = exponent.PowModWithCycle(q, qCycle) == UInt128.One;
+				}
+				else
+				{
+					phi = q - one;
+					reject = false;
+					if (phi <= ulong.MaxValue)
+					{
+						phi64 = (ulong)phi;
+						if (phi64.PowModWithCycle(q, qCycle) == one)
+						{
+							halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - one;
+							if (halfPow.BinaryGcd(q) == one)
+							{
+								div = phi64.FastDiv64(exponent, divMul);
+								divPow = div.PowModWithCycle(q, qCycle) - one;
+								if (divPow.BinaryGcd(q) == one)
+								{
+									reject = true;
+								}
+							}
+						}
+					}
+				}
+
+				if (reject && q.IsPrimeCandidate())
+				{
+					Volatile.Write(ref isPrime, false);
+					break;
+				}
+			}
+
+			k += 1UL;
+			auto.Next();
+			q = auto.CurrentQ();
+		}
+
+		return;
+	}
+}

PerfectNumbers.Core/Cpu/MersenneNumberLucasLehmerCpuTester.cs

@@ -0,0 +1,43 @@
+using PeterO.Numbers;
+
+namespace PerfectNumbers.Core;
+
+public class MersenneNumberLucasLehmerCpuTester
+{
+    public bool IsPrime(ulong exponent)
+	{
+		// Early rejections aligned with incremental/order sieves, but safe for small p:
+		// - If 3 | p and p != 3, then 7 | M_p -> composite.
+		// - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast.
+		// - If p is divisible by divisors specific to numbers ending with 1 or 7, reject fast.
+		if (
+			((exponent % 3UL) == 0UL && exponent != 3UL) ||
+			((exponent % 5UL) == 0UL && exponent != 5UL) ||
+			((exponent % 7UL) == 0UL && exponent != 7UL) ||
+			((exponent % 11UL) == 0UL && exponent != 11UL)
+		)
+		{
+			return false;
+		}
+
+		if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
+		{
+			return false;
+		}
+
+		UInt128 s = UInt128Numbers.Four, two = UInt128Numbers.Two,
+			 m = (UInt128.One << (int)exponent) - UInt128.One;
+
+		UInt128 limit = exponent - UInt128Numbers.Two;
+		for (UInt128 i = UInt128.Zero; i < limit; i++)
+		{
+			s = s.PowModWithCycle(two, MersenneDivisorCycles.GetCycle(m)) - two;
+			if (s < UInt128.Zero)
+			{
+				s += m;
+			}
+		}
+
+		return s == UInt128.Zero;
+	}
+}

PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs

@@ -0,0 +1,74 @@
+namespace PerfectNumbers.Core.Cpu;
+
+public class MersenneNumberOrderCpuTester(GpuKernelType kernelType)
+{
+	private readonly GpuKernelType _kernelType = kernelType;
+
+	public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
+	{
+		UInt128 k = 1UL;
+		var auto = new MersenneResidueAutomaton(exponent);
+		UInt128 qCycle,
+				q = auto.CurrentQ();
+
+		while (k <= maxK && Volatile.Read(ref isPrime))
+		{
+			qCycle = MersenneDivisorCycles.GetCycle(q);
+			int allowMask2 = lastIsSeven ? CpuConstants.LastSevenMask10 : CpuConstants.LastOneMask10;
+			bool shouldCheck = ((allowMask2 >> (int)auto.Mod10R) & 1) != 0;
+			if (shouldCheck)
+			{
+				ulong mod8 = auto.Mod8R;
+				if (mod8 != 1UL && mod8 != 7UL)
+				{
+					shouldCheck = false;
+				}
+				else if (auto.Mod3R == 0UL || auto.Mod5R == 0UL)
+				{
+					shouldCheck = false;
+				}
+			}
+			if (shouldCheck)
+			{
+				if (_kernelType == GpuKernelType.Pow2Mod)
+				{
+					if (exponent.PowModWithCycle(q, qCycle) == 1UL)
+					{
+						Volatile.Write(ref isPrime, false);
+						break;
+					}
+				}
+				else
+				{
+					UInt128 phi = q - 1UL;
+					if (phi <= ulong.MaxValue)
+					{
+						ulong phi64 = (ulong)phi;
+						if (phi64.PowModWithCycle(q, qCycle) == 1UL)
+						{
+							UInt128 halfPow = (phi64 >> 1).PowModWithCycle(q, qCycle) - 1UL;
+							if (halfPow.BinaryGcd(q) == 1UL)
+							{
+								ulong divMul = (ulong)((((UInt128)1 << 64) - 1UL) / exponent) + 1UL;
+								ulong div = phi64.FastDiv64(exponent, divMul);
+								UInt128 divPow = div.PowModWithCycle(q, qCycle) - 1UL;
+								if (divPow.BinaryGcd(q) == 1UL)
+								{
+									Volatile.Write(ref isPrime, false);
+									break;
+								}
+							}
+						}
+					}
+				}
+			}
+
+			k += 1UL;
+			auto.Next();
+			q = auto.CurrentQ();
+		}
+		
+		return;
+	}
+
+}

PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs

@@ -0,0 +1,169 @@
+namespace PerfectNumbers.Core.Cpu;
+
+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
+
+		_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);
+
+		// Start at k = 1: q0 = 2*p*1 + 1
+		UInt128 k = UInt128.One;
+		UInt128 q = twoP + UInt128.One;
+		// Current residues
+		q.Mod10_8_5_3(out ulong r10, out ulong r8, out ulong r5, out ulong r3);
+		// Step residues for q += 2*p
+		ulong step10 = (ulong)(twoP % 10UL) % 10UL;
+		ulong step8 = (ulong)(twoP % 8UL);
+		ulong step3 = (ulong)(twoP % 3UL);
+		ulong step5 = (ulong)(twoP % 5UL);
+
+		// Allowed last-digit sets for q depending on last digit of M_p
+		// lastIsSeven == true  => allow {7,9}
+		// lastIsSeven == false => allow {1,3,7,9}
+		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;
+
+		bool localIsPrime, isDivider;
+
+		maxK++;
+
+		ModResidueTracker tracker = _mersenneResidueTracker!;
+		while (k < maxK && (localIsPrime = Volatile.Read(ref isPrime)))
+		{
+			remaining = maxK - k;
+			if (remaining == UInt128.Zero)
+			{
+				break;
+			}
+
+			lanes = remaining > 4UL ? 4 : (int)remaining;
+
+			// 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;
+
+			// r8
+			r8_0 = r8 & 7UL;
+			r8_1 = (r8_0 + step8) & 7UL;
+			r8_2 = (r8_1 + step8) & 7UL;
+			r8_3 = (r8_2 + step8) & 7UL;
+
+			// 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;
+
+			// 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 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;
+				}
+			}
+
+			// lane 1
+			qIncremental = q + twoP;
+			if (lanes >= 2 &&
+				((allowMask >> (int)r10_1) & 1) != 0 && (r8_1 == 1UL || r8_1 == 7UL) && r3_1 != 0UL && r5_1 != 0UL &&
+				localIsPrime)
+			{
+				if (localIsPrime && tracker.MergeOrAppend(exponent, qIncremental, out isDivider) && isDivider &&
+					(qIncremental <= ulong.MaxValue ? ((ulong)qIncremental).IsPrimeCandidate() : qIncremental.IsPrimeCandidate()))
+				{
+					localIsPrime = false;
+				}
+			}
+
+			// lane 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 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;
+				}
+			}
+
+			if (!localIsPrime)
+			{
+				Volatile.Write(ref isPrime, false);
+				return;
+			}
+			else if (!Volatile.Read(ref isPrime))
+			{
+				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;
+			}
+		}
+
+		return;
+	}
+}

PerfectNumbers.Core/EulerPrimeTester.cs

@@ -0,0 +1,19 @@
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+public sealed class EulerPrimeTester
+{
+    private readonly PrimeTester _primeTester;
+
+    public EulerPrimeTester(PrimeTester primeTester)
+    {
+        _primeTester = primeTester;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public bool IsEulerPrime(ulong n)
+    {
+        return (n & 3UL) == 1UL && _primeTester.IsPrime(n, CancellationToken.None);
+    }
+}

PerfectNumbers.Core/Gpu/GpuConstants.cs

@@ -0,0 +1,7 @@
+namespace PerfectNumbers.Core.Gpu
+{
+    public static class GpuConstants
+    {
+	    public const int ScanBatchSize = 2_097_152;
+    }
+}

PerfectNumbers.Core/Gpu/GpuContextPool.cs

@@ -0,0 +1,152 @@
+using System.Collections.Concurrent;
+using ILGPU;
+using ILGPU.Runtime;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public static class GpuContextPool
+{
+	private static readonly bool PoolingEnabled = false;
+	// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
+	public static bool ForceCpu { get; set; } = false;
+
+	internal sealed class PooledContext
+	{
+		private bool disposed;
+
+		public Context Context { get; }
+		public Accelerator Accelerator { get; }
+		public bool IsCpu { get; }
+
+		public PooledContext(bool preferCpu)
+		{
+			Context = Context.CreateDefault();
+			Accelerator = Context.GetPreferredDevice(preferCpu).CreateAccelerator(Context);
+			IsCpu = Accelerator.AcceleratorType == AcceleratorType.CPU;
+			// NOTE: Avoid loading/compiling any kernel here to prevent implicit
+			// CL stream/queue creation during accelerator construction.
+			// Some OpenCL drivers are fragile when a queue is created immediately
+			// after device init. Real kernels will be JIT-loaded at first use.
+		}
+
+		public void Dispose()
+		{
+			if (disposed)
+			{
+				return;
+			}
+
+			// Ensure all cached GPU buffers for this accelerator are released.
+			NttGpuMath.ClearCaches(Accelerator);
+			// Release PrimeTester GPU state for this accelerator (kernel/device primes)
+			// Clear any per-accelerator cached resources to avoid releasing
+			// after the accelerator is destroyed.
+			PrimeTester.ClearGpuCaches(Accelerator);
+
+			Accelerator.Dispose();
+			Context.Dispose();
+			disposed = true;
+		}
+	}
+
+	private static readonly ConcurrentQueue<PooledContext> CpuPool = new();
+	private static readonly ConcurrentQueue<PooledContext> GpuPool = new();
+	private static readonly object CreationLock = new();
+
+	public static GpuContextLease Rent()
+	{
+		return RentPreferred(ForceCpu);
+	}
+
+	// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
+	public static GpuContextLease RentPreferred(bool preferCpu)
+	{
+		if (PoolingEnabled)
+		{
+			if (preferCpu && CpuPool.TryDequeue(out var cpu))
+			{
+				return new GpuContextLease(cpu, CreationLock);
+			}
+
+			if (!preferCpu && GpuPool.TryDequeue(out var gpu))
+			{
+				return new GpuContextLease(gpu, CreationLock);
+			}
+		}
+
+		// Serialize accelerator creation to avoid concurrent CL device/stream setup across threads.
+		lock (CreationLock)
+		{
+			return new GpuContextLease(new PooledContext(preferCpu), CreationLock);
+		}
+	}
+
+	public static void DisposeAll()
+	{
+		if (!PoolingEnabled)
+		{
+			return;
+		}
+
+		while (CpuPool.TryDequeue(out var cpu))
+		{
+			cpu.Dispose();
+		}
+
+		while (GpuPool.TryDequeue(out var gpu))
+		{
+			gpu.Dispose();
+		}
+	}
+
+	private static void Return(PooledContext ctx)
+	{
+		if (PoolingEnabled)
+		{
+			ctx.Accelerator.Synchronize();
+			if (ctx.IsCpu)
+			{
+				CpuPool.Enqueue(ctx);
+			}
+			else
+			{
+				GpuPool.Enqueue(ctx);
+			}
+		}
+
+		else
+		{
+			ctx.Dispose();
+		}
+	}
+
+    public struct GpuContextLease : IDisposable
+    {
+        private readonly PooledContext _ctx;
+		private bool disposed;
+
+		internal GpuContextLease(PooledContext ctx, object kernelInitLock)
+        {
+            _ctx = ctx;
+			KernelInitLock = kernelInitLock;
+		}
+
+        public Context Context => _ctx.Context;
+
+        public Accelerator Accelerator => _ctx.Accelerator;
+
+		public object KernelInitLock { get; }
+
+
+		public void Dispose()
+		{
+			if (disposed)
+			{
+				return;
+			}
+
+			Return(_ctx);
+			disposed = true;
+        }
+    }
+}

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

@@ -0,0 +1,550 @@
+using System.Collections.Concurrent;
+using ILGPU;
+using ILGPU.Runtime;
+using static PerfectNumbers.Core.Gpu.GpuContextPool;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public readonly struct ResidueAutomatonArgs
+{
+    public readonly ulong Q0M10;
+    public readonly ulong Step10;
+    public readonly ulong Q0M8;
+    public readonly ulong Step8;
+    public readonly ulong Q0M3;
+    public readonly ulong Step3;
+    public readonly ulong Q0M5;
+    public readonly ulong Step5;
+
+    public ResidueAutomatonArgs(ulong q0m10, ulong step10, ulong q0m8, ulong step8, ulong q0m3, ulong step3, ulong q0m5, ulong step5)
+    {
+        Q0M10 = q0m10; Step10 = step10; Q0M8 = q0m8; Step8 = step8; Q0M3 = q0m3; Step3 = step3; Q0M5 = q0m5; Step5 = step5;
+    }
+}
+
+public sealed class KernelContainer
+{
+	// Serializes first-time initialization of kernels/buffers per accelerator.
+	public Action<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>? Order;
+    public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+        ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>>? Incremental;
+    public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+        ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>>? Pow2Mod;
+	public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+		ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>>? IncrementalOrder;
+    public Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+        ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>>? Pow2ModOrder;
+
+    // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
+    public MemoryBuffer1D<uint, Stride1D.Dense>? SmallCycles;
+
+	public static T InitOnce<T>(ref T? slot, Func<T> factory) where T : class
+	{
+		var current = Volatile.Read(ref slot);
+		if (current is not null)
+		{
+			return current;
+		}
+
+		current = Volatile.Read(ref slot);
+		if (current is null)
+		{
+			current = factory();
+			Volatile.Write(ref slot, current);
+		}
+		return current;
+	}
+}
+
+public ref struct GpuKernelLease(IDisposable limiter, GpuContextLease gpu, KernelContainer kernels)
+{
+	private bool disposedValue;
+	private IDisposable? _limiter = limiter;
+	private GpuContextLease? _gpu = gpu;
+	private KernelContainer _kernels = kernels;
+
+	public readonly Accelerator Accelerator => _gpu?.Accelerator ?? throw new NullReferenceException("GPU context is already released");
+
+	public readonly Action<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>> OrderKernel
+	{
+		get
+		{
+			var accel = Accelerator; // avoid capturing 'this' in lambda
+			lock (_gpu!.Value.KernelInitLock)
+			{
+				return KernelContainer.InitOnce(ref _kernels.Order, () => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernelScan));
+			}
+		}
+	}
+
+    public readonly Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>> Pow2ModKernel
+	{
+		get
+		{
+			var accel = Accelerator;
+			lock (_gpu!.Value.KernelInitLock)
+			{
+				return KernelContainer.InitOnce(ref _kernels.Pow2Mod, () => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>>(Pow2ModKernelScan));
+			}
+		}
+	}
+
+    public readonly Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>> IncrementalKernel
+	{
+		get
+		{
+			var accel = Accelerator;
+			lock (_gpu!.Value.KernelInitLock)
+			{
+				return KernelContainer.InitOnce(ref _kernels.Incremental, () => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<uint, Stride1D.Dense>>(IncrementalKernelScan));
+			}
+		}
+	}
+
+	public readonly Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>> IncrementalOrderKernel
+	{
+		get
+		{
+			var accel = Accelerator;
+			lock (_gpu!.Value.KernelInitLock)
+			{
+				return KernelContainer.InitOnce(ref _kernels.IncrementalOrder, () => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>>(IncrementalOrderKernelScan));
+			}
+		}
+	}
+
+	public readonly Action<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>> Pow2ModOrderKernel
+	{
+		get
+		{
+			var accel = Accelerator;
+			lock (_gpu!.Value.KernelInitLock)
+			{
+				return KernelContainer.InitOnce(ref _kernels.Pow2ModOrder, () => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<uint, Stride1D.Dense>>(Pow2ModOrderKernelScan));
+			}
+		}
+	}
+
+	// TODO: Plumb the small-cycles device buffer into all kernels that can benefit
+	// (some already accept it). Consider a compact type (byte/ushort) for memory footprint.
+
+    private static void IncrementalKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul,
+        ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView<ulong> orders,
+        ArrayView1D<uint, Stride1D.Dense> smallCycles)
+    {
+        ulong idx = (ulong)index.X;
+		// residue automaton
+        ulong step10 = ((exponent % 10UL) << 1) % 10UL;
+        ulong r10 = q0m10 + (step10 * idx) % 10UL;
+		r10 -= (r10 >= 10UL) ? 10UL : 0UL;
+		bool shouldCheck = r10 != 5UL; // skip q ≡ 5 (mod 10); other residues handled by r8/r3/r5
+		if (shouldCheck)
+		{
+            ulong step8 = ((exponent & 7UL) << 1) & 7UL;
+            ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL;
+			if (r8 != 1UL && r8 != 7UL)
+			{
+				shouldCheck = false;
+			}
+            else
+            {
+                ulong step3 = ((exponent % 3UL) << 1) % 3UL;
+                ulong r3 = q0m3 + (step3 * (idx % 3UL)); r3 -= (r3 >= 3UL) ? 3UL : 0UL;
+                if (r3 == 0UL)
+                {
+                    shouldCheck = false;
+                }
+                else
+                {
+                    ulong step5Loc = ((exponent % 5UL) << 1) % 5UL;
+                    ulong r5 = q0m5 + (step5Loc * (idx % 5UL)); if (r5 >= 5UL) r5 -= 5UL;
+                    if (r5 == 0UL)
+                    {
+                        shouldCheck = false;
+                    }
+                }
+            }
+		}
+
+		if (!shouldCheck)
+		{
+			orders[index] = 0UL;
+			return;
+		}
+
+        GpuUInt128 k = kStart + (GpuUInt128)idx;
+        GpuUInt128 q = twoP.Mul64(k) + GpuUInt128.One;
+        
+        // Small-cycles in-kernel early rejection from device table
+        if (q.High == 0UL && q.Low <= (ulong)(smallCycles.Length - 1))
+        {
+            uint cyc = smallCycles[(int)q.Low];
+            if (cyc != 0U && (exponent % cyc) == 0UL)
+            {
+                orders[index] = exponent;
+                return;
+            }
+        }
+
+        GpuUInt128 phi = q - GpuUInt128.One;
+        if (phi.High != 0UL)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        ulong phi64 = phi.Low;
+        if (GpuUInt128.Pow2Mod(phi64, q) != GpuUInt128.One)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, q) - GpuUInt128.One;
+        if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        ulong div = phi64.FastDiv64(exponent, divMul);
+        GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, q) - GpuUInt128.One;
+        if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        orders[index] = 1UL;
+    }
+
+    private static void OrderKernelScan(Index1D index, ulong exponent, ulong divMul, ArrayView<GpuUInt128> qs, ArrayView<ulong> orders)
+    {
+        GpuUInt128 q = qs[index];
+        GpuUInt128 phi = q - GpuUInt128.One;
+        if (phi.High != 0UL)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        ulong phi64 = phi.Low;
+        GpuUInt128 pow = GpuUInt128.Pow2Mod(phi64, q);
+        if (pow != GpuUInt128.One)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, q) - GpuUInt128.One;
+        if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        ulong div = phi64.FastDiv64(exponent, divMul);
+        GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, q) - GpuUInt128.One;
+        if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        orders[index] = exponent;
+    }
+
+	private void Dispose(bool disposing)
+	{
+		if (!disposedValue)
+		{
+			if (disposing)
+			{
+				_gpu?.Dispose();
+				_gpu = null;
+
+				_limiter?.Dispose();
+				_limiter = null;
+			}
+
+			// free unmanaged resources (unmanaged objects) and override finalizer
+			// set large fields to null
+			disposedValue = true;
+		}
+	}
+
+    private static void Pow2ModKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _,
+        ulong q0m10, ulong q0m8, ulong q0m3, ulong q0m5, ArrayView<ulong> orders,
+        ArrayView1D<uint, Stride1D.Dense> smallCycles)
+    {
+        ulong idx = (ulong)index.X;
+        ulong step10 = ((exponent % 10UL) << 1) % 10UL;
+        ulong r10 = q0m10 + (step10 * idx) % 10UL; r10 -= (r10 >= 10UL) ? 10UL : 0UL;
+		bool shouldCheck = r10 != 5UL;
+		if (shouldCheck)
+		{
+            ulong step8 = ((exponent & 7UL) << 1) & 7UL;
+            ulong r8 = (q0m8 + ((step8 * idx) & 7UL)) & 7UL;
+			if (r8 != 1UL && r8 != 7UL)
+			{
+				shouldCheck = false;
+			}
+            else
+            {
+                ulong step3 = ((exponent % 3UL) << 1) % 3UL;
+                ulong r3 = q0m3 + (step3 * (idx % 3UL)); r3 -= (r3 >= 3UL) ? 3UL : 0UL;
+                if (r3 == 0UL)
+                {
+                    shouldCheck = false;
+                }
+                else
+                {
+                    ulong step5Loc = ((exponent % 5UL) << 1) % 5UL;
+                    ulong r5 = q0m5 + (step5Loc * (idx % 5UL)); if (r5 >= 5UL) r5 -= 5UL;
+                    if (r5 == 0UL)
+                    {
+                        shouldCheck = false;
+                    }
+                }
+            }
+		}
+		if (!shouldCheck)
+		{
+			orders[index] = 0UL;
+			return;
+		}
+
+        GpuUInt128 k = kStart + (GpuUInt128)idx;
+        GpuUInt128 q = twoP.Mul64(k) + GpuUInt128.One;
+        // Small-cycles in-kernel early rejection from device table
+        if (q.High == 0UL && q.Low <= (ulong)(smallCycles.Length - 1))
+        {
+            uint cyc = smallCycles[(int)q.Low];
+            if (cyc != 0U && (exponent % cyc) == 0UL)
+            {
+                orders[index] = exponent;
+                return;
+            }
+        }
+        if (GpuUInt128.Pow2Mod(exponent, q) != GpuUInt128.One)
+        {
+            orders[index] = 0UL;
+            return;
+        }
+
+        orders[index] = exponent;
+    }
+
+    // Device-friendly small cycle length for 2 mod divisor
+    private static ulong CalculateCycleLengthSmall(ulong divisor)
+    {
+        if ((divisor & (divisor - 1UL)) == 0UL)
+            return 1UL;
+
+        ulong order = 1UL, pow = 2UL;
+        while (pow != 1UL)
+        {
+            pow <<= 1;
+            if (pow > divisor)
+                pow -= divisor;
+
+            order++;
+        }
+
+        return order;
+    }
+
+    private static void IncrementalOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong divMul,
+		ResidueAutomatonArgs ra, ArrayView<int> found, ArrayView1D<uint, Stride1D.Dense> smallCycles)
+    {
+        ulong idx = (ulong)index.X;
+		ulong r10 = ra.Q0M10 + (ra.Step10 * idx) % 10UL; r10 -= (r10 >= 10UL) ? 10UL : 0UL;
+		bool shouldCheck = r10 != 5UL;
+		if (shouldCheck)
+		{
+			ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL;
+			if (r8 != 1UL && r8 != 7UL)
+			{
+				shouldCheck = false;
+			}
+			else
+			{
+				ulong r3 = ra.Q0M3 + (ra.Step3 * (idx % 3UL)); r3 -= (r3 >= 3UL) ? 3UL : 0UL;
+				if (r3 == 0UL)
+				{
+					shouldCheck = false;
+				}
+				else
+				{
+					ulong r5 = ra.Q0M5 + (ra.Step5 * (idx % 5UL)); if (r5 >= 5UL) r5 -= 5UL;
+					if (r5 == 0UL)
+					{
+						shouldCheck = false;
+					}
+				}
+			}
+		}
+		if (!shouldCheck)
+		{
+			return;
+		}
+
+        GpuUInt128 k = kStart + (GpuUInt128)idx;
+        GpuUInt128 q = twoP.Mul64(k) + GpuUInt128.One;
+
+        // Small-cycles in-kernel early rejection from device table
+        if (q.High == 0UL && q.Low <= (ulong)(smallCycles.Length - 1))
+        {
+            uint cyc = smallCycles[(int)q.Low];
+            if (cyc != 0U && (exponent % (ulong)cyc) == 0UL)
+            {
+                return;
+            }
+        }
+        GpuUInt128 phi = q - GpuUInt128.One;
+        if (phi.High != 0UL)
+        {
+            return;
+        }
+
+        ulong phi64 = phi.Low;
+        if (GpuUInt128.Pow2Mod(phi64, q) != GpuUInt128.One)
+        {
+            return;
+        }
+
+        GpuUInt128 halfPow = GpuUInt128.Pow2Mod(phi64 >> 1, q) - GpuUInt128.One;
+        if (GpuUInt128.BinaryGcd(halfPow, q) != GpuUInt128.One)
+        {
+            return;
+        }
+
+        ulong div = phi64.FastDiv64(exponent, divMul);
+        GpuUInt128 divPow = GpuUInt128.Pow2Mod(div, q) - GpuUInt128.One;
+        if (GpuUInt128.BinaryGcd(divPow, q) != GpuUInt128.One)
+        {
+            return;
+        }
+
+        Atomic.Or(ref found[0], 1);
+    }
+
+    private static void Pow2ModOrderKernelScan(Index1D index, ulong exponent, GpuUInt128 twoP, GpuUInt128 kStart, byte lastIsSeven, ulong _,
+		ResidueAutomatonArgs ra, ArrayView<int> found, ArrayView1D<uint, Stride1D.Dense> smallCycles)
+    {
+        ulong idx = (ulong)index.X;
+		ulong r10 = ra.Q0M10 + (ra.Step10 * idx) % 10UL; r10 -= (r10 >= 10UL) ? 10UL : 0UL;
+		bool shouldCheck = r10 != 5UL;
+		if (shouldCheck)
+		{
+			ulong r8 = (ra.Q0M8 + ((ra.Step8 * idx) & 7UL)) & 7UL;
+			if (r8 != 1UL && r8 != 7UL)
+			{
+				shouldCheck = false;
+			}
+			else
+			{
+				ulong r3 = ra.Q0M3 + (ra.Step3 * (idx % 3UL)); r3 -= (r3 >= 3UL) ? 3UL : 0UL;
+				if (r3 == 0UL)
+				{
+					shouldCheck = false;
+				}
+				else
+				{
+					ulong r5 = ra.Q0M5 + (ra.Step5 * (idx % 5UL)); if (r5 >= 5UL) r5 -= 5UL;
+					if (r5 == 0UL)
+					{
+						shouldCheck = false;
+					}
+				}
+			}
+		}
+		if (!shouldCheck)
+		{
+			return;
+		}
+
+        GpuUInt128 k = kStart + (GpuUInt128)idx;
+        GpuUInt128 q = twoP.Mul64(k) + GpuUInt128.One;
+        // Small-cycles in-kernel early rejection from device table
+        if (q.High == 0UL && q.Low <= (ulong)(smallCycles.Length - 1))
+        {
+            uint cyc = smallCycles[(int)q.Low];
+            if (cyc != 0U && (exponent % (ulong)cyc) == 0UL)
+            {
+                return;
+            }
+        }
+        if (GpuUInt128.Pow2Mod(exponent, q) != GpuUInt128.One)
+        {
+            return;
+        }
+
+        Atomic.Or(ref found[0], 1);
+    }
+
+	public void Dispose()
+	{
+		// Do not change this code. Put cleanup code in 'Dispose(bool disposing)' method
+		Dispose(disposing: true);
+		// Enabled this when the struct is converted to a class
+		// GC.SuppressFinalize(this);
+	}
+}
+
+public class GpuKernelPool
+{
+    private static readonly ConcurrentDictionary<Accelerator, KernelContainer> KernelCache = new();
+
+    private static KernelContainer GetKernels(Accelerator accelerator)
+    {
+        return KernelCache.GetOrAdd(accelerator, _ => new KernelContainer());
+    }
+
+    // Ensures the small cycles table is uploaded to the device for the given accelerator.
+    // Returns the ArrayView to pass into kernels (when kernels are extended to accept it).
+    public static ArrayView1D<uint, Stride1D.Dense> EnsureSmallCyclesOnDevice(Accelerator accelerator)
+    {
+        var kernels = GetKernels(accelerator);
+        if (kernels.SmallCycles is { } buffer)
+        {
+            return buffer.View;
+        }
+
+        // Ensure single upload per accelerator even if multiple threads race here.
+        lock (kernels)
+        {
+            if (kernels.SmallCycles is { } existing)
+            {
+                return existing.View;
+            }
+
+            var host = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
+            var device = accelerator.Allocate1D<uint>(host.Length);
+            device.View.CopyFromCPU(host);
+            kernels.SmallCycles = device;
+            return device.View;
+        }
+    }
+
+	public static GpuKernelLease GetKernel(bool useGpuOrder)
+	{
+		var limiter = GpuPrimeWorkLimiter.Acquire();
+		var gpu = RentPreferred(preferCpu: !useGpuOrder);
+		var accelerator = gpu.Accelerator;
+		var kernels = GetKernels(accelerator);
+		return new(limiter, gpu, kernels);
+	}
+
+	/// <summary>
+	/// Runs a GPU action with an acquired accelerator and stream.
+	/// </summary>
+	/// <param name="action">Action to run with (Accelerator, Stream).</param>
+    public static void Run(Action<Accelerator, AcceleratorStream> action)
+    {
+        using var lease = GetKernel(useGpuOrder: true);
+        var accelerator = lease.Accelerator;
+        using var stream = accelerator.CreateStream();
+        action(accelerator, stream);
+    }
+}

PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs

@@ -0,0 +1,57 @@
+namespace PerfectNumbers.Core.Gpu;
+
+public static class GpuPrimeWorkLimiter
+{
+    private static SemaphoreSlim _semaphore = new(1, int.MaxValue);
+    private static readonly object _lock = new();
+    private static int _currentLimit = 1;
+
+    public static IDisposable Acquire()
+    {
+        var sem = _semaphore;
+        sem.Wait();
+        return new Releaser(sem);
+    }
+
+    public static void SetLimit(int value)
+    {
+        if (value < 1)
+        {
+            value = 1;
+        }
+
+        lock (_lock)
+        {
+            if (value == _currentLimit)
+            {
+                return;
+            }
+
+            _semaphore = new SemaphoreSlim(value, value);
+            _currentLimit = value;
+        }
+    }
+
+    private sealed class Releaser : IDisposable
+    {
+        private readonly SemaphoreSlim _sem;
+        private bool _disposed;
+
+        public Releaser(SemaphoreSlim sem)
+        {
+            _sem = sem;
+        }
+
+        public void Dispose()
+        {
+            if (_disposed)
+            {
+                return;
+            }
+
+            _sem.Release();
+            _disposed = true;
+        }
+    }
+}
+

PerfectNumbers.Core/Gpu/GpuUInt128.cs

@@ -0,0 +1,991 @@
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public readonly struct GpuUInt128 : IComparable<GpuUInt128>, IEquatable<GpuUInt128>
+{
+    public readonly ulong High;
+
+    public readonly ulong Low;
+
+    public static readonly GpuUInt128 Zero = new(0UL, 0UL);
+    public static readonly GpuUInt128 One = new(0UL, 1UL);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128(ulong high, ulong low)
+    {
+        High = high;
+        Low = low;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128(UInt128 value)
+    {
+        High = (ulong)(value >> 64);
+        Low = (ulong)value;
+    }
+
+    public readonly bool IsZero
+    {
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        get => High == 0UL && Low == 0UL;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static implicit operator GpuUInt128(UInt128 value) => new(value);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static implicit operator GpuUInt128(ulong value) => new(0UL, value);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static implicit operator UInt128(GpuUInt128 value) =>
+        ((UInt128)value.High << 64) | value.Low;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 operator +(GpuUInt128 left, GpuUInt128 right)
+    {
+        var low = left.Low + right.Low;
+        var high = left.High + right.High + (low < left.Low ? 1UL : 0UL);
+        return new(high, low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 operator -(GpuUInt128 left, GpuUInt128 right)
+    {
+        var borrow = left.Low < right.Low ? 1UL : 0UL;
+        var low = left.Low - right.Low;
+        var high = left.High - right.High - borrow;
+        return new(high, low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 operator ^(GpuUInt128 left, GpuUInt128 right) =>
+        new(left.High ^ right.High, left.Low ^ right.Low);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 operator *(GpuUInt128 left, GpuUInt128 right)
+    {
+        Multiply(left, right, out var high, out var low);
+        return new(high, low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 operator <<(GpuUInt128 value, int shift)
+    {
+        shift &= 127;
+        if (shift == 0)
+        {
+            return value;
+        }
+
+        if (shift >= 64)
+        {
+            ulong high = value.Low << (shift - 64);
+            return new(high, 0UL);
+        }
+
+        ulong highPart = (value.High << shift) | (value.Low >> (64 - shift));
+        ulong lowPart = value.Low << shift;
+        return new(highPart, lowPart);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 operator >>(GpuUInt128 value, int shift)
+    {
+        shift &= 127;
+        if (shift == 0)
+        {
+            return value;
+        }
+
+        if (shift >= 64)
+        {
+            ulong low = value.High >> (shift - 64);
+            return new(0UL, low);
+        }
+
+        ulong highPart = value.High >> shift;
+        ulong lowPart = (value.Low >> shift) | (value.High << (64 - shift));
+        return new(highPart, lowPart);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool operator ==(GpuUInt128 left, GpuUInt128 right) => left.Equals(right);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool operator !=(GpuUInt128 left, GpuUInt128 right) => !left.Equals(right);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool operator <(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) < 0;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool operator >(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) > 0;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool operator <=(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) <= 0;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool operator >=(GpuUInt128 left, GpuUInt128 right) => left.CompareTo(right) >= 0;
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static GpuUInt128 operator ++(GpuUInt128 value) => value.Add(1UL);
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public readonly int CompareTo(GpuUInt128 other)
+    {
+        if (High < other.High)
+        {
+            return -1;
+        }
+
+        if (High > other.High)
+        {
+            return 1;
+        }
+
+        if (Low < other.Low)
+        {
+            return -1;
+        }
+
+        if (Low > other.Low)
+        {
+            return 1;
+        }
+
+        return 0;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public readonly bool Equals(GpuUInt128 other) => High == other.High && Low == other.Low;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public override readonly bool Equals(object? obj) => obj is GpuUInt128 value && Equals(value);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public override readonly int GetHashCode() => HashCode.Combine(High, Low);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 Add(GpuUInt128 other)
+    {
+        ulong low = Low + other.Low;
+        return new(High + other.High + (low < Low ? 1UL : 0UL), low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 Add(ulong value)
+    {
+        var low = Low + value;
+        return new(High + (low < Low ? 1UL : 0UL), low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 AddMod(GpuUInt128 value, GpuUInt128 modulus)
+    {
+        if (modulus.High == 0UL)
+        {
+            return AddMod(value, modulus.Low);
+        }
+
+        value = Add(value);
+        if (value.CompareTo(modulus) >= 0)
+        {
+            value = value.Sub(modulus);
+        }
+
+        return value;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 AddMod(ulong value, GpuUInt128 modulus)
+    {
+        if (modulus.High == 0UL)
+        {
+            return AddMod(value, modulus.Low);
+        }
+
+        var result = Add(value);
+        if (result.CompareTo(modulus) >= 0)
+        {
+            result = result.Sub(modulus);
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 AddMod(GpuUInt128 value, ulong modulus)
+    {
+        ulong a = Low % modulus;
+        ulong b = value.Low % modulus;
+        ulong sum = a + b;
+        if (sum >= modulus || sum < a)
+        {
+            sum -= modulus;
+        }
+
+        return new(0UL, sum);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 AddMod(ulong value, ulong modulus)
+    {
+        ulong a = Low % modulus;
+        ulong b = value % modulus;
+        ulong sum = a + b;
+        if (sum >= modulus || sum < a)
+        {
+            sum -= modulus;
+        }
+
+        return new(0UL, sum);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 Sub(GpuUInt128 other)
+    {
+        ulong borrow = Low < other.Low ? 1UL : 0UL;
+        return new(High - other.High + borrow, Low - other.Low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static ulong MulHigh(ulong x, ulong y)
+    {
+        ulong xLow = (uint)x;
+        ulong xHigh = x >> 32;
+        ulong yLow = (uint)y;
+        ulong yHigh = y >> 32;
+
+        ulong w1 = xLow * yHigh;
+        ulong w2 = xHigh * yLow;
+        ulong carry =
+            (((xLow * yLow) >> 32) +
+            (uint)w1 +
+            (uint)w2) >> 32;
+        return (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32) + carry;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public readonly GpuUInt128 Mul64(GpuUInt128 other)
+    {
+        // Multiply this.Low (assumed 64-bit value) by full 128-bit other
+        ulong a = Low;
+        ulong low = a * other.Low;
+        ulong high = a * other.High + MulHigh(a, other.Low);
+        return new(high, low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 AddMod(GpuUInt128 a, GpuUInt128 b, GpuUInt128 modulus)
+    {
+        a = a.Add(b);
+        if (a >= modulus)
+        {
+            a = a.Sub(modulus);
+        }
+
+        return a;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 MulMod(GpuUInt128 a, GpuUInt128 b, GpuUInt128 modulus)
+    {
+        GpuUInt128 result = new(0UL, 0UL);
+        while (!b.IsZero)
+        {
+            if ((b.Low & 1UL) != 0UL)
+            {
+                result = AddMod(result, a, modulus);
+            }
+
+            // a = (a << 1) % modulus
+            a <<= 1;
+            if (a >= modulus)
+            {
+                a = a.Sub(modulus);
+            }
+
+            // b >>= 1
+            if (b.High == 0UL)
+            {
+                b = new GpuUInt128(0UL, b.Low >> 1);
+            }
+            else
+            {
+                ulong low = (b.Low >> 1) | (b.High << 63);
+                ulong high = b.High >> 1;
+                b = new GpuUInt128(high, low);
+            }
+        }
+
+        return result;
+    }
+
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 Pow2Mod(ulong exponent, GpuUInt128 modulus)
+    {
+        if (modulus.IsZero || modulus == One)
+        {
+            return new(0UL, 0UL);
+        }
+
+        GpuUInt128 result = new(0UL, 1UL);
+		GpuUInt128 baseVal = new(0UL, 2UL);
+
+        ulong e = exponent;
+        while (e != 0UL)
+        {
+            if ((e & 1UL) != 0UL)
+            {
+                result = MulMod(result, baseVal, modulus);
+            }
+
+            baseVal = MulMod(baseVal, baseVal, modulus);
+            e >>= 1;
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static int TrailingZeroCount(GpuUInt128 value)
+    {
+        if (value.Low != 0UL)
+        {
+            return BitOperations.TrailingZeroCount(value.Low);
+        }
+
+        return 64 + BitOperations.TrailingZeroCount(value.High);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static GpuUInt128 BinaryGcd(GpuUInt128 u, GpuUInt128 v)
+    {
+        if (u.IsZero)
+        {
+            return v;
+        }
+
+        if (v.IsZero)
+        {
+            return u;
+        }
+
+        int shift = TrailingZeroCount(new GpuUInt128(u.High | v.High, u.Low | v.Low));
+        // remove factors of 2
+        int zu = TrailingZeroCount(u);
+        u >>= zu;
+
+        do
+        {
+            int zv = TrailingZeroCount(v);
+            v >>= zv;
+            if (u > v)
+            {
+                (u, v) = (v, u);
+            }
+
+            v -= u;
+        }
+        while (!v.IsZero);
+
+        return u << shift;
+    }
+
+
+    public GpuUInt128 Sub(ulong value)
+    {
+        ulong borrow = Low < value ? 1UL : 0UL;
+        return new(High - borrow, Low - value);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 SubMod(GpuUInt128 value, GpuUInt128 modulus)
+    {
+        if (CompareTo(value) >= 0)
+        {
+            return Sub(value);
+        }
+        else
+        {
+            return Add(modulus).Sub(value);
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 SubMod(ulong value, GpuUInt128 modulus)
+    {
+		ulong low, high;
+		if (High == 0UL && Low < value)
+		{
+			low = Low + modulus.Low;
+			ulong carry = low < Low ? 1UL : 0UL;
+			high = High + modulus.High + carry;
+		}
+		else
+		{
+			high = High;
+			low = Low;
+		}
+
+        ulong borrow = low < value ? 1UL : 0UL;
+        return new (high - borrow, low - value);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 Mul(GpuUInt128 other)
+    {
+        Multiply(this, other, out var high, out var low);
+        return new(high, low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 Mul(ulong value)
+    {
+        var (highPart, lowPart) = Mul64(Low, value);
+        return new(highPart + High * value, lowPart);
+    }
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public GpuUInt128 Xor(GpuUInt128 other) => new(High ^ other.High, Low ^ other.Low);
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public GpuUInt128 Xor(ulong value) => new(High, Low ^ value);
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public GpuUInt128 ShiftLeft(int shift) => this << shift;
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public GpuUInt128 ShiftRight(int shift) => this >> shift;
+
+	public GpuUInt128 MulMod(GpuUInt128 other, GpuUInt128 modulus)
+    {
+        if (modulus.High == 0UL)
+        {
+            return MulMod(other, modulus.Low);
+        }
+
+        MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0);
+
+        GpuUInt128 remainder = new(0UL, 0UL);
+        int bit;
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p3 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p2 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p1 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p0 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+        return remainder;
+        // TODO(MOD-OPT): Replace this bitwise long-division style reduction
+        // with a faster algorithm suitable for GPU kernels without '%' support.
+        // Options:
+        // - Montgomery reduction for 128-bit moduli (R=2^128), requires pre-
+        //   computed modulus-dependent constants (n' and R^2 mod n). This is
+        //   ideal for NTT primes (odd modulus). Cache per-modulus constants.
+        // - Barrett reduction with a 256/128 quotient approximation using only
+        //   multiplies and shifts. Cache mu = floor(2^k / n) for k=256.
+        // Implement a fast path for common 64-bit NTT moduli and a separate
+        // path for 128-bit moduli. This will significantly reduce the cost of
+        // each butterfly in NTT and LL steps.
+        // TODO(MOD-OPT): Plumb constants through caches in NttGpuMath.SquareCacheEntry
+        // (e.g., Montgomery n', R2) and provide device-friendly accessors.
+    }
+
+    /// <summary>
+    /// Modular multiplication using <see cref="BigInteger"/> reduction. This
+    /// method is intended for validation only and should not be used inside
+    /// GPU kernels.
+    /// </summary>
+    internal GpuUInt128 MulModBigInteger(GpuUInt128 other, GpuUInt128 modulus)
+    {
+        var left = (BigInteger)(UInt128)this;
+        var right = (BigInteger)(UInt128)other;
+        var mod = (BigInteger)(UInt128)modulus;
+        var reduced = (UInt128)((left * right) % mod);
+        return new((ulong)(reduced >> 64), (ulong)reduced);
+    }
+
+    /// <summary>
+    /// Experimental limb-based reduction. The current implementation performs
+    /// repeated subtractions and becomes extremely slow for large remainders.
+    /// Kept for future optimization work.
+    /// </summary>
+    internal GpuUInt128 MulModByLimb(GpuUInt128 other, GpuUInt128 modulus)
+    {
+        MultiplyFull(this, other, out var p3, out var p2, out var p1, out var p0);
+
+        GpuUInt128 remainder = new(p3, p2);
+        while (remainder.CompareTo(modulus) >= 0)
+        {
+            remainder = remainder.Sub(modulus);
+        }
+
+        ulong limb = p1;
+        for (int i = 0; i < 2; i++)
+        {
+            remainder <<= 64;
+            remainder = new(remainder.High, limb);
+            while (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+
+            limb = p0;
+        }
+
+        return remainder;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 MulMod(ulong value, GpuUInt128 modulus)
+    {
+        if (modulus.High == 0UL)
+        {
+            return MulMod(value, modulus.Low);
+        }
+
+        GpuUInt128 other = new(0UL, value);
+        return MulMod(other, modulus);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 MulMod(GpuUInt128 other, ulong modulus)
+    {
+        // Fast-path: use single 64-bit multiply when it cannot overflow
+        // (enabled by ILGPU support for 64-bit % in kernels). Otherwise fall
+        // back to the generic branchless double-and-add implementation.
+        ulong a = Low % modulus;
+        ulong b = other.Low % modulus;
+        ulong result;
+        if (a == 0UL || b == 0UL)
+        {
+            result = 0UL;
+        }
+        else if (b <= ulong.MaxValue / a)
+        {
+            // Safe to multiply without overflow
+            result = (a * b) % modulus;
+        }
+        else
+        {
+            result = MulMod64(a, b, modulus);
+        }
+
+        return new(0UL, result);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 MulMod(ulong value, ulong modulus)
+    {
+        // Same fast-path as above for 64-bit operands
+        ulong a = Low % modulus;
+        ulong b = value % modulus;
+        ulong result;
+        if (a == 0UL || b == 0UL)
+        {
+            result = 0UL;
+        }
+        else if (b <= ulong.MaxValue / a)
+        {
+            result = (a * b) % modulus;
+        }
+        else
+        {
+            result = MulMod64(a, b, modulus);
+        }
+
+        return new(0UL, result);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 MulModMontgomery64(GpuUInt128 other, ulong modulus, ulong nPrime, ulong r2)
+    {
+        // Montgomery-based 64-bit modular multiplication without division.
+        // Computes (this * other) % modulus using:
+        // aR = MontMul(a, R^2), bR = MontMul(b, R^2), cR = MontMul(aR, bR), c = MontMul(cR, 1)
+        // where R = 2^64, nPrime = -mod^{-1} mod 2^64, r2 = R^2 mod mod.
+        ulong a = Low % modulus;
+        ulong b = other.Low % modulus;
+        if (a == 0UL || b == 0UL)
+        {
+			return Zero;
+        }
+
+        ulong aR = MontMul64(a, r2, modulus, nPrime);
+        ulong bR = MontMul64(b, r2, modulus, nPrime);
+        ulong cR = MontMul64(aR, bR, modulus, nPrime);
+        ulong c = MontMul64(cR, 1UL, modulus, nPrime);
+        return new(0UL, c);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 MulModMontgomery64(ulong other, ulong modulus, ulong nPrime, ulong r2)
+    {
+        ulong a = Low % modulus;
+        ulong b = other % modulus;
+        if (a == 0UL || b == 0UL)
+        {
+			return Zero; 
+        }
+
+        ulong aR = MontMul64(a, r2, modulus, nPrime);
+        ulong bR = MontMul64(b, r2, modulus, nPrime);
+        ulong cR = MontMul64(aR, bR, modulus, nPrime);
+        ulong c = MontMul64(cR, 1UL, modulus, nPrime);
+        return new(0UL, c);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 SquareMod(GpuUInt128 modulus)
+    {
+        if (modulus.High == 0UL)
+        {
+            // Use fast 64-bit path.
+            return MulMod(Low, modulus.Low);
+        }
+
+        // Specialized 128-bit squaring with bitwise reduction (faster than generic multiply).
+        SquareFull(this, out var p3, out var p2, out var p1, out var p0);
+
+        GpuUInt128 remainder = new(0UL, 0UL);
+        int bit;
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p3 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p2 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p1 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+        for (bit = 63; bit >= 0; bit--)
+        {
+            remainder <<= 1;
+            remainder = new(remainder.High, remainder.Low | ((p0 >> bit) & 1UL));
+            if (remainder.CompareTo(modulus) >= 0)
+            {
+                remainder = remainder.Sub(modulus);
+            }
+        }
+
+		return remainder;
+    }
+
+    public GpuUInt128 ModPow(GpuUInt128 exponent, GpuUInt128 modulus)
+    {
+        GpuUInt128 result = new(0UL, 1UL);
+        GpuUInt128 baseValue = this;
+
+        while (!exponent.IsZero)
+        {
+            if ((exponent.Low & 1UL) != 0UL)
+            {
+                result = result.MulMod(baseValue, modulus);
+            }
+
+            exponent >>= 1;
+            baseValue = baseValue.MulMod(baseValue, modulus);
+        }
+
+		return result;
+    }
+
+    public GpuUInt128 ModPow(ulong exponent, GpuUInt128 modulus)
+    {
+        GpuUInt128 result = new(0UL, 1UL);
+        GpuUInt128 baseValue = this;
+
+        while (exponent != 0UL)
+        {
+            if ((exponent & 1UL) != 0UL)
+            {
+                result = result.MulMod(baseValue, modulus);
+            }
+
+            exponent >>= 1;
+            baseValue = baseValue.MulMod(baseValue, modulus);
+        }
+
+        return result;
+    }
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public GpuUInt128 ModInv(ulong modulus) => ModPow(modulus - 2UL, new GpuUInt128(0UL, modulus));
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public GpuUInt128 ModInv(GpuUInt128 modulus)
+    {
+        if (modulus.High == 0UL)
+        {
+            return ModInv(modulus.Low);
+        }
+
+        var exponent = new GpuUInt128(modulus).Sub(2UL);
+        return ModPow(exponent, modulus);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public override readonly string ToString() =>
+        ((UInt128)this).ToString();
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void MultiplyFull(GpuUInt128 left, GpuUInt128 right, out ulong p3, out ulong p2, out ulong p1, out ulong p0)
+    {
+        var (h0, l0) = Mul64(left.Low, right.Low);
+        var (h1, l1) = Mul64(left.Low, right.High);
+        var (h2, l2) = Mul64(left.High, right.Low);
+        var (h3, l3) = Mul64(left.High, right.High);
+
+        p0 = l0;
+        ulong carry = 0UL;
+        ulong sum = h0;
+        sum += l1;
+        if (sum < l1)
+        {
+            carry++;
+        }
+
+        sum += l2;
+        if (sum < l2)
+        {
+            carry++;
+        }
+
+        p1 = sum;
+        sum = h1;
+        sum += h2;
+        ulong carry2 = sum < h2 ? 1UL : 0UL;
+        sum += l3;
+        if (sum < l3)
+        {
+            carry2++;
+        }
+
+        sum += carry;
+        if (sum < carry)
+        {
+            carry2++;
+        }
+
+        p2 = sum;
+        p3 = h3 + carry2;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void Multiply(GpuUInt128 left, GpuUInt128 right, out ulong high, out ulong low)
+    {
+        var (h0, l0) = Mul64(left.Low, right.Low);
+        var (_, l1) = Mul64(left.Low, right.High);
+        var (_, l2) = Mul64(left.High, right.Low);
+
+        low = l0;
+        high = h0 + l1 + l2;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static (ulong High, ulong Low) Mul64(ulong left, ulong right)
+    {
+        ulong a0 = (uint)left;
+        ulong a1 = left >> 32;
+        ulong b0 = (uint)right;
+        ulong b1 = right >> 32;
+
+        ulong lo = a0 * b0;
+        ulong mid1 = a1 * b0;
+        ulong mid2 = a0 * b1;
+        ulong hi = a1 * b1;
+
+        ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2;
+        ulong low = (lo & 0xFFFFFFFFUL) | (carry << 32);
+        hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32);
+
+        return (hi, low);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void SquareFull(GpuUInt128 value, out ulong p3, out ulong p2, out ulong p1, out ulong p0)
+    {
+        // Compute (H*2^64 + L)^2 = H^2*2^128 + 2*H*L*2^64 + L^2
+        ulong L = value.Low;
+        ulong H = value.High;
+        var (hLL, lLL) = Mul64(L, L);      // L^2 -> (hLL,lLL)
+        var (hHH, lHH) = Mul64(H, H);      // H^2 -> (hHH,lHH)
+        var (hLH, lLH) = Mul64(L, H);      // L*H -> (hLH,lLH)
+
+        // double LH: (hLH,lLH) << 1
+        ulong dLH_low = lLH << 1;
+        ulong carry = (lLH >> 63) & 1UL;
+        ulong dLH_high = (hLH << 1) | carry;
+
+        // Assemble 256-bit result parts
+        p0 = lLL;
+
+        // p1 = hLL + dLH_low (with carry)
+        ulong sum1 = hLL + dLH_low;
+        ulong c1 = sum1 < hLL ? 1UL : 0UL;
+        p1 = sum1;
+
+        // p2 = lHH + dLH_high + c1
+        ulong sum2 = lHH + dLH_high;
+        ulong c2 = sum2 < lHH ? 1UL : 0UL;
+        sum2 += c1;
+        if (sum2 < c1)
+        {
+            c2++;
+        }
+        p2 = sum2;
+
+        // p3 = hHH + c2
+        p3 = hHH + c2;
+    }
+
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong MulMod64(ulong a, ulong b, ulong modulus)
+    {
+        // Note: 64-bit % is supported in ILGPU kernels, but the straightforward
+        // (UInt128)a*b % modulus is not, because 128-bit multiply is not
+        // supported in device code. We therefore implement a branch-free
+        // shift-add reduction here based solely on 64-bit ops.
+        ulong result = 0UL;
+        ulong x = a % modulus;
+        ulong y = b;
+        while (y != 0UL)
+        {
+            if ((y & 1UL) != 0UL)
+            {
+                result += x;
+                if (result >= modulus)
+                {
+                    result -= modulus;
+                }
+            }
+
+            x <<= 1;
+            if (x >= modulus)
+            {
+                x -= modulus;
+            }
+
+            y >>= 1;
+        }
+
+        return result;
+    }
+
+    // Montgomery core for 64-bit operands. Returns a*b*R^{-1} mod modulus.
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong MontMul64(ulong aR, ulong bR, ulong modulus, ulong nPrime)
+    {
+        // t = aR * bR (128-bit)
+        ulong tLow, tHigh;
+        Mul64Parts(aR, bR, out tHigh, out tLow);
+        // m = (tLow * nPrime) mod 2^64 (low 64 bits only)
+        ulong mLow, mHigh;
+        Mul64Parts(tLow, nPrime, out mHigh, out mLow);
+        // u = (t + m * modulus) >> 64
+        ulong mmLow, mmHigh;
+        Mul64Parts(mLow, modulus, out mmHigh, out mmLow);
+        ulong carry = 0UL;
+        ulong low = tLow + mmLow;
+        if (low < tLow)
+        {
+            carry = 1UL;
+        }
+
+        ulong high = tHigh + mmHigh + carry;
+        ulong u = high; // (t + m*n) >> 64
+        if (u >= modulus)
+        {
+            u -= modulus;
+        }
+
+        return u;
+    }
+
+    // 64x64 -> 128 multiply into (high, low)
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low)
+    {
+        ulong a0 = (uint)a;
+        ulong a1 = a >> 32;
+        ulong b0 = (uint)b;
+        ulong b1 = b >> 32;
+
+        ulong lo = a0 * b0;
+        ulong mid1 = a1 * b0;
+        ulong mid2 = a0 * b1;
+        ulong hi = a1 * b1;
+
+        ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2;
+        low = (lo & 0xFFFFFFFFUL) | (carry << 32);
+        hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32);
+        high = hi;
+    }
+}
+
+    // TODO(MOD-OPT): Montgomery/Barrett integration plan
+    // - Introduce caches of modulus-dependent constants:
+    //   * Montgomery: n' (-(n^{-1}) mod 2^64 or 2^128), R2 = (R^2 mod n)
+    //   * Barrett: mu = floor(2^k / n) for k ∈ {128, 192, 256}
+    // - Add fast-path for 64-bit NTT primes (modulus.High == 0UL) using pure 64-bit ops.
+    // - Expose helpers to retrieve/calc constants once per modulus and reuse in kernels.
+    // - Wire these into MulMod and SquareMod hot paths under feature toggles.
+    // - Ensure ILGPU compatibility (no BigInteger, no % inside kernels).

PerfectNumbers.Core/Gpu/GpuWorkLimiter.cs

@@ -0,0 +1,60 @@
+namespace PerfectNumbers.Core.Gpu;
+
+public static class GpuWorkLimiter
+{
+    private static SemaphoreSlim _semaphore = new(1, int.MaxValue);
+    private static readonly object _lock = new();
+    private static int _currentLimit = 1;
+
+    public static IDisposable Acquire()
+    {
+        var sem = _semaphore;
+        sem.Wait();
+        return new Releaser(sem);
+    }
+
+    public static void SetLimit(int value)
+    {
+        if (value < 1)
+        {
+            value = 1;
+        }
+
+        lock (_lock)
+        {
+            if (value == _currentLimit)
+            {
+                return;
+            }
+
+            // Swap to a new semaphore for the new limit. We intentionally
+            // do not dispose the old semaphore to avoid releasing handles
+            // that may still be in use by existing Releaser instances.
+            _semaphore = new SemaphoreSlim(value, value);
+            _currentLimit = value;
+        }
+    }
+
+    private sealed class Releaser : IDisposable
+    {
+        private readonly SemaphoreSlim _sem;
+        private bool _disposed;
+
+        public Releaser(SemaphoreSlim sem)
+        {
+            _sem = sem;
+        }
+
+        public void Dispose()
+        {
+            if (_disposed)
+            {
+                return;
+            }
+
+            _sem.Release();
+            _disposed = true;
+        }
+    }
+}
+

PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs

@@ -0,0 +1,77 @@
+using System.Buffers;
+using System.Runtime.InteropServices;
+using ILGPU.Runtime;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool useGpuOrder)
+{
+	private readonly GpuKernelType _kernelType = kernelType;
+	private readonly bool _useGpuOrder = useGpuOrder;
+
+	public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
+	{
+		using var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
+		var accelerator = gpuLease.Accelerator;
+		int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy
+		UInt128 kStart = 1UL;
+		ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
+		byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
+
+        var kernel = (_kernelType == GpuKernelType.Pow2Mod)
+            ? gpuLease.Pow2ModKernel
+            : gpuLease.IncrementalKernel;
+
+		using var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
+		// Avoid giant stack allocations that can trigger StackOverflow when batchSize is large.
+		// Rent a reusable array from the shared pool instead.
+		ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
+		UInt128 remaining;
+		int currentSize;
+		int i;
+		UInt128 q = UInt128.Zero;
+		try
+		{
+			while (kStart <= maxK && Volatile.Read(ref isPrime))
+			{
+				remaining = maxK - kStart + UInt128.One;
+				currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
+				Span<ulong> orders = orderArray.AsSpan(0, currentSize);
+				// Precompute residue automaton bases for this batch
+				UInt128 q0 = twoP * kStart + UInt128.One;
+				ulong q0m10 = q0.Mod10();
+				ulong q0m8 = q0.Mod8();
+				ulong q0m3 = q0.Mod3();
+				ulong q0m5 = q0.Mod5();
+				ulong step10 = ((exponent % 10UL) << 1) % 10UL;
+				ulong step8 = ((exponent & 7UL) << 1) & 7UL;
+				ulong step3 = ((exponent % 3UL) << 1) % 3UL;
+				ulong step5 = ((exponent % 5UL) << 1) % 5UL;
+
+                var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
+                kernel(currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, divMul,
+                    q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView);
+
+				accelerator.Synchronize();
+				orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
+				q = twoP * kStart + UInt128.One;
+				for (i = 0; i < currentSize && Volatile.Read(ref isPrime); i++, q += twoP)
+				{
+					if (orders[i] != 0UL)
+					{
+						if (q.IsPrimeCandidate())
+						{
+							Volatile.Write(ref isPrime, false);
+						}
+					}
+				}
+
+				kStart += (UInt128)currentSize;
+			}
+		}
+		finally
+		{
+			ArrayPool<ulong>.Shared.Return(orderArray);
+		}
+	}
+}

PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs

@@ -0,0 +1,864 @@
+using System.Collections.Concurrent;
+using ILGPU;
+using ILGPU.Runtime;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public class MersenneNumberLucasLehmerGpuTester
+{
+    private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, GpuUInt128, ArrayView<GpuUInt128>>> KernelCache = new();
+    private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<GpuUInt128>>> AddSmallKernelCache = new();
+    private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<GpuUInt128>>> SubSmallKernelCache = new();
+    private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ulong, ArrayView<GpuUInt128>>> ReduceKernelCache = new();
+    private readonly ConcurrentDictionary<Accelerator, Action<Index1D, ArrayView<GpuUInt128>, ArrayView<byte>>> IsZeroKernelCache = new();
+
+    private Action<Index1D, ulong, GpuUInt128, ArrayView<GpuUInt128>> GetKernel(Accelerator accelerator)
+    {
+        return KernelCache.GetOrAdd(accelerator, accel => accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, ArrayView<GpuUInt128>>(Kernel));
+    }
+
+    private Action<Index1D, ulong, ArrayView<GpuUInt128>> GetAddSmallKernel(Accelerator accelerator) =>
+        AddSmallKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<GpuUInt128>>(AddSmallKernel));
+
+    private Action<Index1D, ulong, ArrayView<GpuUInt128>> GetSubSmallKernel(Accelerator accelerator) =>
+        SubSmallKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<GpuUInt128>>(SubtractSmallKernel));
+
+    private Action<Index1D, ulong, ArrayView<GpuUInt128>> GetReduceKernel(Accelerator accelerator) =>
+        ReduceKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView<GpuUInt128>>(ReduceModMersenneKernel));
+
+    private Action<Index1D, ArrayView<GpuUInt128>, ArrayView<byte>> GetIsZeroKernel(Accelerator accelerator) =>
+        IsZeroKernelCache.GetOrAdd(accelerator, acc => acc.LoadAutoGroupedStreamKernel<Index1D, ArrayView<GpuUInt128>, ArrayView<byte>>(IsZeroKernel));
+
+    // Configurable LL slice size to keep kernels short. Default 32.
+    public int SliceSize = 32;
+
+    public bool IsMersennePrime(ulong exponent)
+    {
+        // Default to global kernel preference for backward compatibility.
+        bool runOnGpu = !GpuContextPool.ForceCpu;
+        return IsPrime(exponent, runOnGpu);
+    }
+
+    public bool IsPrime(ulong exponent, bool runOnGpu)
+    {
+        // Early rejections aligned with incremental/order sieves, but safe for small p:
+        // - If 3 | p and p != 3, then 7 | M_p -> composite.
+        // - If p ≡ 1 (mod 4) and p shares a factor with (p-1), reject fast.
+        if ((exponent % 3UL) == 0UL && exponent != 3UL)
+        {
+            return false;
+        }
+
+        if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
+        {
+            return false;
+        }
+
+        using var limiter = GpuPrimeWorkLimiter.Acquire();
+        if (exponent >= 128UL)
+        {
+            if (!TryGetNttParameters(exponent, out var nttMod, out var primitiveRoot))
+            {
+                throw new NotSupportedException("NTT parameters for the given exponent are not available.");
+            }
+
+            return IsMersennePrimeNtt(exponent, nttMod, primitiveRoot, runOnGpu);
+        }
+
+        using var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu);
+        var accelerator = gpu.Accelerator;
+        var kernel = GetKernel(accelerator);
+        var modulus = new GpuUInt128(((UInt128)1 << (int)exponent) - 1UL);
+        using var buffer = accelerator.Allocate1D<GpuUInt128>(1);
+        kernel(1, exponent, modulus, buffer.View);
+        accelerator.Synchronize();
+        return buffer.GetAsArray1D()[0].IsZero;
+    }
+
+    private bool IsMersennePrimeNtt(ulong exponent, GpuUInt128 nttMod, GpuUInt128 primitiveRoot, bool runOnGpu)
+    {
+        using var gpu = GpuContextPool.RentPreferred(preferCpu: !runOnGpu);
+        var accelerator = gpu.Accelerator;
+        var addKernel = GetAddSmallKernel(accelerator);
+        var subKernel = GetSubSmallKernel(accelerator);
+        var reduceKernel = GetReduceKernel(accelerator);
+        var zeroKernel = GetIsZeroKernel(accelerator);
+
+        int limbCount = (int)((exponent + 127UL) / 128UL);
+        int target = limbCount * 2;
+        int transformLength = 1;
+        while (transformLength < target)
+        {
+            transformLength <<= 1;
+        }
+        using var stateBuffer = accelerator.Allocate1D<GpuUInt128>(transformLength);
+        stateBuffer.MemSetToZero();
+        addKernel(1, 4UL, stateBuffer.View);
+
+        // TODO(LL-SLICE): Slice Lucas–Lehmer iterations into short batches to
+        // avoid long-running kernels and TDR. Example: process 8–64 iterations
+        // per slice, with synchronization between slices. Combine this with
+        // stage-wise NTT to ensure each kernel stays < ~0.5–1.0s.
+        ulong i = 0UL;
+        ulong limit = exponent - 2UL;
+        int slice = Math.Max(1, SliceSize); // process LL iterations in small slices
+        while (i < limit)
+        {
+            int iter = (int)Math.Min((ulong)slice, limit - i);
+            for (int s = 0; s < iter; s++)
+            {
+                // Square in NTT domain with staged kernels (short runtime per kernel).
+                NttGpuMath.SquareDevice(accelerator, stateBuffer.View, nttMod, primitiveRoot);
+                subKernel(1, 2UL, stateBuffer.View);
+                reduceKernel(1, exponent, stateBuffer.View);
+            }
+
+            // Synchronize between slices to yield to scheduler and avoid TDR.
+            accelerator.Synchronize();
+            i += (ulong)iter;
+        }
+
+        reduceKernel(1, exponent, stateBuffer.View);
+
+        using var resultBuffer = accelerator.Allocate1D<byte>(1);
+        zeroKernel(1, stateBuffer.View, resultBuffer.View);
+        accelerator.Synchronize();
+        byte[] result = new byte[1];
+        resultBuffer.View.CopyToCPU(result);
+        return result[0] != 0;
+    }
+
+    private static readonly ConcurrentDictionary<int, (GpuUInt128 Modulus, GpuUInt128 PrimitiveRoot)> NttParameterCache = new()
+    {
+        // 2^12 transform length using a 64-bit NTT-friendly prime
+        [4096] = (new GpuUInt128(0UL, 18446744069414584321UL), new GpuUInt128(0UL, 7UL)),
+        // 2^22 transform length precomputed for p ~= 138,000,000
+        [4194304] = (new GpuUInt128(0UL, 104857601UL), new GpuUInt128(0UL, 39193363UL))
+    };
+
+    private readonly object ParameterFileLock = new();
+    private readonly string ParameterFilePath = Path.Combine(AppContext.BaseDirectory, "ntt-params.txt");
+
+    public void WarmUpNttParameters(ulong exponent)
+    {
+        _ = TryGetNttParameters(exponent, out _, out _);
+    }
+
+    private bool TryGetNttParameters(ulong exponent, out GpuUInt128 modulus, out GpuUInt128 primitiveRoot)
+    {
+        int limbCount = (int)((exponent + 127UL) / 128UL);
+        int length = 1;
+        int target = limbCount * 2;
+        while (length < target)
+        {
+            length <<= 1;
+        }
+
+        if (NttParameterCache.TryGetValue(length, out var cached))
+        {
+            modulus = cached.Modulus;
+            primitiveRoot = cached.PrimitiveRoot;
+            return true;
+        }
+
+        if (TryLoadPersistedParameters(length, out modulus, out primitiveRoot))
+        {
+            NttParameterCache[length] = (modulus, primitiveRoot);
+            return true;
+        }
+
+        if (GenerateNttParameters(length, out modulus, out primitiveRoot))
+        {
+            NttParameterCache[length] = (modulus, primitiveRoot);
+            PersistParameters(length, modulus, primitiveRoot);
+            return true;
+        }
+
+        modulus = new GpuUInt128(0UL, 18446744069414584321UL); // 2^64 - 2^32 + 1
+        primitiveRoot = new GpuUInt128(0UL, 7UL);
+        return false;
+    }
+
+    private bool TryLoadPersistedParameters(int length, out GpuUInt128 modulus, out GpuUInt128 primitiveRoot)
+    {
+        if (!File.Exists(ParameterFilePath))
+        {
+            modulus = new GpuUInt128(0UL, 0UL);
+            primitiveRoot = new GpuUInt128(0UL, 0UL);
+            return false;
+        }
+
+        foreach (var line in File.ReadLines(ParameterFilePath))
+        {
+            if (string.IsNullOrWhiteSpace(line))
+            {
+                continue;
+            }
+
+            var parts = line.Split(' ', StringSplitOptions.RemoveEmptyEntries);
+            if (parts.Length != 5)
+            {
+                continue;
+            }
+
+            if (!int.TryParse(parts[0], out var storedLength) || storedLength != length)
+            {
+                continue;
+            }
+
+            ulong modHigh = ulong.Parse(parts[1]);
+            ulong modLow = ulong.Parse(parts[2]);
+            ulong rootHigh = ulong.Parse(parts[3]);
+            ulong rootLow = ulong.Parse(parts[4]);
+            modulus = new GpuUInt128(modHigh, modLow);
+            primitiveRoot = new GpuUInt128(rootHigh, rootLow);
+            return true;
+        }
+
+        modulus = new GpuUInt128(0UL, 0UL);
+        primitiveRoot = new GpuUInt128(0UL, 0UL);
+        return false;
+    }
+
+    private void PersistParameters(int length, GpuUInt128 modulus, GpuUInt128 primitiveRoot)
+    {
+        lock (ParameterFileLock)
+        {
+            using var writer = new StreamWriter(ParameterFilePath, append: true);
+            writer.WriteLine($"{length} {modulus.High} {modulus.Low} {primitiveRoot.High} {primitiveRoot.Low}");
+        }
+    }
+
+    private static readonly ulong[] SmallPrimes = GenerateSmallPrimes(1000);
+
+    private static ulong[] GenerateSmallPrimes(int limit)
+    {
+        var isComposite = new bool[limit + 1];
+        var primes = new List<ulong>();
+        for (int p = 2; p <= limit; p++)
+        {
+            if (isComposite[p])
+            {
+                continue;
+            }
+
+            primes.Add((ulong)p);
+            long square = (long)p * p;
+            for (long multiple = square; multiple <= limit; multiple += p)
+            {
+                isComposite[(int)multiple] = true;
+            }
+        }
+
+        return primes.ToArray();
+    }
+
+    private static bool GenerateNttParameters(int length, out GpuUInt128 modulus, out GpuUInt128 primitiveRoot)
+    {
+        ulong foundCandidate = 0UL;
+        ulong foundRoot = 0UL;
+        bool success = false;
+        var processorCount = Environment.ProcessorCount;
+        using var cts = new System.Threading.CancellationTokenSource();
+        object sync = new();
+
+        System.Threading.Tasks.Parallel.For(0, processorCount, (worker, state) =>
+        {
+            ulong k = (ulong)worker + 1UL;
+            while (!cts.IsCancellationRequested)
+            {
+                ulong candidate;
+                try
+                {
+                    candidate = checked((ulong)length * k + 1UL);
+                }
+                catch (OverflowException)
+                {
+                    break;
+                }
+
+                bool divisible = false;
+                foreach (ulong p in SmallPrimes)
+                {
+                    if (candidate % p == 0UL)
+                    {
+                        divisible = true;
+                        break;
+                    }
+                }
+
+                if (!divisible && IsPrime(candidate))
+                {
+                    ulong root = FindPrimitiveRoot(candidate, (ulong)length);
+                    lock (sync)
+                    {
+                        foundCandidate = candidate;
+                        foundRoot = root;
+                        success = true;
+                        cts.Cancel();
+                        state.Stop();
+                    }
+                    break;
+                }
+
+                k += (ulong)processorCount;
+            }
+        });
+
+        if (success)
+        {
+            modulus = new GpuUInt128(0UL, foundCandidate);
+            primitiveRoot = new GpuUInt128(0UL, foundRoot);
+            return true;
+        }
+
+        modulus = new GpuUInt128(0UL, 18446744069414584321UL);
+        primitiveRoot = new GpuUInt128(0UL, 7UL);
+        return false;
+    }
+
+    private static ulong FindPrimitiveRoot(ulong modulus, ulong order)
+    {
+        ulong phi = modulus - 1UL;
+        var factors = Factorize(phi);
+        for (ulong g = 2UL; g < modulus; g++)
+        {
+            bool ok = true;
+            foreach (ulong f in factors)
+            {
+                if (ModPow(g, phi / f, modulus) == 1UL)
+                {
+                    ok = false;
+                    break;
+                }
+            }
+
+            if (ok)
+            {
+                return ModPow(g, phi / order, modulus);
+            }
+        }
+
+        throw new InvalidOperationException("Primitive root not found.");
+    }
+
+    private static List<ulong> Factorize(ulong n)
+    {
+        var factors = new List<ulong>();
+        for (ulong p = 2UL; p * p <= n; p += p == 2UL ? 1UL : 2UL)
+        {
+            if (n % p == 0UL)
+            {
+                factors.Add(p);
+                while (n % p == 0UL)
+                {
+                    n /= p;
+                }
+            }
+        }
+
+        if (n > 1UL)
+        {
+            factors.Add(n);
+        }
+
+        return factors;
+    }
+
+    private static bool IsPrime(ulong n)
+    {
+        if (n < 2UL)
+        {
+            return false;
+        }
+
+        ulong[] bases = [2UL, 3UL, 5UL, 7UL, 11UL, 13UL];
+        ulong d = n - 1UL;
+        int s = 0;
+        while ((d & 1UL) == 0UL)
+        {
+            d >>= 1;
+            s++;
+        }
+
+        foreach (ulong a in bases)
+        {
+            if (a >= n)
+            {
+                continue;
+            }
+
+            ulong x = ModPow(a, d, n);
+            if (x == 1UL || x == n - 1UL)
+            {
+                continue;
+            }
+
+            bool cont = false;
+            for (int r = 1; r < s; r++)
+            {
+                x = MulMod(x, x, n);
+                if (x == n - 1UL)
+                {
+                    cont = true;
+                    break;
+                }
+            }
+
+            if (!cont)
+            {
+                return false;
+            }
+        }
+
+        return true;
+    }
+
+    private static ulong ModPow(ulong value, ulong exponent, ulong modulus)
+    {
+        ulong result = 1UL;
+        while (exponent > 0UL)
+        {
+            if ((exponent & 1UL) != 0UL)
+            {
+                result = MulMod(result, value, modulus);
+            }
+
+            value = MulMod(value, value, modulus);
+            exponent >>= 1;
+        }
+
+        return result;
+    }
+
+    private static ulong MulMod(ulong a, ulong b, ulong modulus)
+    {
+        return (ulong)(((UInt128)a * b) % modulus);
+    }
+
+    private static void AddSmallKernel(Index1D index, ulong add, ArrayView<GpuUInt128> value)
+    {
+        ulong carry = add;
+        for (int i = 0; i < value.Length && carry != 0UL; i++)
+        {
+            var limb = value[i];
+            ulong low = limb.Low + carry;
+            ulong carryOut = low < carry ? 1UL : 0UL;
+            ulong high = limb.High + carryOut;
+            carry = high < limb.High ? 1UL : 0UL;
+            value[i] = new GpuUInt128(high, low);
+        }
+    }
+
+    private static void SubtractSmallKernel(Index1D index, ulong subtract, ArrayView<GpuUInt128> value)
+    {
+        ulong borrow = subtract;
+        for (int i = 0; i < value.Length && borrow != 0UL; i++)
+        {
+            var limb = value[i];
+            ulong low = limb.Low;
+            ulong high = limb.High;
+            if (low >= borrow)
+            {
+                low -= borrow;
+                borrow = 0UL;
+            }
+            else
+            {
+                low = unchecked(low - borrow);
+                borrow = 1UL;
+                if (high != 0UL)
+                {
+                    high--;
+                    borrow = 0UL;
+                }
+                else
+                {
+                    high = ulong.MaxValue;
+                }
+            }
+
+            value[i] = new GpuUInt128(high, low);
+        }
+    }
+
+    private static void ReduceModMersenneKernel(Index1D index, ulong exponent, ArrayView<GpuUInt128> value)
+    {
+        int limbCount = (int)((exponent + 127UL) / 128UL);
+        for (int i = limbCount; i < value.Length; i++)
+        {
+            var limb = value[i];
+            if (limb.IsZero)
+            {
+                continue;
+            }
+
+            int target = i - limbCount;
+            var sum = value[target];
+            var original = sum;
+            sum = sum.Add(limb);
+            value[target] = sum;
+            value[i] = new GpuUInt128(0UL, 0UL);
+            if (sum.CompareTo(original) < 0)
+            {
+                int j = target + 1;
+                while (j < limbCount)
+                {
+                    var next = value[j];
+                    var nextOriginal = next;
+                    next = next.Add(1UL);
+                    value[j] = next;
+                    if (next.CompareTo(nextOriginal) >= 0)
+                    {
+                        break;
+                    }
+
+                    j++;
+                }
+            }
+        }
+
+        int topBits = (int)(exponent % 128UL);
+        if (topBits == 0)
+        {
+            topBits = 128;
+        }
+
+        int topIndex = limbCount - 1;
+        ulong topLowMask = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL;
+        ulong topHighMask = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL;
+        var top = value[topIndex];
+        var carryBits = top >> topBits;
+        top = new(top.High & topHighMask, top.Low & topLowMask);
+        value[topIndex] = top;
+
+        if (!carryBits.IsZero)
+        {
+            int j = 0;
+            while (j < limbCount && !carryBits.IsZero)
+            {
+                var cur = value[j];
+                var curOriginal = cur;
+                cur = cur.Add(carryBits);
+                value[j] = cur;
+                carryBits = cur.CompareTo(curOriginal) < 0 ? new GpuUInt128(0UL, 1UL) : new GpuUInt128(0UL, 0UL);
+                j++;
+            }
+        }
+
+        bool geq = true;
+        for (int i = topIndex; i >= 0; i--)
+        {
+            var limb = value[i];
+            ulong modHigh = i == topIndex ? topHighMask : ulong.MaxValue;
+            ulong modLow = i == topIndex ? topLowMask : ulong.MaxValue;
+            if (limb.High > modHigh || (limb.High == modHigh && limb.Low > modLow))
+            {
+                geq = true;
+                break;
+            }
+
+            if (limb.High < modHigh || (limb.High == modHigh && limb.Low < modLow))
+            {
+                geq = false;
+                break;
+            }
+        }
+
+        if (geq)
+        {
+            ulong borrow = 0UL;
+            for (int i = 0; i < limbCount - 1; i++)
+            {
+                var limb = value[i];
+                ulong newLow = limb.Low - ulong.MaxValue - borrow;
+                borrow = limb.Low < ulong.MaxValue + borrow ? 1UL : 0UL;
+                ulong newHigh = limb.High - ulong.MaxValue - borrow;
+                borrow = limb.High < ulong.MaxValue + borrow ? 1UL : 0UL;
+                value[i] = new GpuUInt128(newHigh, newLow);
+            }
+
+            var topLimb = value[topIndex];
+            ulong newTopLow = topLimb.Low - topLowMask - borrow;
+            borrow = topLimb.Low < topLowMask + borrow ? 1UL : 0UL;
+            ulong newTopHigh = topLimb.High - topHighMask - borrow;
+            value[topIndex] = new GpuUInt128(newTopHigh, newTopLow);
+        }
+
+        for (int i = limbCount; i < value.Length; i++)
+        {
+            value[i] = new GpuUInt128(0UL, 0UL);
+        }
+    }
+
+    private static void IsZeroKernel(Index1D index, ArrayView<GpuUInt128> value, ArrayView<byte> result)
+    {
+        byte isZero = 1;
+        for (int i = 0; i < value.Length; i++)
+        {
+            var limb = value[i];
+            if (limb.High != 0UL || limb.Low != 0UL)
+            {
+                isZero = 0;
+                break;
+            }
+        }
+
+        result[0] = isZero;
+    }
+
+    private static void ReduceModMersenne(Span<GpuUInt128> value, ulong exponent)
+    {
+        int limbCount = (int)((exponent + 127UL) / 128UL);
+        for (int i = limbCount; i < value.Length; i++)
+        {
+            var limb = value[i];
+            if (limb.IsZero)
+            {
+                continue;
+            }
+
+            int target = i - limbCount;
+            UInt128 original = value[target];
+            UInt128 sum = original + (UInt128)limb;
+            bool carry = sum < original;
+            value[target] = new GpuUInt128(sum);
+            value[i] = new GpuUInt128(0UL, 0UL);
+            if (carry)
+            {
+                int j = target + 1;
+                while (j < limbCount)
+                {
+                    UInt128 next = value[j];
+                    UInt128 nextSum = next + 1UL;
+                    value[j] = new GpuUInt128(nextSum);
+                    if (nextSum != 0UL)
+                    {
+                        break;
+                    }
+
+                    j++;
+                }
+            }
+        }
+
+        int topBits = (int)(exponent % 128UL);
+        if (topBits == 0)
+        {
+            topBits = 128;
+        }
+
+        int topIndex = limbCount - 1;
+        ulong topLowMask = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL;
+        ulong topHighMask = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL;
+        UInt128 mask = ((UInt128)topHighMask << 64) | topLowMask;
+        UInt128 top = value[topIndex];
+        UInt128 carryBits = top >> topBits;
+        top &= mask;
+        value[topIndex] = new GpuUInt128(top);
+        if (carryBits != 0UL)
+        {
+            int j = 0;
+            UInt128 carry = carryBits;
+            while (j < limbCount && carry != 0UL)
+            {
+                UInt128 cur = value[j];
+                UInt128 sum = cur + carry;
+                bool overflow = sum < cur;
+                value[j] = new GpuUInt128(sum);
+                carry = overflow ? 1UL : 0UL;
+                j++;
+            }
+        }
+
+        bool geq = true;
+        for (int i = topIndex; i >= 0; i--)
+        {
+            UInt128 limb = value[i];
+            UInt128 modLimb = i == topIndex ? mask : UInt128.MaxValue;
+            if (limb > modLimb)
+            {
+                geq = true;
+                break;
+            }
+
+            if (limb < modLimb)
+            {
+                geq = false;
+                break;
+            }
+        }
+
+        if (geq)
+        {
+            ulong borrow = 0UL;
+            for (int i = 0; i < limbCount - 1; i++)
+            {
+                ulong low = value[i].Low;
+                ulong high = value[i].High;
+                ulong newLow = low - ulong.MaxValue - borrow;
+                borrow = low < ulong.MaxValue + borrow ? 1UL : 0UL;
+                ulong newHigh = high - ulong.MaxValue - borrow;
+                borrow = high < ulong.MaxValue + borrow ? 1UL : 0UL;
+                value[i] = new GpuUInt128(newHigh, newLow);
+            }
+
+            ulong topLow = value[topIndex].Low;
+            ulong topHigh = value[topIndex].High;
+            ulong newTopLow = topLow - topLowMask - borrow;
+            borrow = topLow < topLowMask + borrow ? 1UL : 0UL;
+            ulong newTopHigh = topHigh - topHighMask - borrow;
+            value[topIndex] = new GpuUInt128(newTopHigh, newTopLow);
+        }
+
+        for (int i = limbCount; i < value.Length; i++)
+        {
+            value[i] = new GpuUInt128(0UL, 0UL);
+        }
+    }
+
+    private static void Kernel(Index1D index, ulong exponent, GpuUInt128 modulus, ArrayView<GpuUInt128> state)
+    {
+        // Lucas–Lehmer iteration in the field GF(2^p-1).
+        // For p < 128, use a Mersenne-specific squaring + reduction to avoid
+        // the generic 128-bit long-division-like reduction.
+        var s = new GpuUInt128(0UL, 4UL);
+        ulong i = 0UL, limit = exponent - 2UL;
+        if (exponent < 128UL)
+        {
+            for (; i < limit; i++)
+            {
+                s = SquareModMersenne128(s, exponent).SubMod(2UL, modulus);
+            }
+        }
+        else
+        {
+            for (; i < limit; i++)
+            {
+                s = s.SquareMod(modulus).SubMod(2UL, modulus);
+            }
+        }
+
+        state[0] = s;
+    }
+
+    // Squares a 128-bit value and reduces modulo M_p = 2^p - 1 using
+    // Mersenne folding. Valid for 0 < p < 128 and input s in [0, M_p].
+    private static GpuUInt128 SquareModMersenne128(GpuUInt128 s, ulong p)
+    {
+        // Full square: (H*2^64 + L)^2 = H^2*2^128 + 2HL*2^64 + L^2
+        SquareFull128(s, out var q3, out var q2, out var q1, out var q0);
+
+        // v = low128 + high128 (fold by 128 bits once)
+        var vHigh = new GpuUInt128(q3, q2);
+        var vLow = new GpuUInt128(q1, q0);
+        vLow = vLow.Add(vHigh);
+
+        // Mask to keep only p low bits and fold the carry bits (v >> p)
+        int topBits = (int)(p % 128UL);
+        if (topBits == 0)
+        {
+            topBits = 128;
+        }
+
+        ulong maskLow = topBits >= 64 ? ulong.MaxValue : (1UL << topBits) - 1UL;
+        ulong maskHigh = topBits > 64 ? (1UL << (topBits - 64)) - 1UL : 0UL;
+
+        var carryBits = vLow >> topBits;
+		vLow = new(vLow.High & maskHigh, vLow.Low & maskLow);
+
+		// Propagate carryBits back into the masked value (single-limb fold)
+		GpuUInt128 before;
+        while (!carryBits.IsZero)
+		{
+			before = vLow;
+			vLow = vLow.Add(carryBits);
+			carryBits = vLow.CompareTo(before) < 0 ? GpuUInt128.One : GpuUInt128.Zero;
+		}
+
+        // Final correction if v >= modulus (mask)
+        if (vLow.High > maskHigh || (vLow.High == maskHigh && vLow.Low >= maskLow))
+        {
+            ulong borrow = 0UL;
+            ulong newLow = vLow.Low - maskLow;
+            borrow = vLow.Low < maskLow ? 1UL : 0UL;
+            ulong newHigh = vLow.High - maskHigh - borrow;
+            vLow = new GpuUInt128(newHigh, newLow);
+        }
+
+        return vLow;
+    }
+
+    // Helper: 128-bit square -> 256-bit product as four 64-bit limbs (q3..q0)
+    private static void SquareFull128(GpuUInt128 value, out ulong q3, out ulong q2, out ulong q1, out ulong q0)
+    {
+        ulong L = value.Low;
+        ulong H = value.High;
+        Mul64Parts(L, L, out var hLL, out var lLL);    // L^2
+        Mul64Parts(H, H, out var hHH, out var lHH);    // H^2
+        Mul64Parts(L, H, out var hLH, out var lLH);    // L*H
+
+        // double LH
+        ulong dLH_low = lLH << 1;
+        ulong carry = (lLH >> 63) & 1UL;
+        ulong dLH_high = (hLH << 1) | carry;
+
+        q0 = lLL;
+        ulong sum1 = hLL + dLH_low;
+        ulong c1 = sum1 < hLL ? 1UL : 0UL;
+        q1 = sum1;
+
+        ulong sum2 = lHH + dLH_high;
+        ulong c2 = sum2 < lHH ? 1UL : 0UL;
+        sum2 += c1;
+        if (sum2 < c1)
+        {
+            c2++;
+        }
+        q2 = sum2;
+        q3 = hHH + c2;
+    }
+
+    // 64x64 -> 128 multiply (high,low)
+    [MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)]
+    private static void Mul64Parts(ulong a, ulong b, out ulong high, out ulong low)
+    {
+        ulong a0 = (uint)a;
+        ulong a1 = a >> 32;
+        ulong b0 = (uint)b;
+        ulong b1 = b >> 32;
+
+        ulong lo = a0 * b0;
+        ulong mid1 = a1 * b0;
+        ulong mid2 = a0 * b1;
+        ulong hi = a1 * b1;
+
+        ulong carry = (lo >> 32) + (uint)mid1 + (uint)mid2;
+        low = (lo & 0xFFFFFFFFUL) | (carry << 32);
+        hi += (mid1 >> 32) + (mid2 >> 32) + (carry >> 32);
+        high = hi;
+    }
+
+    [MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)]
+    private static ulong ModPow64(ulong value, ulong exponent, ulong modulus)
+    {
+        ulong result = 1UL;
+        ulong baseValue = value % modulus;
+        ulong exp = exponent;
+
+        while (exp != 0UL)
+        {
+            if ((exp & 1UL) != 0UL)
+            {
+                result = (ulong)(((UInt128)result * baseValue) % modulus);
+            }
+
+            baseValue = (ulong)(((UInt128)baseValue * baseValue) % modulus);
+            exp >>= 1;
+        }
+
+        return result;
+    }
+}
+

PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs

@@ -0,0 +1,63 @@
+using ILGPU.Runtime;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuOrder)
+{
+	private readonly GpuKernelType _kernelType = kernelType;
+	private readonly bool _useGpuOrder = useGpuOrder;
+
+	public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime)
+	{
+		using var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
+		var accelerator = gpuLease.Accelerator;
+		int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy and avoids TDR
+		UInt128 kStart = UInt128.One;
+		ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL;
+		byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
+
+        var kernel = _kernelType switch
+        {
+            GpuKernelType.Pow2Mod => gpuLease.Pow2ModOrderKernel,
+            _ => gpuLease.IncrementalOrderKernel,
+        };
+
+		using var foundBuffer = accelerator.Allocate1D<int>(1);
+		UInt128 remaining;
+		int currentSize;
+		int found = 0;
+		while (kStart <= maxK && Volatile.Read(ref isPrime))
+		{
+			remaining = maxK - kStart + 1UL;
+			currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
+			foundBuffer.MemSetToZero();
+			// Precompute residue automaton bases for this batch
+			UInt128 q0 = twoP * kStart + 1UL;
+			ulong q0m10 = q0.Mod10();
+			ulong q0m8 = q0.Mod8();
+			ulong q0m3 = q0.Mod3();
+			ulong q0m5 = q0.Mod5();
+			ulong step10 = ((exponent % 10UL) << 1) % 10UL;
+			ulong step8 = ((exponent & 7UL) << 1) & 7UL;
+			ulong step3 = ((exponent % 3UL) << 1) % 3UL;
+			ulong step5 = ((exponent % 5UL) << 1) % 5UL;
+			var ra = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
+
+
+			// Ensure device has small cycles table for early rejections in order kernels
+            var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
+            kernel(currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, divMul,
+                ra, foundBuffer.View, smallCyclesView);
+
+			accelerator.Synchronize();
+			foundBuffer.View.CopyToCPU(ref found, 1);
+			if (found != 0)
+			{
+				Volatile.Write(ref isPrime, false);
+				break;
+			}
+
+			kStart += (UInt128)currentSize;
+		}
+	}
+}

PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs

@@ -0,0 +1,83 @@
+using System.Buffers;
+using System.Runtime.InteropServices;
+using ILGPU.Runtime;
+
+namespace PerfectNumbers.Core.Gpu;
+
+public class MersenneNumberResidueGpuTester(bool useGpuOrder)
+{
+    private readonly bool _useGpuOrder = useGpuOrder;
+    // TODO: Integrate MersenneDivisorCycles.Shared to consult cycle lengths for small q (<= 4M) and fast-reject
+    // 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)
+    {
+        using var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
+        var accelerator = gpuLease.Accelerator;
+        // Ensure device has small cycles table for in-kernel lookup
+        var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
+        int batchSize = GpuConstants.ScanBatchSize;
+		UInt128 kStart = 1UL;
+		byte last = lastIsSeven ? (byte)1 : (byte)0;
+		var kernel = gpuLease.Pow2ModKernel;
+
+		using var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
+		ulong[] orderArray = ArrayPool<ulong>.Shared.Rent(batchSize);
+		UInt128 remaining;
+		int currentSize;
+		int i;
+		UInt128 batchSize128 = (UInt128)batchSize,
+				q = UInt128.Zero;
+		Span<ulong> orders = orderArray.AsSpan(0, batchSize);
+        try
+        {
+            while (kStart < maxK && Volatile.Read(ref isPrime))
+            {
+				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);
+				ulong step10 = ((exponent % 10UL) << 1) % 10UL;
+				ulong step8 = ((exponent & 7UL) << 1) & 7UL;
+				ulong step3 = ((exponent % 3UL) << 1) % 3UL;
+				ulong step5 = ((exponent % 5UL) << 1) % 5UL;
+
+                kernel(currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, 0UL,
+                    q0m10, q0m8, q0m3, q0m5, orderBuffer.View, smallCyclesView);
+
+				accelerator.Synchronize();
+				orderBuffer.View.CopyToCPU(ref MemoryMarshal.GetReference(orders), currentSize);
+				if (!Volatile.Read(ref isPrime))
+				{
+					return;
+				}
+
+                for (i = 0; i < currentSize; i++, q += twoP)
+                {
+                    if (orders[i] != 0UL && q.IsPrimeCandidate())
+                    {
+                        Volatile.Write(ref isPrime, false);
+                        break;
+                    }
+                }
+
+				kStart += batchSize128;
+			}
+		}
+		finally
+		{
+			ArrayPool<ulong>.Shared.Return(orderArray);
+		}
+	}
+}

PerfectNumbers.Core/InputParser.cs

@@ -0,0 +1,26 @@
+using System.Globalization;
+using System.Numerics;
+
+namespace PerfectNumbers.Core;
+
+public static class InputParser
+{
+    public static BigInteger ParseBigInteger(string value)
+    {
+        if (value.Contains('^'))
+        {
+            var parts = value.Split('^');
+            BigInteger b = BigInteger.Parse(parts[0], CultureInfo.InvariantCulture);
+            int e = int.Parse(parts[1], CultureInfo.InvariantCulture);
+            return BigInteger.Pow(b, e);
+        }
+
+        return BigInteger.Parse(value, CultureInfo.InvariantCulture);
+    }
+
+    public static decimal ParseDecimal(string value)
+    {
+        return decimal.Parse(value, CultureInfo.InvariantCulture);
+    }
+}
+

PerfectNumbers.Core/KRangeFinder.cs

@@ -0,0 +1,79 @@
+using PeterO.Numbers;
+
+namespace PerfectNumbers.Core;
+
+public static class KRangeFinder
+{
+    public static (int? Min, int? Max) FindForPrime(ERational alphaM, EInteger p, int kStart, int kEnd)
+    {
+        int? min = null;
+        int? max = null;
+        for (int k = kStart; k <= kEnd; k++)
+        {
+            ERational aP = AlphaCalculations.ComputeAlphaP(p, k);
+            ERational prod = aP.Multiply(alphaM);
+            int cmp = prod.CompareTo(RationalNumbers.Two);
+            if (cmp == 0)
+            {
+                if (min == null || k < min)
+                {
+                    min = k;
+                }
+                if (max == null || k > max)
+                {
+                    max = k;
+                }
+            }
+            if (cmp < 0)
+            {
+                break;
+            }
+        }
+        return (min, max);
+    }
+
+    public static bool TryFindForPrime(ERational alphaM, EInteger p, out int? min, out int? max, int kStart, int kEnd)
+    {
+        (min, max) = FindForPrime(alphaM, p, kStart, kEnd);
+        return min.HasValue && max.HasValue && min.Value <= max.Value;
+    }
+
+    public static (int? Min, int? Max) Find(ERational alphaM, EInteger pMin, EInteger pMax, int kStart, int kEnd, PrimeCache primes)
+    {
+        int? min = null;
+        int? max = null;
+        EInteger[] eulerPrimes = primes.GetEulerPrimes(pMin, pMax);
+        object sync = new();
+        Parallel.For(
+            0,
+            eulerPrimes.Length,
+            new ParallelOptions { MaxDegreeOfParallelism = Environment.ProcessorCount },
+            i =>
+            {
+                (int? localMin, int? localMax) = FindForPrime(alphaM, eulerPrimes[i], kStart, kEnd);
+                if (localMin == null)
+                {
+                    return;
+                }
+                lock (sync)
+                {
+                    if (min == null || localMin < min)
+                    {
+                        min = localMin;
+                    }
+                    if (max == null || localMax > max)
+                    {
+                        max = localMax;
+                    }
+                }
+            });
+        return (min, max);
+    }
+
+    public static bool TryFind(ERational alphaM, EInteger pMin, EInteger pMax, out int? min, out int? max, int kStart, int kEnd, PrimeCache primes)
+    {
+        (min, max) = Find(alphaM, pMin, pMax, kStart, kEnd, primes);
+        return min.HasValue && max.HasValue && min.Value <= max.Value;
+    }
+}
+

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -0,0 +1,381 @@
+using System.Buffers;
+using System.Text;
+using PerfectNumbers.Core.Gpu;
+using ILGPU;
+using ILGPU.Runtime;
+using System.Runtime.CompilerServices;
+using Open.Collections.Synchronized;
+
+namespace PerfectNumbers.Core;
+
+public class MersenneDivisorCycles
+{
+	private List<(ulong divisor, ulong cycleLength)> _table = [];
+	// Lightweight read-mostly cache for small divisors (<= 4,000,000). 0 => unknown
+	private uint[]? _smallCycles;
+
+	public static MersenneDivisorCycles Shared { get; } = new MersenneDivisorCycles();
+
+	public MersenneDivisorCycles()
+	{
+	}
+
+	public void LoadFrom(string path)
+	{
+		EnsureSmallBuffer();
+		List<(ulong divisor, ulong cycle)> cycles = [];
+		using Stream outputStream = new FileStream(path, FileMode.Open, FileAccess.Read, FileShare.Read, BufferSize10M, useAsync: false);
+		foreach (var (d, c) in EnumerateStream(outputStream))
+		{
+			cycles.Add((d, c));
+			if (d <= PerfectNumberConstants.MaxQForDivisorCycles)
+			{
+				_smallCycles![(int)d] = (uint)Math.Min(uint.MaxValue, c == 0UL ? 1UL : c);
+			}
+		}
+
+		cycles.Sort((a, b) => a.divisor < b.divisor ? -1 : (a.divisor > b.divisor ? 1 : 0));
+		_table = cycles;
+	}
+
+	public const int SmallDivisorsMax = PerfectNumberConstants.MaxQForDivisorCycles;
+
+	// Provides a snapshot array of small cycles [0..PerfectNumberConstants.MaxQForDivisorCycles]. Index is the divisor.
+	// The returned array is a copy safe for use across threads and device uploads.
+	public uint[] ExportSmallCyclesSnapshot()
+	{
+		EnsureSmallBuffer();
+		uint[] snapshot = new uint[PerfectNumberConstants.MaxQForDivisorCycles + 1];
+		Array.Copy(_smallCycles!, snapshot, snapshot.Length);
+		return snapshot;
+	}
+
+	public static IEnumerable<(ulong divisor, ulong cycleLength)> EnumerateStream(Stream compressor)
+	{
+		// Binary pairs: (ulong divisor, ulong cycle)
+		using var reader = new BinaryReader(compressor, Encoding.UTF8, leaveOpen: true);
+		while (true)
+		{
+			ulong d, c;
+			try
+			{
+				d = reader.ReadUInt64();
+				c = reader.ReadUInt64();
+			}
+			catch (EndOfStreamException)
+			{
+				yield break;
+			}
+
+			yield return (d, c);
+		}
+	}
+
+	public ulong GetCycle(ulong divisor)
+	{
+		// Fast-path: in-memory array for small divisors
+		if (divisor <= (ulong)PerfectNumberConstants.MaxQForDivisorCycles)
+		{
+			var arr = _smallCycles;
+			if (arr is not null)
+			{
+				uint cached = arr[(int)divisor];
+				if (cached != 0U)
+				{
+					return cached;
+				}
+			}
+		}
+		
+		// Binary search for divisor in sorted _table
+		int mid, right, left = 0;
+		ulong d;
+		List<(ulong divisor, ulong cycleLength)> cycles = _table;
+
+		right = cycles.Count - 1;
+		while (left <= right)
+		{
+			mid = left + ((right - left) >> 1);
+			d = cycles[mid].divisor;
+
+			if (d == divisor)
+			{
+				return cycles[mid].cycleLength;
+			}
+
+			if (d < divisor)
+			{
+				left = mid + 1;
+			}
+			else
+			{
+				right = mid - 1;
+			}
+		}
+
+		return CalculateCycleLength(divisor);
+	}
+
+	public static UInt128 GetCycle(UInt128 divisor)
+	{
+		// For divisor = 2^k, cycle is 1
+		if ((divisor & (divisor - UInt128.One)) == UInt128.Zero)
+		{
+			return UInt128.One;
+		}
+
+		// Otherwise, find order of 2 mod divisor
+		UInt128 order = UInt128.One, pow = UInt128Numbers.Two;
+		while (pow != UInt128.One)
+		{
+			pow <<= 1;
+			if (pow >= divisor)
+				pow -= divisor;
+
+			order++;
+		}
+
+		return order;
+	}
+
+	public static (long nextPosition, long completeCount) FindLast(string path)
+	{
+		using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false);
+		long count = 0;
+		foreach (var _ in EnumerateStream(outputStream))
+		{
+			count++;
+		}
+		return (outputStream.Position, count);
+	}
+
+	private const int BufferSize10M = 10 * 1024 * 1024;
+
+	public static void Generate(string path, ulong maxDivisor, int threads = 16)
+	{
+		// Binary-only generation of (divisor, cycle) pairs. CSV scaffolding removed.
+		ulong start, d = start = 2;
+		ulong blockSize = maxDivisor / (ulong)threads + 1;
+		Task[] tasks = new Task[threads];
+		using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: false);
+		// Start fresh to avoid mixing formats
+		outputStream.SetLength(0L);
+		outputStream.Position = 0L;
+		using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false);
+		for (var taskIndex = 0; taskIndex < threads; taskIndex++)
+		{
+			Task t = Task.Factory.StartNew((taskIndex) =>
+			{
+				ulong cycle, d = start + (ulong)(int)taskIndex! * blockSize;
+				(ulong d, ulong cycle)[] localCycles = ArrayPool<(ulong, ulong)>.Shared.Rent((int)(maxDivisor - d + 1UL));
+				int localCycleIndex = 0;
+				for (; d <= maxDivisor; d++)
+				{
+					if ((d & 1UL) == 0UL)
+					{
+						// skip even divisors
+						continue;
+					}
+
+					cycle = CalculateCycleLength(d);
+
+					localCycles[localCycleIndex++] = (d, cycle);
+				}
+
+				if (localCycleIndex > 0)
+				{
+					lock (writer)
+					{
+						writer.BaseStream.Position = writer.BaseStream.Length;
+						for (int i = 0; i < localCycleIndex; i++)
+						{
+							(d, cycle) = localCycles[i];
+							if (d % 3UL == 0UL || d % 5UL == 0UL || d % 7UL == 0UL || d % 11UL == 0UL)
+								continue;
+
+							writer.Write(d);
+							writer.Write(cycle);
+						}
+
+						writer.Flush();
+					}
+				}
+
+				ArrayPool<(ulong, ulong)>.Shared.Return(localCycles, clearArray: false);
+
+			}, taskIndex);
+
+			tasks[taskIndex] = t;
+		}
+
+		Task.WaitAll(tasks);
+	}
+
+	public static void GenerateGpu(string path, ulong maxDivisor, int batchSize = 1_000_000, long skipCount = 0L, long nextPosition = 0L)
+	{
+		// Prepare output file with header
+
+		ulong start = 3UL;
+		var locker = new object();
+
+		var pool = ArrayPool<ulong>.Shared;
+		ulong batchSizeUL = (ulong)batchSize, d, end;
+		int count = (int)Math.Min(batchSizeUL, maxDivisor), i, idx;
+		ulong[] divisors, outCycles, validDivisors;
+
+		divisors = pool.Rent(count);
+		outCycles = pool.Rent(count);
+		using var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
+		var accelerator = lease.Accelerator;
+		using var stream = accelerator.CreateStream();
+		var kernel = accelerator.LoadAutoGroupedStreamKernel<
+			Index1D,
+			ArrayView1D<ulong, Stride1D.Dense>,
+			ArrayView1D<ulong, Stride1D.Dense>>(GpuDivisorCycleKernel);
+
+
+		using Stream outputStream = new FileStream(path, FileMode.OpenOrCreate, FileAccess.ReadWrite, FileShare.Read, BufferSize10M, useAsync: true);
+		if (nextPosition > 0L)
+		{
+			outputStream.Position = nextPosition;
+		}
+		else
+		{
+			outputStream.SetLength(0L);
+			outputStream.Position = 0L;
+		}
+
+		using var writer = new BinaryWriter(outputStream, Encoding.UTF8, leaveOpen: false);
+		try
+		{
+			while (skipCount > 0L && start <= maxDivisor)
+			{
+				if ((start & 1UL) == 0UL || start % 3UL == 0UL || start % 5UL == 0UL || start % 7UL == 0UL || start % 11UL == 0UL)
+				{
+					start++;
+					continue;
+				}
+
+				start++;
+				skipCount--;
+			}
+
+			while (start <= maxDivisor)
+			{
+				end = Math.Min(start + batchSizeUL - 1UL, maxDivisor);
+				count = checked((int)(end - start + 1UL));
+
+				idx = 0;
+				for (d = start; d <= end; d++)
+				{
+					if ((d & 1UL) == 0UL)
+						continue;
+
+					divisors[idx++] = d;
+				}
+
+				validDivisors = divisors[..idx];
+
+				// Use GpuKernelPool to get a kernel and context
+				// GpuKernelPool.Run((accelerator, stream) =>
+				var bufferDiv = accelerator.Allocate1D(validDivisors);
+				var bufferCycle = accelerator.Allocate1D<ulong>(idx);
+
+				kernel(
+					count,
+					bufferDiv.View,
+					bufferCycle.View);
+
+				accelerator.Synchronize();
+
+				bufferCycle.View.CopyToCPU(ref outCycles[0], idx);
+				bufferCycle.Dispose();
+				bufferDiv.Dispose();
+
+				// Collect results in order
+				for (i = 0; i < idx; i++)
+				{
+					d = divisors[i];
+
+					if (d % 3UL == 0UL || d % 5UL == 0UL || d % 7UL == 0UL || d % 11UL == 0UL)
+						continue;
+
+					writer.Write(d);
+					writer.Write(outCycles[i]);
+				}
+
+				writer.Flush();
+				Console.WriteLine($"...processed divisor = {d}");
+				start = end + 1UL;
+			}
+		}
+		finally
+		{
+			pool.Return(divisors, clearArray: false);
+			pool.Return(outCycles, clearArray: false);
+		}
+	}
+
+	const byte ByteZero = 0;
+	const byte ByteOne = 1;
+
+	// GPU kernel for divisor cycle calculation
+	static void GpuDivisorCycleKernel(
+		Index1D index,
+		ArrayView1D<ulong, Stride1D.Dense> divisors,
+		ArrayView1D<ulong, Stride1D.Dense> outCycles)
+	{
+		int i = index.X;
+		outCycles[i] = CalculateCycleLengthGpu(divisors[i]);
+	}
+
+	// GPU-friendly version of cycle length calculation
+	public static ulong CalculateCycleLengthGpu(ulong divisor)
+	{
+		if ((divisor & (divisor - 1UL)) == 0UL)
+			return 1UL;
+
+		ulong order = 1UL, pow = 2UL;
+		while (pow != 1UL)
+		{
+			pow <<= 1;
+			if (pow >= divisor)
+				pow -= divisor;
+
+			order++;
+		}
+
+		return order;
+	}
+
+	public static ulong CalculateCycleLength(ulong divisor)
+	{
+		// For divisor = 2^k, cycle is 1
+		if ((divisor & (divisor - 1UL)) == 0UL)
+		{
+			return 1UL;
+		}
+
+		// Otherwise, find order of 2 mod divisor
+		ulong order = 1UL, pow = 2UL;
+		while (pow != 1)
+		{
+			pow <<= 1;
+			if (pow >= divisor)
+				pow -= divisor;
+
+			order++;
+		}
+
+		return order;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private void EnsureSmallBuffer()
+	{
+		if (_smallCycles is null)
+		{
+			_smallCycles = new uint[PerfectNumberConstants.MaxQForDivisorCycles + 1];
+		}
+	}
+}

PerfectNumbers.Core/MersenneNumberTester.cs

@@ -0,0 +1,312 @@
+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;
+
+namespace PerfectNumbers.Core;
+
+public enum GpuKernelType
+{
+    Incremental,
+    Pow2Mod,
+}
+
+public sealed class MersenneNumberTester(
+    bool useIncremental = true,
+    bool useOrderCache = false,
+    GpuKernelType kernelType = GpuKernelType.Incremental,
+    bool useModuloWorkaround = false,
+    bool useOrder = false,
+    bool useGpuLucas = true,
+    bool useGpuScan = true,
+    bool useGpuOrder = true,
+	bool useResidue = true,
+    ulong maxK = 5_000_000UL)
+{
+    private readonly bool _useIncremental = useIncremental;
+    private readonly bool _useResidue = useResidue;
+	private readonly ulong _maxK = maxK;
+	private readonly GpuKernelType _kernelType = kernelType;
+    private readonly bool _useModuloWorkaround = useModuloWorkaround;
+    private readonly bool _useGpuLucas = useGpuLucas;
+    private readonly bool _useGpuScan = useGpuScan;     // device for pow2mod/incremental scanning
+    private readonly bool _useGpuOrder = useGpuOrder;   // device for order computations
+    private static readonly ConcurrentDictionary<UInt128, ulong> OrderCache = new();
+	private readonly 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.
+
+    public void WarmUpOrders(ulong exponent, ulong limit = 5_000_000UL)
+	{
+		bool cacheEnabled = useOrderCache;
+		if (!cacheEnabled)
+		{
+			return;
+		}
+
+		UInt128 twoP = (UInt128)exponent << 1; // 2 * p
+		bool lastIsSeven = (exponent & 3UL) == 3UL;
+		// 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();
+				bool shouldCheck = lastIsSeven
+					? 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.CalculateOrder();
+				OrderCache[q] = ord;
+			}
+
+			return;
+		}
+
+		UInt128[] qsBuffer = ArrayPool<UInt128>.Shared.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 = lastIsSeven
+				? 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)
+		{
+			ArrayPool<UInt128>.Shared.Return(qsBuffer);
+			return;
+		}
+
+                UInt128[] qs = qsBuffer;
+
+                using var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
+                var accelerator = gpuLease.Accelerator;
+        var orderKernel = gpuLease.OrderKernel!;
+
+		// 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;
+		while (offset < idx)
+		{
+            int count = Math.Min(batchSize, idx - offset);
+            using var qBuffer = accelerator.Allocate1D<PerfectNumbers.Core.Gpu.GpuUInt128>(count);
+            using var orderBuffer = accelerator.Allocate1D<ulong>(count);
+            // Convert to device-friendly GpuUInt128 on the fly
+            var tmp = System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Rent(count);
+            for (int i = 0; i < count; i++)
+            {
+                tmp[i] = (PerfectNumbers.Core.Gpu.GpuUInt128)qs[offset + i];
+            }
+            // Copy only the populated elements to the device buffer to avoid overrun on pooled arrays.
+            qBuffer.View.CopyFromCPU(ref tmp[0], count);
+            System.Buffers.ArrayPool<PerfectNumbers.Core.Gpu.GpuUInt128>.Shared.Return(tmp);
+            orderKernel(count, exponent, divMul, qBuffer.View, orderBuffer.View);
+            accelerator.Synchronize();
+
+			ulong[] orders = orderBuffer.GetAsArray1D();
+			for (int i = 0; i < count; i++)
+			{
+				ulong order = orders[i];
+				if (order == 0UL)
+				{
+					order = qs[offset + i].CalculateOrder();
+				}
+
+				if (useOrderCache)
+				{
+					OrderCache[qs[offset + i]] = order;
+				}
+			}
+
+			offset += count;
+		}
+
+		ArrayPool<UInt128>.Shared.Return(qs);
+	}
+
+    public bool IsMersennePrime(ulong exponent)
+    {
+        // Safe early rejections using known orders for tiny primes:
+        // 7 | M_p iff 3 | p. Avoid rejecting p == 3 where M_3 == 7 is prime.
+        if ((exponent % 3UL) == 0UL && exponent != 3UL)
+        {
+            return false;
+        }
+
+        if ((exponent & 3UL) == 1UL && exponent.SharesFactorWithExponentMinusOne())
+        {
+            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 ^ exponent) / 2 + 1;
+        bool lastIsSeven = (exponent & 3UL) == 3UL;
+
+		// This is how maxK was for residue test before changes UInt128 maxK = maxKOverride ?? ResidueMaxK;
+		if (_useResidue)
+		{
+			// Disable this is if we don't want to hard limit maxK
+			// maxK = _maxK > maxK ? maxK : _maxK;
+
+			if (_useGpuScan)
+			{
+				_residueGpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime);
+			}
+			else
+			{
+				_residueCpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime);
+			}
+
+			return prePrime;
+
+			// if (!prePrime)
+			// {
+			// 	return false;
+			// }
+		}
+
+		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;
+			bool lastIsSevenPre = (exponent & 3UL) == 3UL;
+			if (_useGpuScan)
+			{
+				_incrementalGpuTester.Scan(exponent, twoPPre, lastIsSevenPre, maxKPre, ref prePrime);
+			}
+			else
+			{
+				_incrementalCpuTester.Scan(exponent, twoPPre, lastIsSevenPre, maxKPre, ref prePrime);
+			}
+
+			if (!prePrime)
+			{
+				return false;
+			}
+
+			return _useGpuLucas
+				? _lucasLehmerGpuTester.IsPrime(exponent, runOnGpu: true)
+				: _lucasLehmerCpuTester.IsPrime(exponent);
+		}
+
+        bool isPrime = true;
+
+        if (useOrder)
+        {
+            if (_useGpuOrder)
+            {
+                _orderGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime);
+            }
+            else
+            {
+                _orderCpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime);
+            }
+        }
+        else
+        {
+            if (_useGpuScan)
+            {
+                _incrementalGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime);
+            }
+            else
+            {
+                _incrementalCpuTester.Scan(exponent, twoP, lastIsSeven, 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/MersennePrimeFactorTester.cs

@@ -0,0 +1,558 @@
+using System.Collections.Concurrent;
+using System.Runtime.CompilerServices;
+using System.Numerics;
+using UInt128 = System.UInt128;
+
+namespace PerfectNumbers.Core;
+
+public static class MersennePrimeFactorTester
+{
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool IsPrimeFactor(ulong p, ulong q, CancellationToken ct)
+    {
+        if (q < 2UL || (q & 1UL) == 0)
+        {
+            return false;
+        }
+
+        if ((q & 7UL) != 1UL && (q & 7UL) != 7UL)
+        {
+            return false;
+        }
+
+        if (ct.IsCancellationRequested)
+        {
+            return false;
+        }
+
+        if (!PrimeTester.IsPrimeInternal(q, ct))
+        {
+            return false;
+        }
+
+        if (p < 2UL || (p & 1UL) == 0)
+        {
+            return false;
+        }
+
+        ulong value = checked(p << 1);
+        if ((q - 1UL) % value != 0UL)
+        {
+            return false;
+        }
+
+        ulong order = GetOrderOf2ModPrime(q, ct);
+        return !ct.IsCancellationRequested && order == p;
+    }
+
+    private static readonly Dictionary<UInt128, UInt128> _orderCache = [];
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static ulong GetOrderOf2ModPrime(ulong q, CancellationToken ct) =>
+        (ulong)GetOrderOf2ModPrime((UInt128)q, ct);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static UInt128 GetOrderOf2ModPrime(UInt128 q, CancellationToken ct) =>
+        OrderOf2ModPrime(q, ct);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static UInt128 OrderOf2ModPrime(UInt128 q, CancellationToken ct)
+    {
+        lock (_orderCache)
+        {
+            if (_orderCache.TryGetValue(q, out var cached))
+            {
+                return cached;
+            }
+        }
+
+        UInt128 phi = q - 1;
+        UInt128[] factors = phi <= ulong.MaxValue
+            ? Array.ConvertAll(Factor64((ulong)phi, ct), x => (UInt128)x)
+            : Factor128(phi, ct);
+        UInt128 order = phi;
+
+        UInt128 candidate, prime;
+        int i, factorsCount = factors.Length;
+        for (i = 0; i < factorsCount; i++)
+        {
+            prime = factors[i];
+            if (ct.IsCancellationRequested)
+            {
+                return 0;
+            }
+
+            while (order % prime == 0)
+            {
+                if (ct.IsCancellationRequested)
+                {
+                    return 0;
+                }
+
+                candidate = order / prime;
+                if (ModPow128(2, candidate, q, ct) == 1)
+                {
+                    order = candidate;
+                }
+                else
+                {
+                    break;
+                }
+            }
+        }
+
+        if (!ct.IsCancellationRequested)
+        {
+            lock (_orderCache)
+            {
+                _orderCache[q] = order;
+            }
+        }
+
+        return order;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong ModPow64(ulong value, ulong exponent, ulong modulus, CancellationToken ct)
+    {
+        ulong result = 1UL;
+        ulong baseValue = value % modulus;
+        ulong exp = exponent;
+
+        while (exp != 0UL)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                return 0UL;
+            }
+
+            if ((exp & 1UL) != 0UL)
+            {
+                result = MulMod64(result, baseValue, modulus);
+            }
+
+            if (ct.IsCancellationRequested)
+            {
+                return 0UL;
+            }
+
+            baseValue = MulMod64(baseValue, baseValue, modulus);
+            exp >>= 1;
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong MulMod64(ulong a, ulong b, ulong mod) => (ulong)((UInt128)a * b % mod);
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static UInt128 ModPow128(UInt128 value, UInt128 exponent, UInt128 modulus, CancellationToken ct)
+    {
+        UInt128 result = 1;
+        UInt128 baseValue = value % modulus;
+        UInt128 exp = exponent;
+
+        while (exp != 0)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                return 0;
+            }
+
+            if ((exp & 1) != 0)
+            {
+                result = MulMod128(result, baseValue, modulus);
+            }
+
+            if (ct.IsCancellationRequested)
+            {
+                return 0;
+            }
+
+            baseValue = MulMod128(baseValue, baseValue, modulus);
+            exp >>= 1;
+        }
+
+        return result;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 mod)
+    {
+        BigInteger product = (BigInteger)a * b;
+        product %= (BigInteger)mod;
+        return (UInt128)product;
+    }
+
+    private static readonly ConcurrentBag<List<ulong>> ListPool = new();
+    private static readonly ConcurrentBag<Dictionary<ulong, byte>> DictPool = new();
+    private static readonly Dictionary<ulong, ulong[]> FactorCache = new();
+    private static readonly object FactorCacheLock = new();
+
+    private static List<ulong> RentList() => ListPool.TryTake(out var list) ? list : new List<ulong>(4);
+
+    private static void ReturnList(List<ulong> list)
+    {
+        list.Clear();
+        ListPool.Add(list);
+    }
+
+    private static Dictionary<ulong, byte> RentDictionary() => DictPool.TryTake(out var dict) ? dict : new Dictionary<ulong, byte>();
+
+    private static void ReturnDictionary(Dictionary<ulong, byte> dict)
+    {
+        dict.Clear();
+       DictPool.Add(dict);
+    }
+
+    private static readonly ConcurrentBag<Dictionary<UInt128, byte>> DictPool128 = new();
+
+    private static Dictionary<UInt128, byte> RentDictionary128() =>
+        DictPool128.TryTake(out var dict) ? dict : new Dictionary<UInt128, byte>();
+
+    private static void ReturnDictionary128(Dictionary<UInt128, byte> dict)
+    {
+        dict.Clear();
+        DictPool128.Add(dict);
+    }
+
+    private static ulong[] Factor64(ulong n, CancellationToken ct)
+    {
+        // lock (FactorCacheLock)
+        // {
+            // if (FactorCache.TryGetValue(n, out var cached))
+            // {
+            //     return cached;
+            // }
+        // }
+
+        if (n == 1UL)
+        {
+            return [];
+        }
+
+        var dict = RentDictionary();
+        var queue = new Queue<ulong>();
+        queue.Enqueue(n);
+
+        while (queue.Count != 0)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                break;
+            }
+
+            ulong m = queue.Dequeue();
+            if (m == 1UL)
+            {
+                continue;
+            }
+
+            if (PrimeTester.IsPrimeInternal(m, ct))
+            {
+                dict[m] = 1;
+                continue;
+            }
+
+            if (ct.IsCancellationRequested)
+            {
+                break;
+            }
+
+            ulong d = PollardRho(m, ct);
+            if (ct.IsCancellationRequested)
+            {
+                break;
+            }
+
+            queue.Enqueue(d);
+            queue.Enqueue(m / d);
+        }
+
+        ulong[] result = [.. dict.Keys];
+        // lock (FactorCacheLock)
+        // {
+        //     FactorCache[n] = result;
+        // }
+
+        ReturnDictionary(dict);
+        return result;
+    }
+
+    private static UInt128[] Factor128(UInt128 n, CancellationToken ct)
+    {
+        if (n == 1)
+        {
+            return [];
+        }
+
+        var dict = RentDictionary128();
+        var queue = new Queue<UInt128>();
+        queue.Enqueue(n);
+
+        while (queue.Count != 0)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                break;
+            }
+
+            UInt128 m = queue.Dequeue();
+            if (m == 1)
+            {
+                continue;
+            }
+
+            if (IsPrime128(m, ct))
+            {
+                dict[m] = 1;
+                continue;
+            }
+
+            if (ct.IsCancellationRequested)
+            {
+                break;
+            }
+
+            UInt128 d = PollardRho128(m, ct);
+            if (ct.IsCancellationRequested)
+            {
+                break;
+            }
+
+            queue.Enqueue(d);
+            queue.Enqueue(m / d);
+        }
+
+        var result = new UInt128[dict.Count];
+        int index = 0;
+        foreach (var key in dict.Keys)
+        {
+            result[index++] = key;
+        }
+
+        ReturnDictionary128(dict);
+        return result;
+    }
+
+    private static ulong PollardRho(ulong n, CancellationToken ct)
+    {
+        if ((n & 1UL) == 0UL)
+        {
+            return 2UL;
+        }
+
+        var rng = new Random(1234567);
+
+        ulong c, d, x, y, diff;
+        while (true)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                return 0UL;
+            }
+
+            c = (ulong)rng.Next(2, int.MaxValue);
+            x = (ulong)rng.Next(2, int.MaxValue);
+            y = x;
+            d = 1UL;
+
+            while (d == 1UL)
+            {
+                if (ct.IsCancellationRequested)
+                {
+                    return 0UL;
+                }
+
+                x = (MulMod64(x, x, n) + c) % n;
+                if (ct.IsCancellationRequested)
+                {
+                    return 0UL;
+                }
+
+                y = (MulMod64(y, y, n) + c) % n;
+                if (ct.IsCancellationRequested)
+                {
+                    return 0UL;
+                }
+
+                y = (MulMod64(y, y, n) + c) % n;
+                diff = x > y ? x - y : y - x;
+                if (ct.IsCancellationRequested)
+                {
+                    return 0UL;
+                }
+
+                d = Gcd(diff, n, ct);
+            }
+
+            if (d != n)
+            {
+                return d;
+            }
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Gcd(ulong a, ulong b, CancellationToken ct)
+    {
+        ulong t;
+        while (b != 0UL)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                return 1UL; // We're returning 1 to avoid division by 0 in any place
+            }
+
+            t = a % b;
+            a = b;
+            b = t;
+        }
+
+        return a;
+    }
+
+    private static bool IsPrime128(UInt128 n, CancellationToken ct)
+    {
+        if (n <= 3)
+        {
+            return n >= 2;
+        }
+
+        if ((n & 1) == 0)
+        {
+            return false;
+        }
+
+        BigInteger bn = (BigInteger)n;
+        BigInteger d = bn - 1;
+        int s = 0;
+        while ((d & 1) == 0)
+        {
+            d >>= 1;
+            s++;
+        }
+
+        BigInteger[] bases = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];
+        for (int i = 0; i < bases.Length; i++)
+        {
+            BigInteger a = bases[i];
+            if (a >= bn)
+            {
+                continue;
+            }
+
+            if (ct.IsCancellationRequested)
+            {
+                return false;
+            }
+
+            BigInteger x = BigInteger.ModPow(a, d, bn);
+            if (x == 1 || x == bn - 1)
+            {
+                continue;
+            }
+
+            bool cont = false;
+            for (int r = 1; r < s; r++)
+            {
+                if (ct.IsCancellationRequested)
+                {
+                    return false;
+                }
+
+                x = BigInteger.ModPow(x, 2, bn);
+                if (x == bn - 1)
+                {
+                    cont = true;
+                    break;
+                }
+            }
+
+            if (!cont)
+            {
+                return false;
+            }
+        }
+
+        return true;
+    }
+
+    private static UInt128 PollardRho128(UInt128 n, CancellationToken ct)
+    {
+        if ((n & 1) == 0)
+        {
+            return 2;
+        }
+
+        var rng = new Random(1234567);
+
+        UInt128 c, d, x, y, diff;
+        while (true)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                return 0;
+            }
+
+            c = (UInt128)rng.NextInt64(2, long.MaxValue);
+            x = (UInt128)rng.NextInt64(2, long.MaxValue);
+            y = x;
+            d = 1;
+
+            while (d == 1)
+            {
+                if (ct.IsCancellationRequested)
+                {
+                    return 0;
+                }
+
+                x = (MulMod128(x, x, n) + c) % n;
+                if (ct.IsCancellationRequested)
+                {
+                    return 0;
+                }
+
+                y = (MulMod128(y, y, n) + c) % n;
+                if (ct.IsCancellationRequested)
+                {
+                    return 0;
+                }
+
+                y = (MulMod128(y, y, n) + c) % n;
+                diff = x > y ? x - y : y - x;
+                if (ct.IsCancellationRequested)
+                {
+                    return 0;
+                }
+
+                d = Gcd128(diff, n, ct);
+            }
+
+            if (d != n)
+            {
+                return d;
+            }
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static UInt128 Gcd128(UInt128 a, UInt128 b, CancellationToken ct)
+    {
+        UInt128 t;
+        while (b != 0)
+        {
+            if (ct.IsCancellationRequested)
+            {
+                return 1;
+            }
+
+            t = a % b;
+            a = b;
+            b = t;
+        }
+
+        return a;
+    }
+}

PerfectNumbers.Core/MersenneResidueStepper.cs

@@ -0,0 +1,69 @@
+namespace PerfectNumbers.Core;
+
+public class MersenneResidueStepper
+{
+    private UInt128 m;
+    private ulong p;
+    private UInt128 pow2_mod;
+    private UInt128 M_mod;
+
+    public MersenneResidueStepper(UInt128 modulus, ulong p0, UInt128 M_mod_m_at_p0)
+    {
+        if (modulus == 0)
+            throw new ArgumentException("Modulus must be positive.");
+        m = modulus;
+        p = p0;
+        M_mod = M_mod_m_at_p0 % m;
+        pow2_mod = (M_mod + 1) % m;
+    }
+
+    public (ulong, UInt128) State() => (p, M_mod);
+
+    public (ulong, UInt128) Step(long delta)
+    {
+        if (delta == 0)
+            return State();
+        UInt128 g = UInt128.One;
+        ulong absDelta = (ulong)Math.Abs(delta);
+        for (ulong i = 0; i < absDelta; i++)
+        {
+            g <<= 1;
+            if (g >= m)
+                g -= m;
+        }
+        if (delta > 0)
+            pow2_mod = (pow2_mod * g) % m;
+        else
+        {
+            // Modular inverse for negative delta
+            pow2_mod = (pow2_mod * ModInverse(g, m)) % m;
+        }
+        p = (ulong)((long)p + delta);
+        M_mod = (pow2_mod + m - 1) % m;
+        return State();
+    }
+
+    public (ulong, UInt128) To(ulong p_target)
+    {
+        long delta = (long)p_target - (long)p;
+        return Step(delta);
+    }
+
+    private static UInt128 ModInverse(UInt128 a, UInt128 m)
+    {
+        UInt128 m0 = m, t, q;
+        UInt128 x0 = 0, x1 = 1;
+        if (m == 1) return 0;
+        while (a > 1)
+        {
+            q = a / m;
+            t = m;
+            m = a % m; a = t;
+            t = x0;
+            x0 = x1 - q * x0;
+            x1 = t;
+        }
+        if (x1 < 0) x1 += m0;
+        return x1;
+    }
+}

PerfectNumbers.Core/ModResidueTracker.cs

@@ -0,0 +1,313 @@
+using System.Runtime.CompilerServices;
+
+using UInt128 = System.UInt128;
+
+namespace PerfectNumbers.Core;
+
+public enum ResidueModel
+{
+	Identity,   // residue of N modulo d
+	Mersenne,   // residue of (2^N - 1) modulo d
+}
+
+/// Tracks residues for a monotonically increasing sequence of numbers against
+/// a dynamic set of divisors. Designed to support repeated divisor checks
+/// while the tested number grows step-by-step without ever decreasing.
+public sealed class ModResidueTracker
+{
+	public static readonly ModResidueTracker Shared = new(ResidueModel.Mersenne, initialNumber: PerfectNumberConstants.BiggestKnownEvenPerfectP, initialized: true);
+
+	private readonly ResidueModel _model;
+	private UInt128 _currentNumber;
+	private bool _initialized;
+	// Per-divisor sorted state at _currentNumber. Aligned arrays for speed.
+	private readonly List<UInt128> _divisors = new(1024);   // sorted ascending
+	private readonly List<UInt128> _residues = new(1024);   // same index as _divisors
+	private int _cursor;                                    // merge pointer for ascending candidate stream
+	private bool _hasLastAppended;
+	private UInt128 _lastAppended;
+
+	public ModResidueTracker(ResidueModel model, UInt128 initialNumber = default, bool initialized = false)
+	{
+		_model = model;
+		_currentNumber = initialNumber;
+		_initialized = initialized;
+	}
+
+	/// Returns true if the provided number (either p or M_p depending on model)
+	/// is divisible by the divisor. Advances internal state forward if needed
+	/// (number must be non-decreasing across calls). Adds and initializes new
+	/// divisors on first use.
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public bool IsDivisible(UInt128 number, UInt128 divisor)
+	{
+		if (divisor == UInt128.Zero)
+		{
+			throw new ArgumentOutOfRangeException(nameof(divisor));
+		}
+
+		EnsureAt(number);
+
+		List<UInt128> divisors = _divisors;
+		int idx = LowerBound(divisors, divisor);
+		if (idx == divisors.Count || divisors[idx] != divisor)
+		{
+			UInt128 residue = ComputeInitialResidue(number, divisor);
+			divisors.Insert(idx, divisor);
+			_residues.Insert(idx, residue);
+			if (_hasLastAppended && idx == divisors.Count - 1 && divisor >= _lastAppended)
+			{
+				_lastAppended = divisor;
+			}
+			else if (!_hasLastAppended)
+			{
+				_hasLastAppended = true;
+				_lastAppended = divisor;
+			}
+
+			return residue == UInt128.Zero;
+		}
+
+		return _residues[idx] == UInt128.Zero;
+	}
+
+	// Merge-walk support: prepare for increasing sequence of candidate divisors.
+	public void BeginMerge(UInt128 number)
+	{
+		EnsureAt(number);
+		_cursor = 0;
+	}
+
+	// Merge or append current candidate divisor (ascending order required across calls).
+	// Returns true and sets 'divisible' based on residue at 'number'.
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public bool MergeOrAppend(UInt128 number, UInt128 divisor, out bool divisible)
+	{
+		// Assume number is non-decreasing; advance residues lazily when number changes.
+		List<UInt128> divisors = _divisors;
+		if (number != _currentNumber)
+		{
+			EnsureAt(number);
+			// After EnsureAt we can keep _cursor (still valid monotonic pointer)
+			if (_cursor > divisors.Count) _cursor = divisors.Count;
+		}
+
+		UInt128 residue;
+		// Fast path: append to tail when candidates are strictly increasing.
+		List<UInt128> residues = _residues;
+		if (divisors.Count == 0 || (_hasLastAppended && divisor >= _lastAppended && divisors[^1] <= divisor))
+		{
+			residue = ComputeInitialResidue(number, divisor);
+			divisors.Add(divisor);
+			residues.Add(residue);
+			_hasLastAppended = true;
+			_lastAppended = divisor;
+			_cursor = divisors.Count; // position after last
+			divisible = residue == UInt128.Zero;
+			return true;
+		}
+
+		// Merge pointer: advance until current >= divisor
+		var cursor = _cursor;
+		while (cursor < divisors.Count && divisors[cursor] < divisor)
+		{
+			cursor++;
+		}
+
+		if (cursor < divisors.Count && divisors[cursor] == divisor)
+		{
+			divisible = residues[cursor] == UInt128.Zero;
+			_cursor = cursor + 1; // position after matched
+			return true;
+		}
+
+
+		// Insert at cursor to keep order
+		residue = ComputeInitialResidue(number, divisor);
+		divisors.Insert(cursor, divisor);
+		residues.Insert(cursor, residue);
+		divisible = residue == UInt128.Zero;
+		_cursor = cursor + 1;
+
+		// update last appended hint
+		if (!_hasLastAppended || divisor >= _lastAppended)
+		{
+			_hasLastAppended = true;
+			_lastAppended = divisor;
+		}
+
+		return true;
+	}
+
+	/// Ensures the internal state corresponds to the provided number by moving
+	/// forward from the current number (never backwards). O(Count(divisors)).
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public void EnsureAt(UInt128 number)
+	{
+		if (!_initialized)
+		{
+			_currentNumber = number;
+			// Initialize residues lazily upon first check for each divisor.
+			_initialized = true;
+			return;
+		}
+
+		if (number == _currentNumber)
+		{
+			return;
+		}
+
+		if (number < _currentNumber)
+		{
+			throw new InvalidOperationException("Number must be non-decreasing.");
+		}
+
+		UInt128 delta = number - _currentNumber;
+		List<UInt128> divisors = _divisors;
+		int divisorsCount = divisors.Count;
+		if (divisorsCount == 0)
+		{
+			_currentNumber = number;
+			return;
+		}
+
+		// Advance all tracked residues from _currentNumber to 'number'.
+		// For Identity: r' = (r + delta) mod d
+		// For Mersenne: r' = (2^delta * r + (2^delta - 1)) mod d
+		// Optimize small deltas with simple forms when possible.
+		int i = 0;
+		UInt128 d, value;
+		List<UInt128> residues = _residues;
+
+		if (_model == ResidueModel.Identity)
+		{
+			// Avoid 128-bit division: delta is small in scanner (<= 6) so use add/sub.
+			for (; i < divisorsCount; i++)
+			{
+				d = divisors[i];
+				value = residues[i] + delta;
+				// For tiny moduli (2..6) we may need a few subtractions; keep the loop tight.
+				while (value >= d)
+				{
+					value -= d;
+				}
+
+				residues[i] = value;
+			}
+		}
+		else
+		{
+			for (; i < divisors.Count; i++)
+			{
+				d = divisors[i];
+				value = PowMod128(UInt128Numbers.Two, delta, d);
+				// r' = pow2Delta*r + (pow2Delta - 1) mod d
+				value = MulMod128(value, residues[i], d) + value - UInt128.One;
+				if (value >= d)
+				{
+					value -= d;
+				}
+
+				residues[i] = value;
+			}
+		}
+
+		_currentNumber = number;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static int LowerBound(List<UInt128> list, UInt128 value)
+	{
+		int lo = 0, hi = list.Count, mid;
+		while (lo < hi)
+		{
+			mid = (int)(((uint)lo + (uint)hi) >> 1);
+			if (list[mid] < value)
+			{
+				lo = mid + 1;
+			}
+			else
+			{
+				hi = mid;
+			}
+		}
+
+		return lo;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private UInt128 ComputeInitialResidue(UInt128 number, UInt128 divisor)
+	{
+		if (_model == ResidueModel.Identity)
+		{
+			// Fast path for 64-bit inputs
+			if (number <= ulong.MaxValue && divisor <= ulong.MaxValue)
+			{
+				return (ulong)number % (ulong)divisor;
+			}
+
+			return number % divisor;
+		}
+
+		// residue of Mersenne number M_number = 2^number - 1 modulo divisor
+		UInt128 pow = PowMod128(UInt128Numbers.Two, number, divisor);
+
+		// q divides M_p ⇔ 2^p ≡ 1 (mod q)  AND  ord_q(2) | p
+		if (pow != UInt128.One)
+			return (pow - UInt128.One) % divisor;
+
+		ulong ord = divisor.CalculateOrder();
+		if (ord == 0UL || divisor % ord != UInt128.Zero)
+			return UInt128.One; // does NOT divide
+
+		return UInt128.Zero; // q really divides M_p
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static UInt128 PowMod128(UInt128 baseValue, UInt128 exponent, UInt128 modulus)
+	{
+		if (modulus == UInt128.One)
+		{
+			return UInt128.Zero;;
+		}
+
+		UInt128 result = UInt128.One;
+		UInt128 value = baseValue % modulus;
+
+		while (exponent != UInt128.Zero)
+		{
+			if ((exponent & UInt128.One) != UInt128.Zero)
+			{
+				result = MulMod128(result, value, modulus);
+			}
+
+			exponent >>= 1;
+			if (exponent != UInt128.Zero)
+			{
+				value = MulMod128(value, value, modulus);
+			}
+		}
+
+		return result;
+	}
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static UInt128 MulMod128(UInt128 a, UInt128 b, UInt128 modulus)
+	{
+		// Same approach as MulMod in MersenneNumberTester: double-and-add to avoid overflow.
+		UInt128 result = UInt128.Zero;
+		a %= modulus;
+
+		while (b != UInt128.Zero)
+		{
+			if ((b & UInt128.One) != UInt128.Zero)
+			{
+				result += result < modulus ? a : a - modulus;
+			}
+
+			a = a < modulus ? a << 1 : (a << 1) - modulus;
+			b >>= 1;
+		}
+
+		return result;
+	}
+}

PerfectNumbers.Core/ModuloAutomata.cs

@@ -0,0 +1,215 @@
+using System.Runtime.CompilerServices;
+using UInt128 = System.UInt128;
+
+namespace PerfectNumbers.Core;
+
+/// <summary>
+/// Tracks residues for a fixed-step arithmetic progression without divisions.
+/// Optimized variant for the Mersenne scan q = 2*p*k + 1 that needs mod 10, 8, 3 and 5.
+/// </summary>
+public sealed class MersenneResidueAutomaton
+{
+    private readonly UInt128 _step;               // 2*p
+    private UInt128 _currentQ;                     // current q value
+
+    // Cached residues for current q
+    public ulong Mod10R { get; private set; }
+    public ulong Mod8R  { get; private set; }
+    public ulong Mod3R  { get; private set; }
+    public ulong Mod5R  { get; private set; }
+
+    // Cached residues for step (2*p mod m)
+    private readonly ulong _step10;
+    private readonly ulong _step8;
+    private readonly ulong _step3;
+    private readonly ulong _step5;
+
+    public MersenneResidueAutomaton(ulong exponent)
+    {
+        _step = (UInt128)exponent << 1; // 2 * p
+        _currentQ = _step + 1UL;        // start at k = 1
+
+        // init residues
+        Mod10R = Mod10(_currentQ);
+        Mod8R  = Mod8(_currentQ);
+        Mod3R  = Mod3(_currentQ);
+        Mod5R  = Mod5(_currentQ);
+
+        _step10 = Mod10((ulong)(_step % 10UL));
+        _step8  = (ulong)(_step & 7UL); // % 8
+        _step3  = Mod3(_step);
+        _step5  = Mod5(_step);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public UInt128 CurrentQ() => _currentQ;
+
+    /// <summary>
+    /// Advance to next k (q += 2*p), updating residues branchlessly.
+    /// </summary>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public void Next()
+    {
+        _currentQ += _step;
+
+        ulong r10 = Mod10R + _step10;
+        r10 -= (r10 >= 10UL) ? 10UL : 0UL;
+        Mod10R = r10;
+
+        ulong r8 = Mod8R + _step8;
+        r8 &= 7UL; // % 8
+        Mod8R = r8;
+
+        ulong r3 = Mod3R + _step3;
+        r3 -= (r3 >= 3UL) ? 3UL : 0UL;
+        Mod3R = r3;
+
+        ulong r5 = Mod5R + _step5;
+        if (r5 >= 5UL)
+        {
+            r5 -= 5UL;
+        }
+        Mod5R = r5;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Mod10(UInt128 value)
+    {
+        // Split and fold under mod 10: 2^64 ≡ 6 (mod 10)
+        ulong high = (ulong)(value >> 64);
+        ulong low = (ulong)value;
+        ulong highRem = Mod10(high);
+        ulong lowRem = Mod10(low);
+        ulong combined = lowRem + highRem * 6UL;
+        return Mod10(combined);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Mod10(ulong value) => value % 10UL;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Mod8(UInt128 value) => (ulong)value & 7UL;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Mod3(UInt128 value)
+    {
+        // 2^64 ≡ 1 (mod 3)
+        ulong high = (ulong)(value >> 64);
+        ulong low = (ulong)value;
+        ulong rem = (low % 3UL) + (high % 3UL);
+        return rem >= 3UL ? rem - 3UL : rem;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Mod5(UInt128 value)
+    {
+        // 2^64 ≡ 1 (mod 5)
+        ulong high = (ulong)(value >> 64);
+        ulong low = (ulong)value;
+        ulong rem = (low % 5UL) + (high % 5UL);
+        return rem >= 5UL ? rem - 5UL : rem;
+    }
+}
+
+/// <summary>
+/// Generic automaton for scanning integers ending with digit 7: n = 10*t + 7.
+/// Maintains residues r_m without divisions while advancing by +10 each step.
+/// </summary>
+public sealed class Ending7Automaton
+{
+    private ulong _n; // current candidate n, ends with 7
+
+    private readonly ulong[] _moduli;
+    private readonly ulong[] _steps;   // 10 % m
+    private readonly ulong[] _res;     // n % m
+    private readonly bool[] _reachable;
+
+    public Ending7Automaton(ulong start, params ulong[] moduli)
+    {
+        if (moduli == null || moduli.Length == 0)
+        {
+            _moduli = Array.Empty<ulong>();
+            _steps = Array.Empty<ulong>();
+            _res = Array.Empty<ulong>();
+            _reachable = Array.Empty<bool>();
+            _n = AlignTo7(start);
+            return;
+        }
+
+        _n = AlignTo7(start);
+        int len = moduli.Length;
+        _moduli = new ulong[len];
+        _steps = new ulong[len];
+        _res = new ulong[len];
+        _reachable = new bool[len];
+        for (int i = 0; i < len; i++)
+        {
+            ulong m = moduli[i];
+            _moduli[i] = m;
+            ulong step = 10UL % m;
+            _steps[i] = step;
+            _res[i] = _n % m; // one-time division at construction
+            // Zero is reachable iff gcd(step, m) | res
+            ulong g = Gcd(step, m);
+            _reachable[i] = (g == 0UL) ? false : (_res[i] % g == 0UL);
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public ulong Current() => _n;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public bool DivisibleBy(ulong modulus)
+    {
+        for (int i = 0; i < _moduli.Length; i++)
+        {
+            if (_moduli[i] == modulus)
+            {
+                return _reachable[i] && _res[i] == 0UL;
+            }
+        }
+
+        return (_n % modulus) == 0UL;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public void Next()
+    {
+        _n += 10UL;
+        for (int i = 0; i < _moduli.Length; i++)
+        {
+            if (!_reachable[i])
+            {
+                continue;
+            }
+
+            ulong m = _moduli[i];
+            ulong r = _res[i] + _steps[i];
+            if (r >= m)
+            {
+                r -= m;
+            }
+            _res[i] = r;
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong AlignTo7(ulong start)
+    {
+        ulong r = start % 10UL;
+        ulong add = (7UL + 10UL - r) % 10UL;
+        return start + add;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong Gcd(ulong a, ulong b)
+    {
+        while (b != 0UL)
+        {
+            ulong t = a % b;
+            a = b;
+            b = t;
+        }
+        return a;
+    }
+}

PerfectNumbers.Core/PerfectNumberConstants.cs

@@ -0,0 +1,11 @@
+namespace PerfectNumbers.Core;
+
+public static class PerfectNumberConstants
+{
+	public const ulong BiggestKnownEvenPerfectP = 136279841UL;
+	public static readonly ulong PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000;
+
+	// Maximum K for residue-based Mersenne test (q = 2*p*k + 1 up to this k)
+	public static readonly UInt128 ResidueMaxK = 5_000_000UL;
+    public const int MaxQForDivisorCycles = 4_000_000;
+}

PerfectNumbers.Core/PerfectNumbers.Core.csproj

@@ -0,0 +1,23 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <PropertyGroup>
+    <TargetFramework>net8.0</TargetFramework>
+    <ImplicitUsings>enable</ImplicitUsings>
+    <Nullable>enable</Nullable>
+    <LangVersion>12.0</LangVersion>
+  </PropertyGroup>
+
+  <ItemGroup>
+    <PackageReference Include="PeterO.Numbers" Version="1.8.1" />
+    <PackageReference Include="Open.Numeric.Primes" Version="4.0.4" />
+    <PackageReference Include="ILGPU" Version="1.5.3" />
+    <PackageReference Include="ILGPU.Algorithms" Version="1.5.3" />
+    <PackageReference Include="lz4" Version="1.3.1.2" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <InternalsVisibleTo Include="EvenPerfectBitScanner.Tests" />
+    <InternalsVisibleTo Include="PerfectNumbers.Core.Tests" />
+  </ItemGroup>
+
+</Project>

PerfectNumbers.Core/PowerCache.cs

@@ -0,0 +1,65 @@
+using System.Collections.Concurrent;
+using System.Numerics;
+using PeterO.Numbers;
+
+namespace PerfectNumbers.Core;
+
+public static class PowerCache
+{
+    private static readonly ConcurrentDictionary<ulong, BigInteger[]> BigCache = new();
+    private static readonly ConcurrentDictionary<EInteger, EInteger[]> ECache = new();
+
+    public static BigInteger Get(ulong baseVal, int exponent)
+    {
+        if (!BigCache.TryGetValue(baseVal, out var powers))
+        {
+            powers = new BigInteger[exponent + 1];
+            powers[0] = BigInteger.One;
+            BigCache[baseVal] = powers;
+        }
+
+        if (powers.Length <= exponent)
+        {
+            int oldLength = powers.Length;
+            Array.Resize(ref powers, exponent + 1);
+            for (int i = oldLength; i <= exponent; i++)
+                powers[i] = BigInteger.Zero;
+
+            BigCache[baseVal] = powers;
+        }
+
+        if (powers[exponent].IsZero)
+        {
+            for (int i = 1; i <= exponent; i++)
+                if (powers[i].IsZero)
+                    powers[i] = powers[i - 1] * baseVal;
+        }
+        
+        return powers[exponent];
+    }
+
+    public static EInteger Get(EInteger baseVal, int exponent)
+    {
+        if (!ECache.TryGetValue(baseVal, out var powers))
+        {
+            powers = new EInteger[exponent + 1];
+            powers[0] = EInteger.One;
+            ECache[baseVal] = powers;
+        }
+        if (powers.Length <= exponent)
+        {
+            int oldLength = powers.Length;
+            Array.Resize(ref powers, exponent + 1);
+            for (int i = oldLength; i <= exponent; i++)
+                powers[i] = null;
+            ECache[baseVal] = powers;
+        }
+        if (powers[exponent] == null)
+        {
+            for (int i = 1; i <= exponent; i++)
+                if (powers[i] == null)
+                    powers[i] = powers[i - 1]! * baseVal;
+        }
+        return powers[exponent]!;
+    }
+}

PerfectNumbers.Core/PrimeCache.cs

@@ -0,0 +1,96 @@
+using System.Buffers;
+using PeterO.Numbers;
+using Open.Numeric.Primes;
+
+namespace PerfectNumbers.Core;
+
+public sealed class PrimeCache
+{
+    private readonly List<EInteger> _primes = new();
+    private readonly List<byte> _primeMod4 = new();
+    private IEnumerator<ulong> _enumerator = Prime.Numbers.GetEnumerator();
+    private static readonly EInteger[] EmptyEulerPrimesArray = Array.Empty<EInteger>();
+
+    private void EnsurePrimesUpTo(EInteger limit)
+    {
+        EInteger last = _primes.Count > 0 ? _primes[^1] : EInteger.Zero;
+        if (last.CompareTo(limit) >= 0)
+        {
+            return;
+        }
+
+        while (last.CompareTo(limit) < 0 && _enumerator.MoveNext())
+        {
+            ulong current = _enumerator.Current;
+            last = current;
+            _primes.Add(last);
+            _primeMod4.Add((byte)(current & 3));
+        }
+    }
+
+    public EInteger[] GetEulerPrimes(EInteger start, EInteger end)
+    {
+        EnsurePrimesUpTo(end);
+
+        int index = _primes.BinarySearch(start);
+        if (index < 0)
+        {
+            index = ~index;
+        }
+
+        EInteger[] primesArray = ArrayPool<EInteger>.Shared.Rent(_primes.Count);
+        try
+        {
+            Span<EInteger> primes = new(primesArray);
+            int count = 0;
+            for (; index < _primes.Count; index++)
+            {
+                EInteger val = _primes[index];
+                if (val.CompareTo(end) > 0)
+                {
+                    break;
+                }
+                if (_primeMod4[index] == 1)
+                {
+                    primes[count++] = val;
+                }
+            }
+
+            return count > 0 ? primes[..count].ToArray() : EmptyEulerPrimesArray;
+        }
+        finally
+        {
+            ArrayPool<EInteger>.Shared.Return(primesArray);
+        }
+    }
+
+    public IEnumerable<EInteger> EnumeratePrimes(EInteger start)
+    {
+        ulong s = (ulong)start;
+        EnsurePrimesUpTo(start);
+
+        int index = _primes.BinarySearch(s);
+        if (index < 0)
+        {
+            index = ~index;
+        }
+
+        for (; index < _primes.Count; index++)
+        {
+            yield return _primes[index];
+        }
+
+        while (_enumerator.MoveNext())
+        {
+            ulong current = _enumerator.Current;
+            EInteger val = EInteger.FromUInt64(current);
+            _primes.Add(val);
+            _primeMod4.Add((byte)(current & 3));
+            if (current >= s)
+            {
+                yield return val;
+            }
+        }
+    }
+}
+

PerfectNumbers.Core/PrimeTester.cs

@@ -0,0 +1,371 @@
+using System.Buffers;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using ILGPU;
+using ILGPU.Algorithms;
+using ILGPU.Runtime;
+using PerfectNumbers.Core.Gpu;
+
+namespace PerfectNumbers.Core;
+
+public sealed class PrimeTester(bool useInternal = false)
+{
+    public bool IsPrime(ulong n, CancellationToken ct) => IsPrimeInternal(n, ct);
+
+    // Optional GPU-assisted primality: batched small-prime sieve on device.
+    // TODO: Consider extending with deterministic Miller–Rabin bases for 64-bit range when stable on ILGPU.
+    public bool IsPrimeGpu(ulong n, CancellationToken ct)
+    {
+        if (GpuContextPool.ForceCpu)
+        {
+            return IsPrimeInternal(n, ct);
+        }
+
+        Span<ulong> one = stackalloc ulong[1];
+        Span<byte> outFlags = stackalloc byte[1];
+        one[0] = n;
+        IsPrimeBatchGpu(one, outFlags);
+        return outFlags[0] != 0;
+    }
+
+    // Note: GPU-backed primality path is implemented via IsPrimeGpu/IsPrimeBatchGpu and is routed
+    // from EvenPerfectBitScanner based on --primes-device. Follow-up: add deterministic MR rounds.
+
+    internal static bool IsPrimeInternal(ulong n, CancellationToken ct)
+    {
+        bool result = true;
+        if (n <= 3UL)
+        {
+            result = n >= 2UL;
+        }
+        else if ((n & 1UL) == 0)
+        {
+            result = false;
+        }
+        else if (n % 5 == 0)
+        {
+            result = false;
+        }
+        else
+        {
+            if ((n % 10UL) == 1UL && SharesFactorWithMaxExponent(n))
+            {
+                result = false;
+            }
+            else
+            {
+                var smallPrimeDivisorsLength = PrimesGenerator.SmallPrimes.Length;
+				ulong[] smallPrimeDivisors = PrimesGenerator.SmallPrimes;
+				ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2;
+                for (int i = 0; i < smallPrimeDivisorsLength; i++)
+                {
+                    if (smallPrimeDivisorsMul[i] > n)
+                    {
+                        break;
+                    }
+
+					if (n % smallPrimeDivisors[i] == 0)
+                    {
+                        result = false;
+                        break;
+                    }
+                }
+            }
+        }
+
+        return result;
+    }
+
+    public static int GpuBatchSize { get; set; } = 262_144;
+
+    public static void IsPrimeBatchGpu(ReadOnlySpan<ulong> values, Span<byte> results)
+    {
+        // Limit concurrency and declare variables outside loops for performance and reuse
+        using var limiter = GpuPrimeWorkLimiter.Acquire();
+        using var gpu = GpuContextPool.Rent();
+        var accelerator = gpu.Accelerator;
+
+        // Get per-accelerator cached kernel and device primes buffer.
+        var state = GpuKernelState.GetOrCreate(accelerator);
+
+        MemoryBuffer1D<ulong, Stride1D.Dense> input;
+        MemoryBuffer1D<byte, Stride1D.Dense> output;
+        int totalLength, batchSize, pos, remaining, count;
+        ulong[] temp;
+
+        totalLength = values.Length;
+        batchSize = Math.Max(1, GpuBatchSize);
+
+        // Rent per-call scratch buffers sized to current batch (thread-safe)
+        var scratch = state.RentScratch(batchSize, accelerator);
+        input = scratch.Input;
+        output = scratch.Output;
+
+        // Reusable host buffer per call
+        temp = ArrayPool<ulong>.Shared.Rent(batchSize);
+
+        pos = 0;
+        while (pos < totalLength)
+        {
+            remaining = totalLength - pos;
+            count = remaining > batchSize ? batchSize : remaining;
+
+            // Copy chunk to host temp then to device
+            values.Slice(pos, count).CopyTo(temp);
+            input.View.CopyFromCPU(ref temp[0], count);
+
+            state.Kernel(count, input.View, state.DevicePrimes.View, output.View);
+            accelerator.Synchronize();
+            output.View.CopyToCPU(ref results[pos], count);
+
+            pos += count;
+        }
+
+        ArrayPool<ulong>.Shared.Return(temp);
+        state.ReturnScratch(scratch);
+    }
+
+    // Per-accelerator GPU state for prime sieve (kernel + uploaded primes).
+    private sealed class KernelState
+    {
+        public Action<Index1D, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> Kernel { get; }
+        public MemoryBuffer1D<ulong, Stride1D.Dense> DevicePrimes { get; }
+        private readonly Accelerator _accel;
+        private readonly System.Collections.Concurrent.ConcurrentBag<ScratchBuffers> _scratchPool = new();
+
+        public KernelState(Accelerator accelerator)
+        {
+            _accel = accelerator;
+            // Compile once per accelerator and upload primes once.
+            Kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(SmallPrimeSieveKernel);
+
+            var primes = PrimesGenerator.SmallPrimes;
+            DevicePrimes = accelerator.Allocate1D<ulong>(primes.Length);
+            DevicePrimes.View.CopyFromCPU(primes);
+        }
+
+        internal sealed class ScratchBuffers : IDisposable
+        {
+            public MemoryBuffer1D<ulong, Stride1D.Dense> Input { get; private set; }
+            public MemoryBuffer1D<byte, Stride1D.Dense> Output { get; private set; }
+            public int Capacity { get; private set; }
+
+            public ScratchBuffers(Accelerator accel, int capacity)
+            {
+                Capacity = Math.Max(1, capacity);
+                Input = accel.Allocate1D<ulong>(Capacity);
+                Output = accel.Allocate1D<byte>(Capacity);
+            }
+
+            public void Dispose()
+            {
+                Output.Dispose();
+                Input.Dispose();
+            }
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public ScratchBuffers RentScratch(int minCapacity, Accelerator accel)
+        {
+            minCapacity = Math.Max(1, minCapacity);
+            while (_scratchPool.TryTake(out var sb))
+            {
+                if (sb.Capacity >= minCapacity)
+                {
+                    return sb;
+                }
+
+                sb.Dispose();
+            }
+
+            return new ScratchBuffers(accel, minCapacity);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void ReturnScratch(ScratchBuffers scratch)
+        {
+            _scratchPool.Add(scratch);
+        }
+
+        public void Clear()
+        {
+            while (_scratchPool.TryTake(out var sb))
+            {
+                sb.Dispose();
+            }
+            DevicePrimes.Dispose();
+        }
+    }
+
+    private static class GpuKernelState
+    {
+        // Map accelerator to cached state; use Lazy to serialize kernel creation
+        private static readonly System.Collections.Concurrent.ConcurrentDictionary<Accelerator, Lazy<KernelState>> States = new();
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static KernelState GetOrCreate(Accelerator accelerator)
+        {
+            var lazy = States.GetOrAdd(
+                accelerator,
+                acc => new Lazy<KernelState>(() => new KernelState(acc), System.Threading.LazyThreadSafetyMode.ExecutionAndPublication));
+            return lazy.Value;
+        }
+
+        public static void Clear(Accelerator accelerator)
+        {
+            if (States.TryRemove(accelerator, out var lazy))
+            {
+                if (lazy.IsValueCreated)
+                {
+                    var state = lazy.Value;
+                    state.Clear();
+                }
+            }
+        }
+    }
+
+    // Expose cache clearing for accelerator disposal coordination
+    public static void ClearGpuCaches(Accelerator accelerator)
+    {
+        GpuKernelState.Clear(accelerator);
+    }
+
+    // GPU kernel: small-prime sieve only. Returns 1 if passes sieve (probable prime), 0 otherwise.
+    // TODO: Follow-up: consider integration with MR rounds for bases {2,3,5,7,11,13} to make it deterministic for 64-bit.
+    private static void SmallPrimeSieveKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<ulong> smallPrimes, ArrayView<byte> results)
+    {
+        ulong n = numbers[index];
+        if (n <= 3UL)
+        {
+            results[index] = (byte)(n >= 2UL ? 1 : 0);
+            return;
+        }
+
+        if ((n & 1UL) == 0UL || (n % 5UL) == 0UL)
+        {
+            results[index] = 0;
+            return;
+        }
+
+        // Early reject special gcd heuristic with floor(log2 n)
+        ulong m = 63UL - (ulong)ILGPU.Algorithms.XMath.LeadingZeroCount(n);
+        if (BinaryGcdGpu(n, m) != 1UL)
+        {
+            results[index] = 0;
+            return;
+        }
+
+        long len = smallPrimes.Length;
+        for (long i = 0; i < len; i++)
+        {
+            ulong p = smallPrimes[(int)i];
+			if (p * p > n)
+			{
+				break;
+			}
+			
+			if ((n % p) == 0UL)
+			{
+				results[index] = 0;
+				return;
+			}
+        }
+
+        results[index] = 1;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    internal static bool SharesFactorWithMaxExponent(ulong n)
+    {
+        ulong m = (ulong)BitOperations.Log2(n);
+        return BinaryGcd(n, m) != 1UL;
+    }
+
+    internal static void SharesFactorWithMaxExponentBatch(ReadOnlySpan<ulong> values, Span<byte> results)
+    {
+        using var gpu = GpuContextPool.Rent();
+        var accelerator = gpu.Accelerator;
+        var kernel = accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<byte>>(SharesFactorKernel);
+
+        int length = values.Length;
+        using var inputBuffer = accelerator.Allocate1D<ulong>(length);
+        using var resultBuffer = accelerator.Allocate1D<byte>(length);
+
+        ulong[] temp = ArrayPool<ulong>.Shared.Rent(length);
+        values.CopyTo(temp);
+        inputBuffer.View.CopyFromCPU(ref temp[0], length);
+        kernel(length, inputBuffer.View, resultBuffer.View);
+        accelerator.Synchronize();
+        resultBuffer.View.CopyToCPU(ref results[0], length);
+        ArrayPool<ulong>.Shared.Return(temp);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    private static ulong BinaryGcd(ulong u, ulong v)
+    {
+        if (u == 0UL)
+        {
+            return v;
+        }
+
+        if (v == 0UL)
+        {
+            return u;
+        }
+
+        int shift = BitOperations.TrailingZeroCount(u | v);
+        u >>= BitOperations.TrailingZeroCount(u);
+
+        do
+        {
+            v >>= BitOperations.TrailingZeroCount(v);
+            if (u > v)
+            {
+                (u, v) = (v, u);
+            }
+
+            v -= u;
+        }
+        while (v != 0UL);
+
+        return u << shift;
+    }
+
+    private static void SharesFactorKernel(Index1D index, ArrayView<ulong> numbers, ArrayView<byte> results)
+    {
+        ulong n = numbers[index];
+        ulong m = 63UL - (ulong)XMath.LeadingZeroCount(n);
+        ulong gcd = BinaryGcdGpu(n, m);
+        results[index] = gcd == 1UL ? (byte)0 : (byte)1;
+    }
+
+    private static ulong BinaryGcdGpu(ulong u, ulong v)
+    {
+        if (u == 0UL)
+        {
+            return v;
+        }
+
+        if (v == 0UL)
+        {
+            return u;
+        }
+
+        int shift = XMath.TrailingZeroCount(u | v);
+        u >>= XMath.TrailingZeroCount(u);
+
+        do
+        {
+            v >>= XMath.TrailingZeroCount(v);
+            if (u > v)
+            {
+				(u, v) = (v, u);
+            }
+
+            v -= u;
+        }
+        while (v != 0UL);
+
+        return u << shift;
+    }
+}

PerfectNumbers.Core/PrimesGenerator.cs

@@ -0,0 +1,106 @@
+namespace PerfectNumbers.Core;
+
+public static class PrimesGenerator
+{
+	private static ulong[]? _primes;
+	private static ulong[]? _primesPow2;	
+	// private static UInt128[]? _primes128;
+	// private static UInt128[]? _primesPow2_128;	
+
+	public static readonly ulong[] SmallPrimes = BuildSmallPrimesInternal();
+	public static readonly ulong[] SmallPrimesPow2 = BuildSmallPrimesPow2Internal();
+	// public static readonly UInt128[] SmallPrimes128 = BuildSmallPrimes128Internal();
+	// public static readonly UInt128[] SmallPrimesPow2_128 = BuildSmallPrimesPow2_128Internal();
+	
+
+	private static ulong[] BuildSmallPrimesInternal()
+	{
+		if (_primes == null)
+		{
+			BuildSmallPrimes();
+		}
+
+		return _primes!;
+	}
+
+	private static ulong[] BuildSmallPrimesPow2Internal()
+	{
+		if (_primesPow2 == null)
+		{
+			BuildSmallPrimes();
+		}
+
+		return _primesPow2!;
+	}
+
+	// private static UInt128[] BuildSmallPrimes128Internal()
+	// {
+	// 	if (_primes == null)
+	// 	{
+	// 		BuildSmallPrimes();
+	// 	}
+
+	// 	return _primes128!;
+	// }
+
+	// private static UInt128[] BuildSmallPrimesPow2_128Internal()
+	// {
+	// 	if (_primesPow2 == null)
+	// 	{
+	// 		BuildSmallPrimes();
+	// 	}
+
+	// 	return _primesPow2_128!;
+	// }
+
+	private static void BuildSmallPrimes()
+	{
+               ulong limit = PerfectNumberConstants.PrimesLimit;
+               var primes = new ulong[limit];
+               primes[0] = 2UL;
+               var primesPow2 = new ulong[limit];
+               primesPow2[0] = 4UL;
+
+               int j;
+               ulong i = 3UL, p, sqrt;
+               bool isPrime;
+               int primeIndex = 1;
+
+               for (; i <= limit; i += 2UL)
+               {
+                       isPrime = true;
+                       sqrt = (ulong)Math.Sqrt(i);
+                       for (j = 0; j < primeIndex; j++)
+                       {
+                               p = primes[j];
+                               if (p > sqrt)
+                               {
+                                       break;
+                               }
+
+                               if (i % p == 0UL)
+                               {
+                                       isPrime = false;
+                                       break;
+                               }
+                       }
+
+                       if (isPrime)
+                       {
+                               primes[primeIndex] = i;
+                               primesPow2[primeIndex] = i * i;
+                               primeIndex++;
+                       }
+               }
+
+               Array.Resize(ref primes, primeIndex);
+               Array.Resize(ref primesPow2, primeIndex);
+
+               _primes = primes;
+               _primesPow2 = primesPow2;
+               // _primes128 = Array.ConvertAll(_primes, p => (UInt128)p);
+               // _primesPow2_128 = Array.ConvertAll(_primesPow2, p => (UInt128)p);
+
+               return;
+       }
+}

PerfectNumbers.Core/RationalHelper.cs

@@ -0,0 +1,35 @@
+using System.Numerics;
+using PeterO.Numbers;
+
+namespace PerfectNumbers.Core;
+
+public static class RationalHelper
+{
+    public static ERational FromDoubleRounded(double value, int scale = 100)
+    {
+        long numerator = checked((long)Math.Round(value * scale));
+        long denominator = scale;
+        long gcd = (long)BigInteger.GreatestCommonDivisor(numerator, denominator);
+        return ERational.Create(EInteger.FromInt64(numerator / gcd), EInteger.FromInt64(denominator / gcd));
+    }
+
+    public static EInteger ToEInteger(BigInteger value)
+    {
+        return EInteger.FromBytes(value.ToByteArray(), true);
+    }
+
+    public static ulong FloorToUInt64(ERational value)
+    {
+        EInteger ei = value.ToEInteger();
+        return ei.ToUInt64Checked();
+    }
+
+    public static ulong CeilingToUInt64(ERational value)
+    {
+        EInteger ei = value.ToEInteger();
+        if (value.CompareTo(ERational.Create(ei, EInteger.One)) > 0)
+            ei += EInteger.One;
+        return ei.ToUInt64Checked();
+    }
+}
+

PerfectNumbers.Core/RationalNumbers.cs

@@ -0,0 +1,8 @@
+using PeterO.Numbers;
+
+namespace PerfectNumbers.Core;
+
+public static class RationalNumbers
+{
+    public static ERational Two { get; } = ERational.FromInt32(2);
+}

PerfectNumbers.Core/RleBlacklist.cs

@@ -0,0 +1,361 @@
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+public static class RleBlacklist
+{
+    private static readonly object Sync = new();
+    private static volatile bool _loaded;
+    private static HashSet<string>? _patterns;
+
+    public static void Load(string path)
+    {
+        if (string.IsNullOrWhiteSpace(path))
+        {
+            return;
+        }
+
+        if (_loaded)
+        {
+            return;
+        }
+
+        lock (Sync)
+        {
+            if (_loaded)
+            {
+                return;
+            }
+
+            var set = new HashSet<string>(StringComparer.Ordinal);
+            foreach (var line in File.ReadLines(path))
+            {
+                if (string.IsNullOrWhiteSpace(line))
+                {
+                    continue;
+                }
+
+                var span = line.AsSpan().Trim();
+                if (span.Length == 0 || span[0] == '#')
+                {
+                    continue;
+                }
+
+                // CSV support: Results/rle_blacklist_composite_only_end7_base100k_validated_primes1M.csv
+                // Header: rle_pattern,label,count_in_composites_end7_le_100k,appears_in_primes_le_1M
+                // Row example: "[1, 1, 3]",composite,12,False
+                // We only include rows with label==composite and appears_in_primes_le_1M==False.
+                bool handledCsv = false;
+                int commaIdx = span.IndexOf(',');
+                if (commaIdx > 0)
+                {
+                    ReadOnlySpan<char> firstField;
+                    int startNext;
+                    if (span[0] == '"')
+                    {
+                        int endQuote = span.Slice(1).IndexOf('"');
+                        if (endQuote > 0)
+                        {
+                            endQuote += 1; // adjust to original span
+                            firstField = span.Slice(1, endQuote - 1);
+                            // move to char after closing quote, expect comma
+                            int after = endQuote + 1;
+                            if (after < span.Length && span[after] == ',')
+                            {
+                                startNext = after + 1;
+                            }
+                            else
+                            {
+                                startNext = commaIdx + 1;
+                            }
+                        }
+                        else
+                        {
+                            firstField = span.Slice(0, commaIdx);
+                            startNext = commaIdx + 1;
+                        }
+                    }
+                    else
+                    {
+                        firstField = span.Slice(0, commaIdx);
+                        startNext = commaIdx + 1;
+                    }
+
+                    // Parse remaining CSV fields (label, count, appears_in_primes)
+                    // label
+                    int comma2 = span.Slice(startNext).IndexOf(',');
+                    if (comma2 > 0)
+                    {
+                        ReadOnlySpan<char> label = span.Slice(startNext, comma2).Trim();
+                        int startCount = startNext + comma2 + 1;
+                        int comma3 = span.Slice(startCount).IndexOf(',');
+                        if (comma3 > 0)
+                        {
+                            // Read appears_in_primes flag
+                            int startFlag = startCount + comma3 + 1;
+                            ReadOnlySpan<char> appears = span.Slice(startFlag).Trim().Trim('"');
+                            bool composite = label.Equals("composite", StringComparison.OrdinalIgnoreCase);
+                            bool appearsInPrimes = appears.Equals("True", StringComparison.OrdinalIgnoreCase);
+                            if (composite && !appearsInPrimes)
+                            {
+                                var key = NormalizeBracketList(firstField);
+                                if (key != null)
+                                {
+                                    set.Add(key);
+                                }
+                            }
+
+                            handledCsv = true;
+                        }
+                    }
+                }
+
+                if (handledCsv)
+                {
+                    continue;
+                }
+
+                // Accept patterns as space/comma/dash-separated integers; optionally prefixed with firstbit|.
+                var normalized = NormalizePattern(span);
+                if (normalized != null)
+                {
+                    set.Add(normalized);
+                }
+            }
+
+            _patterns = set;
+            _loaded = true;
+        }
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool IsLoaded() => _loaded && _patterns != null && _patterns.Count != 0;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static bool Matches(ulong p)
+    {
+        if (!IsLoaded())
+        {
+            return false;
+        }
+
+        string key = BuildRleKey(p);
+        return _patterns!.Contains(key);
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static string BuildRleKey(ulong value)
+    {
+        // Build RLE from MSB to LSB runs.
+        int bitLen = 64 - int.CreateChecked(ulong.LeadingZeroCount(value));
+        if (bitLen <= 0)
+        {
+            return "";
+        }
+
+        Span<int> runs = stackalloc int[64];
+        int runCount = 0;
+        int currentBit = (int)((value >> (bitLen - 1)) & 1UL);
+        int count = 1;
+        for (int i = bitLen - 2; i >= 0; i--)
+        {
+            int b = (int)((value >> i) & 1UL);
+            if (b == currentBit)
+            {
+                count++;
+            }
+            else
+            {
+                runs[runCount++] = count;
+                currentBit = b;
+                count = 1;
+            }
+        }
+
+        runs[runCount++] = count;
+
+        // Key format: first bit value; then run lengths joined by '-'
+        // Example: for binary 101100 -> first=1, runs: 1,1,2,2 -> "1|1-1-2-2"
+        // This preserves starting bit which some datasets may rely on.
+        var sb = StringBuilderPool.Rent();
+        try
+        {
+            _ = sb.Append(((value >> (bitLen - 1)) & 1UL) == 0 ? '0' : '1');
+            _ = sb.Append('|');
+            for (int i = 0; i < runCount; i++)
+            {
+                if (i != 0)
+                {
+                    _ = sb.Append('-');
+                }
+
+                _ = sb.Append(runs[i]);
+            }
+
+            return sb.ToString();
+        }
+        finally
+        {
+            StringBuilderPool.Return(sb);
+        }
+    }
+
+    private static string? NormalizePattern(ReadOnlySpan<char> span)
+    {
+        // Accept optional "firstbit|a-b-c" or just "a b c" / "a,b,c".
+        int pipe = span.IndexOf('|');
+        int start = 0;
+        char first = '1';
+        if (pipe >= 0)
+        {
+            if (pipe > 0 && (span[0] == '0' || span[0] == '1'))
+            {
+                first = span[0];
+            }
+
+            start = pipe + 1;
+        }
+
+        Span<int> tmp = stackalloc int[128];
+        int count = 0;
+        int i = start;
+        while (i <= span.Length)
+        {
+            // Read token until separator: space/comma/dash or end
+            int j = i;
+            while (j < span.Length)
+            {
+                char c = span[j];
+                if (c == ' ' || c == '\t' || c == ',' || c == ';' || c == '-')
+                {
+                    break;
+                }
+
+                j++;
+            }
+
+            if (j > i)
+            {
+                if (int.TryParse(span.Slice(i, j - i), out int v) && v > 0)
+                {
+                    tmp[count++] = v;
+                }
+                else
+                {
+                    return null;
+                }
+            }
+
+            if (j >= span.Length)
+            {
+                break;
+            }
+
+            i = j + 1;
+        }
+
+        if (count == 0)
+        {
+            return null;
+        }
+
+        // Build normalized key
+        var sb = StringBuilderPool.Rent();
+        try
+        {
+            _ = sb.Append(first);
+            _ = sb.Append('|');
+            for (int k = 0; k < count; k++)
+            {
+                if (k != 0)
+                {
+                    _ = sb.Append('-');
+                }
+
+                _ = sb.Append(tmp[k]);
+            }
+
+            return sb.ToString();
+        }
+        finally
+        {
+            StringBuilderPool.Return(sb);
+        }
+    }
+
+    private static string? NormalizeBracketList(ReadOnlySpan<char> span)
+    {
+        // Expect a JSON-like list: [1, 2, 3]
+        int l = span.IndexOf('[');
+        int r = span.LastIndexOf(']');
+        if (l < 0 || r <= l)
+        {
+            return null;
+        }
+
+        ReadOnlySpan<char> inner = span.Slice(l + 1, r - l - 1);
+        Span<int> tmp = stackalloc int[128];
+        int count = 0;
+        int i = 0;
+        while (i <= inner.Length)
+        {
+            while (i < inner.Length && (inner[i] == ' ' || inner[i] == '\t' || inner[i] == ','))
+            {
+                i++;
+            }
+
+            if (i >= inner.Length)
+            {
+                break;
+            }
+
+            int j = i;
+            while (j < inner.Length && inner[j] != ',' && inner[j] != ' ' && inner[j] != '\t')
+            {
+                j++;
+            }
+
+            if (j > i)
+            {
+                if (int.TryParse(inner.Slice(i, j - i), out int v) && v > 0)
+                {
+                    tmp[count++] = v;
+                }
+                else
+                {
+                    return null;
+                }
+            }
+
+            i = j + 1;
+        }
+
+        if (count == 0)
+        {
+            return null;
+        }
+
+        // Assume first bit is 1 (MSB is 1 for any positive number), so starting run is ones.
+        var sb = StringBuilderPool.Rent();
+        try
+        {
+            _ = sb.Append('1');
+            _ = sb.Append('|');
+            for (int k = 0; k < count; k++)
+            {
+                if (k != 0)
+                {
+                    _ = sb.Append('-');
+                }
+
+                _ = sb.Append(tmp[k]);
+            }
+
+            return sb.ToString();
+        }
+        finally
+        {
+            StringBuilderPool.Return(sb);
+        }
+    }
+}

PerfectNumbers.Core/StringBuilderPool.cs

@@ -0,0 +1,17 @@
+using System.Collections.Concurrent;
+using System.Text;
+
+namespace PerfectNumbers.Core;
+
+public static class StringBuilderPool
+{
+    private static readonly ConcurrentQueue<StringBuilder> _stringBuilderPool = new();
+
+    public static StringBuilder Rent() => _stringBuilderPool.TryDequeue(out var sb) ? sb : new();
+
+    public static void Return(StringBuilder sb)
+    {
+        _ = sb.Clear();
+        _stringBuilderPool.Enqueue(sb);
+    }
+}

PerfectNumbers.Core/TextFileWriter.cs

@@ -0,0 +1,77 @@
+namespace PerfectNumbers.Core;
+
+/// <summary>
+/// Provides thread-safe file writing for both text and binary data.
+/// </summary>
+public sealed class TextFileWriter : IDisposable
+{
+    private readonly FileStream _stream;
+    private StreamWriter? _textWriter;
+    private BinaryWriter? _binaryWriter;
+    private readonly object _lock = new();
+
+    public TextFileWriter(string path, FileMode mode)
+    {
+        var options = new FileStreamOptions
+        {
+            Access = FileAccess.Write,
+            Mode = mode,
+            Share = FileShare.Read
+        };
+        _stream = new FileStream(path, options);
+    }
+
+    public long Position
+    {
+        get { lock (_lock) { return _stream.Position; } }
+    }
+
+    private StreamWriter TextWriter => _textWriter ??= new StreamWriter(_stream, leaveOpen: true);
+    private BinaryWriter BinaryWriter => _binaryWriter ??= new BinaryWriter(_stream, System.Text.Encoding.UTF8, leaveOpen: true);
+
+    /// <summary>
+    /// Clears the file by truncating it to zero length.
+    /// </summary>
+    public void Clear()
+    {
+        lock (_lock)
+        {
+            _stream.SetLength(0);
+            _stream.Position = 0;
+            _textWriter?.Flush();
+            _binaryWriter?.Flush();
+        }
+    }
+
+    /// <summary>
+    /// Writes a line of text to the file.
+    /// </summary>
+    public void WriteLine(string line)
+    {
+        lock (_lock)
+        {
+            TextWriter.WriteLine(line);
+            TextWriter.Flush();
+        }
+    }
+
+    /// <summary>
+    /// Executes a binary writing action within the internal lock.
+    /// </summary>
+    public void WriteBinary(Action<BinaryWriter> action)
+    {
+        lock (_lock)
+        {
+            action(BinaryWriter);
+            BinaryWriter.Flush();
+        }
+    }
+
+    public void Dispose()
+    {
+        _binaryWriter?.Dispose();
+        _textWriter?.Dispose();
+        _stream.Dispose();
+    }
+}
+

PerfectNumbers.Core/UInt128Extensions.cs

@@ -0,0 +1,251 @@
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+public static class UInt128Extensions
+{
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static UInt128 BinaryGcd(this UInt128 u, UInt128 v)
+	{
+		UInt128 zero = UInt128.Zero;
+		if (u == zero)
+		{
+			return v;
+		}
+
+		if (v == zero)
+		{
+			return u;
+		}
+
+		int shift = CountTrailingZeros(u | v);
+		u >>= CountTrailingZeros(u);
+
+		do
+		{
+			v >>= CountTrailingZeros(v);
+			if (u > v)
+			{
+				(u, v) = (v, u);
+			}
+
+			v -= u;
+		}
+		while (v != zero);
+
+		return u << shift;
+	}
+
+	public static ulong CalculateOrder(this UInt128 q)
+	{
+		if (q <= UInt128Numbers.Two)
+		{
+			return 0UL;
+		}
+
+		UInt128 one = UInt128.One;
+		UInt128 phi = q - one;
+		if (phi > ulong.MaxValue)
+		{
+			throw new NotImplementedException("Such big values are not yet supported");
+		}
+
+		ulong order = (ulong)phi;
+		ulong[] smallPrimes = PrimesGenerator.SmallPrimes,
+				smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
+
+		int i = 0, primesLength = smallPrimes.Length;
+		UInt128 cycle = MersenneDivisorCycles.GetCycle(q);
+		ulong prime, temp;
+		for (; i < primesLength; i++)
+		{
+			if (smallPrimesPow2[i] > order)
+			{
+				break;
+			}
+
+			prime = smallPrimes[i];
+			while (order % prime == 0UL)
+			{
+				temp = order / prime;
+				if (temp.PowModWithCycle(q, cycle) == one)
+				{
+					order = temp;
+				}
+				else
+				{
+					break;
+				}
+			}
+		}
+
+		return order;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static int CountTrailingZeros(this UInt128 value)
+	{
+		ulong valuePart = (ulong)value;
+		if (valuePart != 0UL)
+		{
+			return BitOperations.TrailingZeroCount(valuePart);
+		}
+
+		return 64 + BitOperations.TrailingZeroCount((ulong)(value >> 64));
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static bool IsPrimeCandidate(this UInt128 n)
+	{
+		UInt128 p, zero = UInt128.Zero;
+
+		ulong[] smallPrimes = PrimesGenerator.SmallPrimes;
+		ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
+		ulong i, smallPrimesCount = (ulong)smallPrimes.Length;
+		for (i = 0UL; i < smallPrimesCount; i++)
+		{
+			p = smallPrimes[i];
+			if (smallPrimesPow2[i] > n)
+			{
+				break;
+			}
+
+			if (n % p == zero)
+			{
+				return n == p;
+			}
+		}
+
+		return true;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static void Mod10_8_5_3(this UInt128 value, out ulong mod10, out ulong mod8, out ulong mod5, out ulong mod3)
+	{
+		UInt128 zero = UInt128.Zero;
+		if (value == zero)
+		{
+			mod3 = 0UL;
+			mod5 = 0UL;
+			mod8 = 0UL;
+			mod10 = 0UL;
+			return;
+		}
+
+		UInt128 high = value >> 64;
+		ulong result = (ulong)value;
+		mod8 = result & 7UL;
+		// 2^64 ≡ 1 (mod 3)
+		ulong modRem = (result % 3UL) + ((ulong)high % 3UL);
+		mod3 = modRem >= 3UL ? modRem - 3UL : modRem;
+		// 2^64 ≡ 1 (mod 5)
+		modRem = (result % 5UL) + ((ulong)high % 5UL);
+		mod5 = modRem >= 5UL ? modRem - 5UL : modRem;
+
+		while (high != zero)
+		{
+			// 2^64 ≡ 6 (mod 10)
+			result = (result + (ulong)high * 6UL) % 10UL;
+			high >>= 64;
+		}
+
+		mod10 = result % 10UL;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static UInt128 ModPow(this UInt128 value, UInt128 exponent, UInt128 modulus)
+	{
+		UInt128 result, one = result = UInt128.One,
+				zero = UInt128.Zero,
+				baseValue = value % modulus;
+
+		while (exponent != zero)
+		{
+			if ((exponent & one) != zero)
+			{
+				result = MulMod(result, baseValue, modulus);
+			}
+
+			baseValue = MulMod(baseValue, baseValue, modulus);
+			exponent >>= 1;
+		}
+
+		return result;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong Mod10(this UInt128 value)
+	{
+		UInt128 zero = UInt128.Zero;
+		if (value == zero)
+			return 0UL;
+
+		ulong result = (ulong)value;
+		UInt128 high = value >> 64;
+
+		while (high != zero)
+		{
+			// 2^64 ≡ 6 (mod 10)
+			result = (result + (ulong)high * 6UL) % 10UL;
+			high >>= 64;
+		}
+
+		return result % 10UL;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong Mod8(this UInt128 value) => (ulong)value & 7UL;
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong Mod3(this UInt128 value)
+	{
+		// 2^64 ≡ 1 (mod 3)
+		ulong rem = ((ulong)value % 3UL) + ((ulong)(value >> 64) % 3UL);
+		return rem >= 3UL ? rem - 3UL : rem;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong Mod5(this UInt128 value)
+	{
+		// 2^64 ≡ 1 (mod 5)
+		ulong rem = ((ulong)value % 5UL) + ((ulong)(value >> 64) % 5UL);
+		return rem >= 5UL ? rem - 5UL : rem;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static UInt128 Mul64(this UInt128 a, UInt128 b)
+	{
+		ulong aLow = (ulong)a;
+		ulong bLow = (ulong)b;
+		return ((UInt128)(aLow * (ulong)(b >> 64) + aLow.MulHigh(bLow)) << 64) | (aLow * bLow);
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static UInt128 MulMod(this UInt128 a, UInt128 b, UInt128 modulus)
+	{
+		UInt128 zero, result = zero = UInt128.Zero, one = UInt128.One;
+
+		while (b != zero)
+		{
+			if ((b & one) != zero)
+			{
+				result += a;
+				if (result >= modulus)
+				{
+					result -= modulus;
+				}
+			}
+
+			a <<= 1;
+			if (a >= modulus)
+			{
+				a -= modulus;
+			}
+			
+			b >>= 1;
+		}
+
+		return result;
+	}
+}

PerfectNumbers.Core/UInt128Numbers.cs

@@ -0,0 +1,10 @@
+namespace PerfectNumbers.Core
+{
+	public static class UInt128Numbers
+	{
+		public static readonly UInt128 Two = 2UL;
+		public static readonly UInt128 Four = 4UL;
+		public static readonly UInt128 SixtyFour = 64UL;
+		public static readonly UInt128 OneHundredTwentySeven = 127UL;
+    }
+}

PerfectNumbers.Core/ULongExtensions.cs

@@ -0,0 +1,423 @@
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace PerfectNumbers.Core;
+
+public static class ULongExtensions
+{
+	public static ulong CalculateOrder(this ulong q)
+	{
+		if (q <= 2UL)
+		{
+			return 0UL;
+		}
+
+		ulong order = q - 1UL;
+		ulong[] smallPrimes = PrimesGenerator.SmallPrimes;
+		ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
+
+		int i = 0, primesLength = smallPrimes.Length;
+		ulong prime, temp;
+		UInt128 q128 = q,
+				cycle = MersenneDivisorCycles.GetCycle(q128);
+
+		for (; i < primesLength; i++)
+		{
+			prime = smallPrimes[i];
+			if (smallPrimesPow2[i] > order)
+			{
+				break;
+			}
+
+			while (order % prime == 0UL)
+			{
+				temp = order / prime;
+				if (temp.PowModWithCycle(q128, cycle) == UInt128.One)
+				{
+					order = temp;
+				}
+				else
+				{
+					break;
+				}
+			}
+		}
+
+		return order;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong FastDiv64(this ulong value, ulong divisor, ulong mul)
+	{
+		ulong q = value.MulHigh(mul);
+		if ((UInt128)value - q.Mul64(divisor) >= divisor)
+		{
+			q++;
+		}
+
+		return q;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static bool IsPrimeCandidate(this ulong n)
+	{
+		int i = 0;
+		ulong p;
+		ulong[] smallPrimes = PrimesGenerator.SmallPrimes;
+		ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
+		int len = smallPrimes.Length;
+		for (; i < len; i++)
+		{
+			if (smallPrimesPow2[i] > n)
+			{
+				break;
+			}
+
+			p = smallPrimes[i];
+			if ((n % p) == 0UL)
+			{
+				return n == p;
+			}
+		}
+
+		return true;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong Mod10(ulong value)
+	{
+		ulong quotient = (ulong)(((UInt128)value * 0xCCCCCCCCCCCCCCCDUL) >> 67);
+		return value - quotient * 10UL;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static UInt128 Mul64(this ulong a, ulong b) => ((UInt128)a.MulHigh(b) << 64) | (UInt128)(a * b);
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong MulHigh(this ulong x, ulong y)
+	{
+		ulong xLow = (uint)x;
+		ulong xHigh = x >> 32;
+		ulong yLow = (uint)y;
+		ulong yHigh = y >> 32;
+
+		ulong w1 = xLow * yHigh;
+		ulong w2 = xHigh * yLow;
+		ulong carry =
+			(((xLow * yLow) >> 32) +
+			(uint)w1 +
+			(uint)w2) >> 32;
+
+		return (xHigh * yHigh) + (w1 >> 32) + (w2 >> 32) + carry;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong ModPow64(this ulong value, ulong exponent, ulong modulus)
+	{
+		ulong result = 1UL;
+		ulong baseValue = value % modulus;
+
+		while (exponent != 0UL)
+		{
+			if ((exponent & 1UL) != 0UL)
+			{
+				result = MulMod64(result, baseValue, modulus);
+			}
+
+			baseValue = MulMod64(baseValue, baseValue, modulus);
+			exponent >>= 1;
+		}
+
+		return result;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus);
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static UInt128 PowMod(this ulong exponent, UInt128 modulus)
+	{
+		// Return 1 because 2^0 = 1
+		if (exponent == 0UL)
+			return UInt128.One;
+
+		// Any number mod 1 is 0
+		if (modulus == UInt128.One)
+			return UInt128.Zero;
+
+		UInt128 result;
+		// For small exponents, do classic method
+		if (exponent < 64 || modulus < 4)
+		{
+			result = UInt128.One;
+			for (ulong i = 0; i < exponent; i++)
+			{
+				result <<= 1;
+				if (result >= modulus)
+					result -= modulus;
+			}
+
+			return result;
+		}
+
+		// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
+		if ((modulus & (modulus - 1)) == 0)
+		{
+			result = UInt128.One << (int)(exponent & 127);
+			return result & (modulus - 1);
+		}
+
+		// For general modulus, use classic left-to-right binary method for 2^exponent mod modulus
+		UInt128 resultBin = UInt128.One;
+		for (ulong i = 0; i < exponent; i++)
+		{
+			resultBin <<= 1;
+			if (resultBin >= modulus)
+				resultBin -= modulus;
+		}
+
+		return resultBin;
+	}
+
+	/// <summary>
+	/// Computes 2^exponent mod modulus using a known cycle length.
+	/// </summary>
+	public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, ulong cycleLength)
+	{
+		// Return 1 because 2^0 = 1
+		if (exponent == 0UL)
+			return UInt128.One;
+
+		// Any number mod 1 is 0
+		if (modulus == UInt128.One)
+			return UInt128.Zero;
+
+		// For small exponents, do classic method
+		UInt128 result;
+		if (exponent < 64 || modulus < 4)
+		{
+			result = UInt128.One;
+			for (ulong i = 0; i < exponent; i++)
+			{
+				result <<= 1;
+				if (result >= modulus)
+					result -= modulus;
+			}
+
+			return result;
+		}
+
+		// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
+		if ((modulus & (modulus - 1)) == 0)
+		{
+			result = UInt128.One << (int)(exponent & 127);
+			return result & (modulus - 1);
+		}
+
+		// For general modulus, use classic left-to-right binary method for 2^exponent mod modulus
+		result = UInt128.One;
+		UInt128 rotations = exponent % cycleLength;
+		for (ulong i = 0; i < rotations; i++)
+		{
+			result <<= 1;
+			if (result >= modulus)
+				result -= modulus;
+		}
+
+		return result;
+	}
+
+	/// <summary>
+	/// Computes 2^exponent mod modulus using a known cycle length.
+	/// </summary>
+	public static UInt128 PowModWithCycle(this ulong exponent, UInt128 modulus, UInt128 cycleLength)
+	{
+		// Return 1 because 2^0 = 1
+		if (exponent == UInt128.Zero)
+			return UInt128.One;
+			
+		// Any number mod 1 is 0
+		if (modulus == UInt128.One)
+			return UInt128.Zero;
+
+		// For small exponents, do classic method
+		UInt128 result;
+		if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four)
+		{
+			result = UInt128.One;
+			for (ulong i = 0; i < exponent; i++)
+			{
+				result <<= 1;
+				if (result >= modulus)
+					result -= modulus;
+			}
+
+			return result;
+		}
+
+		// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
+		if ((modulus & (modulus - UInt128.One)) == UInt128.Zero)
+		{
+			result = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
+			return result & (modulus - 1);
+		}
+
+		// For general modulus, use classic left-to-right binary method for 2^exponent mod modulus
+		result = UInt128.One;
+		UInt128 rotations = exponent % cycleLength;
+		for (UInt128 i = 0; i < rotations; i++)
+		{
+			result <<= 1;
+			if (result >= modulus)
+				result -= modulus;
+		}
+
+		return result;
+	}
+
+	/// <summary>
+	/// Computes 2^exponent mod modulus using a known cycle length.
+	/// </summary>
+	public static UInt128 PowModWithCycle(this UInt128 exponent, UInt128 modulus, UInt128 cycleLength)
+	{
+		// Return 1 because 2^0 = 1
+		if (exponent == UInt128.Zero)
+			return UInt128.One;
+
+		// Any number mod 1 is 0
+		if (modulus == UInt128.One)
+			return UInt128.Zero;
+
+		UInt128 result;
+
+		// For small exponents, do classic method
+		if (exponent < UInt128Numbers.SixtyFour || modulus < UInt128Numbers.Four)
+		{
+			result = UInt128.One;
+			for (ulong i = 0; i < exponent; i++)
+			{
+				result <<= 1;
+				if (result >= modulus)
+					result -= modulus;
+			}
+
+			return result;
+		}
+
+		// Special case: if modulus is a power of two, use bitmasking for efficiency and correctness
+		if ((modulus & (modulus - UInt128.One)) == UInt128.Zero)
+		{
+			result = UInt128.One << (int)(exponent & UInt128Numbers.OneHundredTwentySeven);
+			return result & (modulus - 1);
+		}
+
+		// For general modulus, use classic left-to-right binary method for 2^exponent mod modulus
+		result = UInt128.One;
+		UInt128 rotations = exponent % cycleLength;
+		for (UInt128 i = 0; i < rotations; i++)
+		{
+			result <<= 1;
+			if (result >= modulus)
+				result -= modulus;
+		}
+
+		return result;
+	}
+
+	/// <summary>
+	/// Computes 2^exponent mod modulus using iterative CRT composition from mod 10 up to modulus.
+	/// Only for modulus >= 10 and reasonable size.
+	/// </summary>
+	public static UInt128 PowModCrt(this ulong exponent, UInt128 modulus, MersenneDivisorCycles cycles)
+	{
+		if (modulus < 10)
+			return PowMod(exponent, modulus); // fallback to classic
+
+		// Use cycle length 4 for mod 10
+		UInt128 result = PowModWithCycle(exponent, 10, 4);
+
+		UInt128 currentModulus = 10;
+		for (UInt128 m = 11; m <= modulus; m++)
+		{
+			UInt128 cycle = MersenneDivisorCycles.GetCycle(m);
+			UInt128 remM = cycle > 0
+				? PowModWithCycle(exponent, m, cycle)
+				: PowMod(exponent, m);
+
+			// Solve x ≡ result mod currentModulus
+			//      x ≡ remM   mod m
+			// Find x mod (currentModulus * m)
+			// Since currentModulus and m are coprime, use CRT:
+			// x = result + currentModulus * t, where t ≡ (remM - result) * inv(currentModulus, m) mod m
+
+			UInt128 inv = ModInverse(currentModulus, m);
+			UInt128 t = ((remM + m - (result % m)) * inv) % m;
+			result += currentModulus * t;
+			currentModulus *= m;
+
+			if (currentModulus >= modulus)
+				break;
+		}
+
+		return result % modulus;
+	}
+
+	// Helper: modular inverse (extended Euclidean algorithm)
+	private static UInt128 ModInverse(UInt128 a, UInt128 m)
+	{
+		UInt128 m0 = m, t, q;
+		UInt128 x0 = 0, x1 = 1;
+		if (m == 1) return 0;
+		while (a > 1)
+		{
+			q = a / m;
+			t = m;
+			m = a % m; a = t;
+			t = x0;
+			x0 = x1 - q * x0;
+			x1 = t;
+		}
+
+		if (x1 < 0) x1 += m0;
+
+		return x1;
+	}
+
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	public static bool SharesFactorWithExponentMinusOne(this ulong exponent)
+	{
+		ulong value = exponent - 1UL;
+		value >>= BitOperations.TrailingZeroCount(value);
+		ulong[] smallPrimes = PrimesGenerator.SmallPrimes;
+		ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2;
+		int i = 1, smallPrimesLength = smallPrimes.Length;
+		ulong prime;
+
+		for (; i < smallPrimesLength && smallPrimesPow2[i] <= value; i++)
+		{
+			prime = smallPrimes[i];
+			if (value % prime != 0UL)
+			{
+				continue;
+			}
+
+			if (exponent % prime.CalculateOrder() == 0UL)
+			{
+				return true;
+			}
+
+			do
+			{
+				value /= prime;
+			}
+			while (value % prime == 0UL);
+		}
+
+		if (value > 1UL && exponent % value.CalculateOrder() == 0UL)
+		{
+			return true;
+		}
+
+		return false;
+	}
+}