Various improvements to make sure scratch buffers don't override each others memory

Add GPU lease warm-up and shared limiter

EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs

@@ -1,6 +1,3 @@
-using System;
-using System.Collections.Generic;
-using System.Threading;
 using BenchmarkDotNet.Attributes;
 using BenchmarkDotNet.Jobs;
 using Open.Numeric.Primes;
@@ -12,7 +9,6 @@ namespace EvenPerfectBitScanner.Benchmarks;
 [MemoryDiagnoser]
 public class PrimeTesterBenchmarks
 {
-    private static readonly HeuristicCombinedPrimeTester _heuristicCombinedTester = new();
     private static readonly HeuristicPrimeTester _heuristicTester = new();
     private ulong[] _candidates = [];
 
@@ -34,7 +30,7 @@ public class PrimeTesterBenchmarks
         _ = Prime.Numbers.IsPrime(_candidates[0]);
     }
 
-    [Benchmark(Baseline = true)]
+    [Benchmark]
     public int HeuristicCpu()
     {
         HeuristicPrimeTester tester = _heuristicTester;
@@ -104,18 +100,15 @@ public class PrimeTesterBenchmarks
         return primeCount;
     }
 
-	private static readonly PrimeTester _nonHeuristicTester = new();
-
     [Benchmark]
     public int NonHeuristicGpu()
 	{
-		PrimeTester tester = _nonHeuristicTester;
         ulong[] values = _candidates;
         int primeCount = 0;
 
         for (int i = 0; i < values.Length; i++)
         {
-            if (tester.IsPrimeGpu(values[i]))
+            if (PrimeTester.IsPrimeGpu(values[i]))
             {
                 primeCount++;
             }
@@ -124,7 +117,7 @@ public class PrimeTesterBenchmarks
         return primeCount;
     }
 
-    [Benchmark]
+    [Benchmark(Baseline = true)]
     public int NonHeuristicCpu()
     {
         ulong[] values = _candidates;

EvenPerfectBitScanner/Program.cs

@@ -12,7 +12,7 @@ namespace EvenPerfectBitScanner;
 
 internal static class Program
 {
-	private static ThreadLocal<PrimeTester> PrimeTesters = null!;
+	// private static ThreadLocal<PrimeTester> PrimeTesters = null!;
 	private static ThreadLocal<MersenneNumberTester> MersenneTesters = null!;
 	private const ulong InitialP = PerfectNumberConstants.BiggestKnownEvenPerfectP;
 	private static MersenneNumberDivisorGpuTester? _divisorTester;
@@ -58,14 +58,15 @@ internal static class Program
 			bool startPrimeProvided = _cliArguments.StartPrimeProvided;
 			int threadCount = Math.Max(1, _cliArguments.ThreadCount);
 			int gpuPrimeBatch = Math.Max(1, _cliArguments.GpuPrimeBatch);
+			int gpuPrimeThreads = Math.Max(1, _cliArguments.GpuPrimeThreads);
 			UnboundedTaskScheduler.ConfigureThreadCount(threadCount);
 			PrimeTester.GpuBatchSize = gpuPrimeBatch;
+			GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads);
 			GpuContextPool.WarmUpPool(threadCount);
-			PrimeTester.WarmUpGpuKernels(threadCount);
+			PrimeTester.WarmUpGpuKernels(gpuPrimeThreads);
 			Console.WriteLine("Starting up threads...");
 			_ = UnboundedTaskScheduler.Instance;
 			int blockSize = Math.Max(1, _cliArguments.BlockSize);
-			int gpuPrimeThreads = Math.Max(1, _cliArguments.GpuPrimeThreads);
 			GpuKernelType kernelType = _cliArguments.KernelType;
 			bool useBitTransform = _cliArguments.UseBitTransform;
 			bool useOrder = _cliArguments.UseOrder;
@@ -117,9 +118,6 @@ internal static class Program
 				return;
 			}
 
-			// Apply GPU prime sieve runtime configuration
-			GpuPrimeWorkLimiter.SetLimit(gpuPrimeThreads);
-			PrimeTester.GpuBatchSize = gpuPrimeBatch;
 
 			MersenneDivisorCycles mersenneDivisorCycles = new();
 			if (!File.Exists(cyclesPath) || continueCyclesGeneration)
@@ -178,7 +176,7 @@ internal static class Program
 			blockSize = useByDivisor ? 1 : blockSize;
 
 			CandidatesCalculator.Initialize(currentP);
-			PrimeTesters = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
+			// PrimeTesters = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
 			// Note: --primes-device controls default device for library kernels; p primality remains CPU here.
 			// Initialize per-thread p residue tracker (Identity model) at currentP
 			if (!useDivisor && !_cliArguments.UseByDivisor)
@@ -628,7 +626,7 @@ internal static class Program
 		{
 			if (!(_runPrimesOnCpu
 					? PrimeTester.IsPrime(p)
-					: PrimeTesters.Value!.IsPrimeGpu(p)))
+					: PrimeTester.IsPrimeGpu(p)))
 			{
 				_lastCompositeP = true;
 				return false;

PerfectNumbers.Core.Tests/Gpu/GpuPrimeLimiterTests.cs

@@ -18,19 +18,19 @@ public class GpuPrimeLimiterTests
         var startedSecond = new TaskCompletionSource<bool>();
         var sw = new Stopwatch();
 
-        var first = GpuPrimeWorkLimiter.Acquire();
+        GpuPrimeWorkLimiter.Acquire();
         sw.Start();
 
         var secondTask = Task.Run(() =>
         {
-            var second = GpuPrimeWorkLimiter.Acquire();
+            GpuPrimeWorkLimiter.Acquire();
             try
             {
                 startedSecond.TrySetResult(true);
             }
             finally
             {
-                second.Dispose();
+				GpuPrimeWorkLimiter.Release();
             }
         });
 
@@ -43,7 +43,7 @@ public class GpuPrimeLimiterTests
         finally
         {
             // Release first after the initial wait so the second task can proceed.
-            first.Dispose();
+			GpuPrimeWorkLimiter.Release();
         }
 
         await secondTask.WaitAsync(TimeSpan.FromSeconds(1));

PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs

@@ -16,7 +16,7 @@ public class PrimeTesterGpuTests
             ulong[] primes = [31UL, 61UL, 89UL, 107UL, 127UL, 521UL];
             foreach (ulong prime in primes)
             {
-                tester.IsPrimeGpu(prime).Should().BeTrue();
+				PrimeTester.IsPrimeGpu(prime).Should().BeTrue();
             }
         }
         finally
@@ -34,7 +34,7 @@ public class PrimeTesterGpuTests
             ulong[] composites = [33UL, 39UL, 51UL, 77UL, 91UL];
             foreach (ulong composite in composites)
             {
-                tester.IsPrimeGpu(composite).Should().BeFalse();
+				PrimeTester.IsPrimeGpu(composite).Should().BeFalse();
             }
         }
         finally

PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs

@@ -39,7 +39,7 @@ public class PrimeTesterIsPrimeTests
     {
         var tester = new PrimeTester();
 
-        tester.IsPrimeGpu(n).Should().Be(expected);
+		PrimeTester.IsPrimeGpu(n).Should().Be(expected);
     }
 
     [Fact]
@@ -64,7 +64,7 @@ public class PrimeTesterIsPrimeTests
         var testers = new ThreadLocal<PrimeTester>(() => new PrimeTester(), trackAllValues: true);
         Parallel.For(0, primes.Length, i =>
         {
-            results[i] = testers.Value!.IsPrimeGpu(primes[i]);
+            results[i] = PrimeTester.IsPrimeGpu(primes[i]);
         });
 
         results.Should().AllBeEquivalentTo(true);
@@ -73,7 +73,7 @@ public class PrimeTesterIsPrimeTests
         {
             foreach (ulong prime in primes)
             {
-                testers.Value!.IsPrimeGpu(prime).Should().BeTrue();
+				PrimeTester.IsPrimeGpu(prime).Should().BeTrue();
             }
         });
     }

PerfectNumbers.Core/Gpu/DivisorGeneratorGpuCache.cs

@@ -1,142 +0,0 @@
-using System;
-using System.Collections.Concurrent;
-using System.Runtime.CompilerServices;
-using ILGPU;
-using ILGPU.Runtime;
-
-namespace PerfectNumbers.Core.Gpu;
-
-internal static class DivisorGeneratorGpuCache
-{
-        private sealed class AcceleratorReferenceComparer : IEqualityComparer<Accelerator>
-        {
-                public static AcceleratorReferenceComparer Instance { get; } = new();
-
-                public bool Equals(Accelerator? x, Accelerator? y)
-                {
-                        return ReferenceEquals(x, y);
-                }
-
-                public int GetHashCode(Accelerator obj)
-                {
-                        return RuntimeHelpers.GetHashCode(obj);
-                }
-        }
-
-        private sealed class DivisorGeneratorGpuTables : IDisposable
-        {
-                internal DivisorGeneratorGpuTables(Accelerator accelerator)
-                {
-                        SmallPrimes = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimes);
-                        SmallPrimesPow2 = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2);
-                        SmallPrimesLastOne = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastOne);
-                        SmallPrimesPow2LastOne = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastOne);
-                        SmallPrimesLastSeven = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastSeven);
-                        SmallPrimesPow2LastSeven = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastSeven);
-                        SmallPrimesLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastThree);
-                        SmallPrimesPow2LastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastThree);
-                        SmallPrimesLastNine = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastNine);
-                        SmallPrimesPow2LastNine = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastNine);
-                        SmallPrimesLastOneWithoutLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastOneWithoutLastThree);
-                        SmallPrimesPow2LastOneWithoutLastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree);
-                        SmallPrimesLastSevenWithoutLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastSevenWithoutLastThree);
-                        SmallPrimesPow2LastSevenWithoutLastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree);
-                        SmallPrimesLastNineWithoutLastThree = AllocateAndCopyUInt(accelerator, DivisorGenerator.SmallPrimesLastNineWithoutLastThree);
-                        SmallPrimesPow2LastNineWithoutLastThree = AllocateAndCopyULong(accelerator, DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree);
-                }
-
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimes { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2 { get; }
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastOne { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastOne { get; }
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastSeven { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastSeven { get; }
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastThree { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastThree { get; }
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastNine { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastNine { get; }
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastOneWithoutLastThree { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastOneWithoutLastThree { get; }
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastSevenWithoutLastThree { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastSevenWithoutLastThree { get; }
-                internal MemoryBuffer1D<uint, Stride1D.Dense> SmallPrimesLastNineWithoutLastThree { get; }
-                internal MemoryBuffer1D<ulong, Stride1D.Dense> SmallPrimesPow2LastNineWithoutLastThree { get; }
-
-                public void Dispose()
-                {
-                        SmallPrimes.Dispose();
-                        SmallPrimesPow2.Dispose();
-                        SmallPrimesLastOne.Dispose();
-                        SmallPrimesPow2LastOne.Dispose();
-                        SmallPrimesLastSeven.Dispose();
-                        SmallPrimesPow2LastSeven.Dispose();
-                        SmallPrimesLastThree.Dispose();
-                        SmallPrimesPow2LastThree.Dispose();
-                        SmallPrimesLastNine.Dispose();
-                        SmallPrimesPow2LastNine.Dispose();
-                        SmallPrimesLastOneWithoutLastThree.Dispose();
-                        SmallPrimesPow2LastOneWithoutLastThree.Dispose();
-                        SmallPrimesLastSevenWithoutLastThree.Dispose();
-                        SmallPrimesPow2LastSevenWithoutLastThree.Dispose();
-                        SmallPrimesLastNineWithoutLastThree.Dispose();
-                        SmallPrimesPow2LastNineWithoutLastThree.Dispose();
-                }
-
-                private static MemoryBuffer1D<uint, Stride1D.Dense> AllocateAndCopyUInt(Accelerator accelerator, uint[] source)
-                {
-                        var buffer = accelerator.Allocate1D<uint>(source.Length);
-                        if (source.Length > 0)
-                        {
-                                buffer.View.CopyFromCPU(source);
-                        }
-
-                        return buffer;
-                }
-
-                private static MemoryBuffer1D<ulong, Stride1D.Dense> AllocateAndCopyULong(Accelerator accelerator, ulong[] source)
-                {
-                        var buffer = accelerator.Allocate1D<ulong>(source.Length);
-                        if (source.Length > 0)
-                        {
-                                buffer.View.CopyFromCPU(source);
-                        }
-
-                        return buffer;
-                }
-        }
-
-        private static readonly ConcurrentDictionary<Accelerator, DivisorGeneratorGpuTables> Cache = new(AcceleratorReferenceComparer.Instance);
-
-        public static void WarmUpByDivisorTables(Accelerator accelerator)
-        {
-                if (accelerator is null || accelerator.AcceleratorType == AcceleratorType.CPU)
-                {
-                        return;
-                }
-
-                Cache.GetOrAdd(accelerator, acc => new DivisorGeneratorGpuTables(acc));
-        }
-
-        public static void Clear(Accelerator accelerator)
-        {
-                if (accelerator is null)
-                {
-                        return;
-                }
-
-                if (Cache.TryRemove(accelerator, out var tables))
-                {
-                        tables.Dispose();
-                }
-        }
-
-        public static void DisposeAll()
-        {
-                foreach (var entry in Cache)
-                {
-                        entry.Value.Dispose();
-                }
-
-                Cache.Clear();
-        }
-}

PerfectNumbers.Core/Gpu/GpuContextPool.cs

@@ -1,4 +1,3 @@
-using System;
 using System.Collections.Concurrent;
 using ILGPU;
 using ILGPU.Runtime;
@@ -12,14 +11,15 @@ public static class GpuContextPool
 	// Default device preference for generic GPU kernels (prime scans, NTT, etc.)
 	internal sealed class PooledContext
 	{
-		public Context Context { get; }
-		public Accelerator Accelerator { get; }
+		public readonly Context Context;
+		public readonly Accelerator Accelerator;
 
-		public PooledContext()
+		public PooledContext(KernelContainer kernels)
 		{
-			Context = Context.CreateDefault();
-			Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context);
-			DivisorGeneratorGpuCache.WarmUpByDivisorTables(Accelerator);
+			Context = SharedGpuContext.Context;
+			Accelerator = SharedGpuContext.Accelerator;
+			GpuStaticTableInitializer.EnsureStaticTables(kernels, Accelerator);
+			kernels.Dispose();
 			// 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
@@ -33,11 +33,11 @@ public static class GpuContextPool
 			// Release PrimeTester GPU state for this accelerator (kernel/device primes)
 			// Clear any per-accelerator cached resources to avoid releasing
 			// after the accelerator is destroyed.
-			DivisorGeneratorGpuCache.Clear(Accelerator);
 			PrimeTester.ClearGpuCaches(Accelerator);
 
-			Accelerator.Dispose();
-			Context.Dispose();
+			// These resources are shared between GPU leases
+			//  Accelerator.Dispose();
+			// 	Context.Dispose();
 		}
 	}
 
@@ -53,9 +53,10 @@ public static class GpuContextPool
 		}
 
 		// Create a new accelerator when the pool does not have one available.
-		return new GpuContextLease(new PooledContext());
+		return new GpuContextLease(new PooledContext(GpuKernelPool.GetKernels()));
 	}
 
+
 	public static void WarmUpPool(int threadCount)
 	{
 		if (threadCount <= WarmedGpuContextCount)
@@ -64,15 +65,16 @@ public static class GpuContextPool
 		}
 
 		int toCreate = threadCount - WarmedGpuContextCount;
+		Console.WriteLine($"Warming up GPU context pool to {toCreate} contexts...");
 		var contexts = new PooledContext[toCreate];
+		var pool = GpuPool;
+		PooledContext context;
 		for (int i = 0; i < toCreate; i++)
 		{
-			contexts[i] = new PooledContext();
-		}
-
-		for (int i = 0; i < toCreate; i++)
-		{
-			GpuPool.Enqueue(contexts[i]);
+			KernelContainer kernels = GpuKernelPool.GetKernels();
+			context = new PooledContext(kernels);
+			contexts[i] = context;
+			pool.Enqueue(context);
 		}
 
 		WarmedGpuContextCount = threadCount;
@@ -85,7 +87,6 @@ public static class GpuContextPool
 			gpu.Dispose();
 		}
 
-		DivisorGeneratorGpuCache.DisposeAll();
 		PrimeTester.DisposeGpuContexts();
 		WarmedGpuContextCount = 0;
 	}
@@ -96,21 +97,30 @@ public static class GpuContextPool
 		GpuPool.Enqueue(ctx);
 	}
 
-	public readonly struct GpuContextLease
+	public struct GpuContextLease
 	{
-		private readonly PooledContext _ctx;
+		private PooledContext? _ctx;
+		private bool _disposed;
 
 		internal GpuContextLease(PooledContext ctx)
 		{
 			_ctx = ctx;
+			_disposed = false;
 		}
 
-		public readonly Context Context => _ctx.Context;
-		public readonly Accelerator Accelerator => _ctx.Accelerator;
+		public Context Context => _ctx!.Context;
+		public Accelerator Accelerator => _ctx!.Accelerator;
 
-		public readonly void Dispose()
+		public void Dispose()
 		{
-			Return(_ctx);
+			if (_disposed)
+			{
+				return;
+			}
+
+			_disposed = true;
+			Return(_ctx!);
+			_ctx = null;
 		}
 	}
 }

PerfectNumbers.Core/Gpu/GpuKernelLease.cs

@@ -2,6 +2,7 @@ using System;
 using System.Collections.Concurrent;
 using ILGPU;
 using ILGPU.Runtime;
+using PerfectNumbers.Core;
 using static PerfectNumbers.Core.Gpu.GpuContextPool;
 
 namespace PerfectNumbers.Core.Gpu;
@@ -10,25 +11,24 @@ public sealed class GpuKernelLease
 {
     private static readonly ConcurrentQueue<GpuKernelLease> Pool = new();
 
-    private Action? _releaseLimiter;
     private GpuContextLease? _gpu;
-    private KernelContainer? _kernels;
+    public KernelContainer Kernels;
     private AcceleratorStream? _stream;
 
-    private GpuKernelLease()
-    {
-    }
+    private GpuKernelLease(KernelContainer kernels)
+	{
+		Kernels = kernels;
+	}
 
-    internal static GpuKernelLease Rent(GpuPrimeWorkLimiter.Lease limiter, GpuContextLease gpu, KernelContainer kernels)
+    internal static GpuKernelLease Rent(GpuContextLease gpu, KernelContainer kernels)
     {
         if (!Pool.TryDequeue(out var lease))
         {
-            lease = new GpuKernelLease();
+            lease = new GpuKernelLease(kernels);
         }
 
-        lease._releaseLimiter = limiter.Dispose;
         lease._gpu = gpu;
-        lease._kernels = kernels;
+        lease.Kernels = kernels;
         lease._stream = null;
         return lease;
     }
@@ -56,7 +56,7 @@ public sealed class GpuKernelLease
         get
         {
             var accel = Accelerator; // avoid capturing 'this' in lambda
-            return KernelContainer.InitOnce(ref _kernels!.Order, () =>
+            return KernelContainer.InitOnce(ref Kernels!.Order, () =>
             {
                 var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>(OrderKernels.OrderKernelScan);
                 var kernel = KernelUtil.GetKernel(loaded);
@@ -70,7 +70,7 @@ public sealed class GpuKernelLease
         get
         {
             var accel = Accelerator;
-            return KernelContainer.InitOnce(ref _kernels!.Pow2Mod, () =>
+            return KernelContainer.InitOnce(ref Kernels!.Pow2Mod, () =>
             {
                 var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernels.Pow2ModKernelScan);
                 var kernel = KernelUtil.GetKernel(loaded);
@@ -84,7 +84,7 @@ public sealed class GpuKernelLease
         get
         {
             var accel = Accelerator;
-            return KernelContainer.InitOnce(ref _kernels!.Incremental, () =>
+            return KernelContainer.InitOnce(ref Kernels!.Incremental, () =>
             {
                 var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernels.IncrementalKernelScan);
                 var kernel = KernelUtil.GetKernel(loaded);
@@ -98,7 +98,7 @@ public sealed class GpuKernelLease
         get
         {
             var accel = Accelerator;
-            return KernelContainer.InitOnce(ref _kernels!.IncrementalOrder, () =>
+            return KernelContainer.InitOnce(ref Kernels!.IncrementalOrder, () =>
             {
                 var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(IncrementalKernels.IncrementalOrderKernelScan);
                 var kernel = KernelUtil.GetKernel(loaded);
@@ -112,7 +112,7 @@ public sealed class GpuKernelLease
         get
         {
             var accel = Accelerator;
-            return KernelContainer.InitOnce(ref _kernels!.Pow2ModOrder, () =>
+            return KernelContainer.InitOnce(ref Kernels!.Pow2ModOrder, () =>
             {
                 var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>(Pow2ModKernels.Pow2ModOrderKernelScan);
                 var kernel = KernelUtil.GetKernel(loaded);
@@ -121,12 +121,26 @@ public sealed class GpuKernelLease
         }
     }
 
+    public Action<AcceleratorStream, Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ushort, Stride1D.Dense>> SpecialMaxKernel
+    {
+        get
+        {
+            var accel = Accelerator;
+            return KernelContainer.InitOnce(ref Kernels!.SpecialMax, () =>
+            {
+                var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ushort, Stride1D.Dense>>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel);
+                var kernel = KernelUtil.GetKernel(loaded);
+                return kernel.CreateLauncherDelegate<Action<AcceleratorStream, Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ushort, Stride1D.Dense>>>();
+            });
+        }
+    }
+
     public Action<AcceleratorStream, Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> SmallPrimeFactorKernel
     {
         get
         {
             var accel = Accelerator;
-            return KernelContainer.InitOnce(ref _kernels!.SmallPrimeFactor, () =>
+            return KernelContainer.InitOnce(ref Kernels!.SmallPrimeFactor, () =>
             {
                 var loaded = accel.LoadAutoGroupedStreamKernel<Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan);
                 var kernel = KernelUtil.GetKernel(loaded);
@@ -143,10 +157,8 @@ public sealed class GpuKernelLease
         _gpu?.Dispose();
         _gpu = null;
 
-        _releaseLimiter?.Invoke();
-        _releaseLimiter = null;
-
-        _kernels = null;
+		Kernels?.Dispose();
+		GpuPrimeWorkLimiter.Release();
         Pool.Enqueue(this);
     }
 

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

@@ -1,3 +1,4 @@
+using PerfectNumbers.Core;
 using System.Collections.Concurrent;
 using ILGPU;
 using ILGPU.Runtime;
@@ -7,297 +8,315 @@ 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 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 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<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>? Order;
-    public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
-    ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>? Incremental;
-    public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
-        ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>,
-        ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>,
-        ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2Mod;
-    public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
-            ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? IncrementalOrder;
-    public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
-        ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2ModOrder;
-    public Action<AcceleratorStream, Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>? SmallPrimeFactor;
+	// Serializes first-time initialization of kernels/buffers per accelerator.
+	public Action<AcceleratorStream, Index1D, ulong, ulong, ArrayView<GpuUInt128>, ArrayView<ulong>>? Order;
+	public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+	ulong, ulong, ulong, ulong, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>>? Incremental;
+	public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+		ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>,
+		ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>,
+		ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2Mod;
+	public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+			ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? IncrementalOrder;
+	public Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong,
+		ResidueAutomatonArgs, ArrayView<int>, ArrayView1D<ulong, Stride1D.Dense>>? Pow2ModOrder;
+	public Action<AcceleratorStream, Index1D, ulong, uint, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, int, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<int, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>>? SmallPrimeFactor;
+	public Action<AcceleratorStream, Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ushort, Stride1D.Dense>>? SpecialMax;
 
-    // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
-    public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallCycles;
-    public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimeFactorsPrimes;
-    public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimeFactorsSquares;
-    public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimeFactorPrimeSlots;
-    public MemoryBuffer1D<int, Stride1D.Dense>? SmallPrimeFactorExponentSlots;
-    public MemoryBuffer1D<int, Stride1D.Dense>? SmallPrimeFactorCountSlot;
-    public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimeFactorRemainingSlot;
-    public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastOne;
-    public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastOne;
-    public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastSeven;
-    public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastSeven;
+	// Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallCycles;
+	public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimeFactorsPrimes;
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimeFactorsSquares;
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimeFactorPrimeSlots;
+	public MemoryBuffer1D<int, Stride1D.Dense>? SmallPrimeFactorExponentSlots;
+	public MemoryBuffer1D<int, Stride1D.Dense>? SmallPrimeFactorCountSlot;
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimeFactorRemainingSlot;
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SpecialMaxFactors;
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SpecialMaxCandidates;
+	public MemoryBuffer1D<ushort, Stride1D.Dense>? SpecialMaxResult;
+	public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastOne;
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastOne;
+	public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastSeven;
+	public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastSeven;
 
-    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;
-        }
+	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;
-    }
+		current = Volatile.Read(ref slot);
+		if (current is null)
+		{
+			current = factory();
+			Volatile.Write(ref slot, current);
+		}
+		return current;
+	}
+
+	public void Dispose() => GpuKernelPool.Return(this);
 }
 
 public readonly struct ResiduePrimeViews(
-    ArrayView1D<uint, Stride1D.Dense> lastOne,
-    ArrayView1D<uint, Stride1D.Dense> lastSeven,
-    ArrayView1D<ulong, Stride1D.Dense> lastOnePow2,
-    ArrayView1D<ulong, Stride1D.Dense> lastSevenPow2)
+	ArrayView1D<uint, Stride1D.Dense> lastOne,
+	ArrayView1D<uint, Stride1D.Dense> lastSeven,
+	ArrayView1D<ulong, Stride1D.Dense> lastOnePow2,
+	ArrayView1D<ulong, Stride1D.Dense> lastSevenPow2)
 {
-    public readonly ArrayView1D<uint, Stride1D.Dense> LastOne = lastOne;
-    public readonly ArrayView1D<uint, Stride1D.Dense> LastSeven = lastSeven;
-    public readonly ArrayView1D<ulong, Stride1D.Dense> LastOnePow2 = lastOnePow2;
-    public readonly ArrayView1D<ulong, Stride1D.Dense> LastSevenPow2 = lastSevenPow2;
+	public readonly ArrayView1D<uint, Stride1D.Dense> LastOne = lastOne;
+	public readonly ArrayView1D<uint, Stride1D.Dense> LastSeven = lastSeven;
+	public readonly ArrayView1D<ulong, Stride1D.Dense> LastOnePow2 = lastOnePow2;
+	public readonly ArrayView1D<ulong, Stride1D.Dense> LastSevenPow2 = lastSevenPow2;
 }
 
-public readonly struct SmallPrimeFactorTables
+public readonly struct SmallPrimeFactorTables(
+	MemoryBuffer1D<uint, Stride1D.Dense> primes,
+	MemoryBuffer1D<ulong, Stride1D.Dense> squares,
+	int count)
 {
-    public readonly MemoryBuffer1D<uint, Stride1D.Dense> Primes;
-    public readonly MemoryBuffer1D<ulong, Stride1D.Dense> Squares;
-    public readonly int Count;
+	public readonly MemoryBuffer1D<uint, Stride1D.Dense> Primes = primes;
+	public readonly MemoryBuffer1D<ulong, Stride1D.Dense> Squares = squares;
+	public readonly int Count = count;
 
-    public SmallPrimeFactorTables(
-        MemoryBuffer1D<uint, Stride1D.Dense> primes,
-        MemoryBuffer1D<ulong, Stride1D.Dense> squares,
-        int count)
-    {
-        Primes = primes;
-        Squares = squares;
-        Count = count;
-    }
-
-    public ArrayView1D<uint, Stride1D.Dense> PrimesView => Primes.View;
-    public ArrayView1D<ulong, Stride1D.Dense> SquaresView => Squares.View;
+	public ArrayView1D<uint, Stride1D.Dense> PrimesView => Primes.View;
+	public ArrayView1D<ulong, Stride1D.Dense> SquaresView => Squares.View;
 }
 
-public readonly struct SmallPrimeFactorScratch
+public readonly struct SmallPrimeFactorScratch(
+		MemoryBuffer1D<ulong, Stride1D.Dense> primeSlots,
+		MemoryBuffer1D<int, Stride1D.Dense> exponentSlots,
+		MemoryBuffer1D<int, Stride1D.Dense> countSlot,
+		MemoryBuffer1D<ulong, Stride1D.Dense> remainingSlot)
 {
-    public readonly MemoryBuffer1D<ulong, Stride1D.Dense> PrimeSlots;
-    public readonly MemoryBuffer1D<int, Stride1D.Dense> ExponentSlots;
-    public readonly MemoryBuffer1D<int, Stride1D.Dense> CountSlot;
-    public readonly MemoryBuffer1D<ulong, Stride1D.Dense> RemainingSlot;
+	public readonly MemoryBuffer1D<ulong, Stride1D.Dense> PrimeSlots = primeSlots;
+	public readonly MemoryBuffer1D<int, Stride1D.Dense> ExponentSlots = exponentSlots;
+	public readonly MemoryBuffer1D<int, Stride1D.Dense> CountSlot = countSlot;
+	public readonly MemoryBuffer1D<ulong, Stride1D.Dense> RemainingSlot = remainingSlot;
 
-    public SmallPrimeFactorScratch(
-        MemoryBuffer1D<ulong, Stride1D.Dense> primeSlots,
-        MemoryBuffer1D<int, Stride1D.Dense> exponentSlots,
-        MemoryBuffer1D<int, Stride1D.Dense> countSlot,
-        MemoryBuffer1D<ulong, Stride1D.Dense> remainingSlot)
-    {
-        PrimeSlots = primeSlots;
-        ExponentSlots = exponentSlots;
-        CountSlot = countSlot;
-        RemainingSlot = remainingSlot;
-    }
+	public void Clear()
+	{
+		PrimeSlots.MemSetToZero();
+		ExponentSlots.MemSetToZero();
+		CountSlot.MemSetToZero();
+		RemainingSlot.MemSetToZero();
+	}
+}
 
-    public void Clear()
-    {
-        PrimeSlots.MemSetToZero();
-        ExponentSlots.MemSetToZero();
-        CountSlot.MemSetToZero();
-        RemainingSlot.MemSetToZero();
-    }
+public readonly struct SpecialMaxScratch(
+	MemoryBuffer1D<ulong, Stride1D.Dense> factorValues,
+	MemoryBuffer1D<ulong, Stride1D.Dense> candidateValues,
+	MemoryBuffer1D<ushort, Stride1D.Dense> resultSlot)
+{
+	public readonly MemoryBuffer1D<ulong, Stride1D.Dense> FactorValues = factorValues;
+	public readonly MemoryBuffer1D<ulong, Stride1D.Dense> CandidateValues = candidateValues;
+	public readonly MemoryBuffer1D<ushort, Stride1D.Dense> ResultSlot = resultSlot;
+
+	public ArrayView1D<ulong, Stride1D.Dense> FactorsView => FactorValues.View;
+	public ArrayView1D<ulong, Stride1D.Dense> CandidatesView => CandidateValues.View;
+	public ArrayView1D<ushort, Stride1D.Dense> ResultView => ResultSlot.View;
 }
 
 public class GpuKernelPool
 {
-    private static readonly ConcurrentDictionary<Accelerator, KernelContainer> KernelCache = new(); // TODO: Replace this concurrent map with a simple accelerator-indexed lookup once kernel launchers are prewarmed during startup so we can drop the thread-safe wrapper entirely.
+	private static readonly ConcurrentQueue<KernelContainer> KernelCache = new(); // TODO: Replace this concurrent map with a simple accelerator-indexed lookup once kernel launchers are prewarmed during startup so we can drop the thread-safe wrapper entirely.
 
-    private static KernelContainer GetKernels(Accelerator accelerator)
-    {
-        return KernelCache.GetOrAdd(accelerator, _ => new KernelContainer());
-    }
+	public static KernelContainer GetKernels() => KernelCache.TryDequeue(out var container)
+		? container
+		: 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<ulong, Stride1D.Dense> EnsureSmallCyclesOnDevice(Accelerator accelerator)
-    {
-        var kernels = GetKernels(accelerator);
-        if (kernels.SmallCycles is { } buffer)
-        {
-            return buffer.View;
-        }
+	public static void Return(KernelContainer container) => KernelCache.Enqueue(container);
 
-        // Ensure single upload per accelerator even if multiple threads race here.
-        lock (kernels) // TODO: Remove this lock by pre-uploading the immutable small-cycle snapshot during initialization; once no mutation happens at runtime, the pool must expose a simple reference without synchronization.
-        {
-            if (kernels.SmallCycles is { } existing)
-            {
-                return existing.View;
-            }
+	internal static void PreloadStaticTables(KernelContainer kernels, Accelerator accelerator)
+	{
+		EnsureSmallCyclesOnDevice(kernels, accelerator);
+		EnsureSmallPrimesOnDevice(kernels, accelerator);
+		EnsureSmallPrimeFactorTables(kernels, accelerator);
+	}
 
-            var host = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); // TODO: Preload this device buffer during startup and keep it immutable so we can delete the lock above in favor of the preloaded snapshot.
-            var device = accelerator.Allocate1D<ulong>(host.Length);
-            device.View.CopyFromCPU(host);
-            kernels.SmallCycles = device;
-            return device.View;
-        }
-    }
+	// 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<ulong, Stride1D.Dense> EnsureSmallCyclesOnDevice(KernelContainer kernels, Accelerator accelerator)
+	{
+		if (kernels.SmallCycles is { } buffer)
+		{
+			return buffer.View;
+		}
 
-    public static ResiduePrimeViews EnsureSmallPrimesOnDevice(Accelerator accelerator)
-    {
-        var kernels = GetKernels(accelerator);
-        if (kernels.SmallPrimesLastOne is { } lastOne &&
-            kernels.SmallPrimesLastSeven is { } lastSeven &&
-            kernels.SmallPrimesPow2LastOne is { } lastOnePow2 &&
-            kernels.SmallPrimesPow2LastSeven is { } lastSevenPow2)
-        {
-            return new ResiduePrimeViews(lastOne.View, lastSeven.View, lastOnePow2.View, lastSevenPow2.View);
-        }
+		var host = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); // TODO: Preload this device buffer during startup and keep it immutable so we can delete the lock above in favor of the preloaded snapshot.
+		var device = accelerator.Allocate1D<ulong>(host.Length);
+		device.View.CopyFromCPU(host);
+		kernels.SmallCycles = device;
+		return device.View;
+	}
 
-        lock (kernels) // TODO: Inline these small-prime uploads into startup initialization alongside the small-cycle snapshot so we can drop runtime locking and keep the GPU pool free of synchronization.
-        {
-            if (kernels.SmallPrimesLastOne is { } existingLastOne &&
-                kernels.SmallPrimesLastSeven is { } existingLastSeven &&
-                kernels.SmallPrimesPow2LastOne is { } existingLastOnePow2 &&
-                kernels.SmallPrimesPow2LastSeven is { } existingLastSevenPow2)
-            {
-                return new ResiduePrimeViews(existingLastOne.View, existingLastSeven.View, existingLastOnePow2.View, existingLastSevenPow2.View);
-            }
+	public static ResiduePrimeViews EnsureSmallPrimesOnDevice(KernelContainer kernels, Accelerator accelerator)
+	{
+		if (kernels.SmallPrimesLastOne is { } lastOne &&
+			kernels.SmallPrimesLastSeven is { } lastSeven &&
+			kernels.SmallPrimesPow2LastOne is { } lastOnePow2 &&
+			kernels.SmallPrimesPow2LastSeven is { } lastSevenPow2)
+		{
+			return new ResiduePrimeViews(lastOne.View, lastSeven.View, lastOnePow2.View, lastSevenPow2.View);
+		}
 
-            var hostLastOne = PrimesGenerator.SmallPrimesLastOne;
-            var hostLastSeven = PrimesGenerator.SmallPrimesLastSeven;
-            var hostLastOnePow2 = PrimesGenerator.SmallPrimesPow2LastOne;
-            var hostLastSevenPow2 = PrimesGenerator.SmallPrimesPow2LastSeven;
+		if (kernels.SmallPrimesLastOne is { } existingLastOne &&
+			kernels.SmallPrimesLastSeven is { } existingLastSeven &&
+			kernels.SmallPrimesPow2LastOne is { } existingLastOnePow2 &&
+			kernels.SmallPrimesPow2LastSeven is { } existingLastSevenPow2)
+		{
+			return new ResiduePrimeViews(existingLastOne.View, existingLastSeven.View, existingLastOnePow2.View, existingLastSevenPow2.View);
+		}
 
-            var deviceLastOne = accelerator.Allocate1D<uint>(hostLastOne.Length);
-            deviceLastOne.View.CopyFromCPU(hostLastOne);
-            var deviceLastSeven = accelerator.Allocate1D<uint>(hostLastSeven.Length);
-            deviceLastSeven.View.CopyFromCPU(hostLastSeven);
-            var deviceLastOnePow2 = accelerator.Allocate1D<ulong>(hostLastOnePow2.Length);
-            deviceLastOnePow2.View.CopyFromCPU(hostLastOnePow2);
-            var deviceLastSevenPow2 = accelerator.Allocate1D<ulong>(hostLastSevenPow2.Length);
-            deviceLastSevenPow2.View.CopyFromCPU(hostLastSevenPow2);
+		var hostLastOne = PrimesGenerator.SmallPrimesLastOne;
+		var hostLastSeven = PrimesGenerator.SmallPrimesLastSeven;
+		var hostLastOnePow2 = PrimesGenerator.SmallPrimesPow2LastOne;
+		var hostLastSevenPow2 = PrimesGenerator.SmallPrimesPow2LastSeven;
 
-            kernels.SmallPrimesLastOne = deviceLastOne;
-            kernels.SmallPrimesLastSeven = deviceLastSeven;
-            kernels.SmallPrimesPow2LastOne = deviceLastOnePow2;
-            kernels.SmallPrimesPow2LastSeven = deviceLastSevenPow2;
+		var deviceLastOne = accelerator.Allocate1D<uint>(hostLastOne.Length);
+		deviceLastOne.View.CopyFromCPU(hostLastOne);
+		var deviceLastSeven = accelerator.Allocate1D<uint>(hostLastSeven.Length);
+		deviceLastSeven.View.CopyFromCPU(hostLastSeven);
+		var deviceLastOnePow2 = accelerator.Allocate1D<ulong>(hostLastOnePow2.Length);
+		deviceLastOnePow2.View.CopyFromCPU(hostLastOnePow2);
+		var deviceLastSevenPow2 = accelerator.Allocate1D<ulong>(hostLastSevenPow2.Length);
+		deviceLastSevenPow2.View.CopyFromCPU(hostLastSevenPow2);
 
-            return new ResiduePrimeViews(deviceLastOne.View, deviceLastSeven.View, deviceLastOnePow2.View, deviceLastSevenPow2.View);
-        }
-    }
+		kernels.SmallPrimesLastOne = deviceLastOne;
+		kernels.SmallPrimesLastSeven = deviceLastSeven;
+		kernels.SmallPrimesPow2LastOne = deviceLastOnePow2;
+		kernels.SmallPrimesPow2LastSeven = deviceLastSevenPow2;
 
-    public static SmallPrimeFactorTables EnsureSmallPrimeFactorTables(Accelerator accelerator)
-    {
-        var kernels = GetKernels(accelerator);
-        if (kernels.SmallPrimeFactorsPrimes is { } primeBuffer && kernels.SmallPrimeFactorsSquares is { } squareBuffer)
-        {
-            return new SmallPrimeFactorTables(primeBuffer, squareBuffer, (int)primeBuffer.Length);
-        }
+		return new ResiduePrimeViews(deviceLastOne.View, deviceLastSeven.View, deviceLastOnePow2.View, deviceLastSevenPow2.View);
+	}
 
-        lock (kernels)
-        {
-            if (kernels.SmallPrimeFactorsPrimes is { } existingPrimes && kernels.SmallPrimeFactorsSquares is { } existingSquares)
-            {
-                return new SmallPrimeFactorTables(existingPrimes, existingSquares, (int)existingPrimes.Length);
-            }
+	public static SmallPrimeFactorTables EnsureSmallPrimeFactorTables(KernelContainer kernels, Accelerator accelerator)
+	{
+		if (kernels.SmallPrimeFactorsPrimes is { } primeBuffer && kernels.SmallPrimeFactorsSquares is { } squareBuffer)
+		{
+			return new SmallPrimeFactorTables(primeBuffer, squareBuffer, (int)primeBuffer.Length);
+		}
 
-            var hostPrimes = PrimesGenerator.SmallPrimes;
-            var hostSquares = PrimesGenerator.SmallPrimesPow2;
+		if (kernels.SmallPrimeFactorsPrimes is { } existingPrimes && kernels.SmallPrimeFactorsSquares is { } existingSquares)
+		{
+			return new SmallPrimeFactorTables(existingPrimes, existingSquares, (int)existingPrimes.Length);
+		}
 
-            var devicePrimes = accelerator.Allocate1D<uint>(hostPrimes.Length);
-            devicePrimes.View.CopyFromCPU(hostPrimes);
-            var deviceSquares = accelerator.Allocate1D<ulong>(hostSquares.Length);
-            deviceSquares.View.CopyFromCPU(hostSquares);
+		var hostPrimes = PrimesGenerator.SmallPrimes;
+		var hostSquares = PrimesGenerator.SmallPrimesPow2;
 
-            kernels.SmallPrimeFactorsPrimes = devicePrimes;
-            kernels.SmallPrimeFactorsSquares = deviceSquares;
+		var devicePrimes = accelerator.Allocate1D<uint>(hostPrimes.Length);
+		devicePrimes.View.CopyFromCPU(hostPrimes);
+		var deviceSquares = accelerator.Allocate1D<ulong>(hostSquares.Length);
+		deviceSquares.View.CopyFromCPU(hostSquares);
 
-            return new SmallPrimeFactorTables(devicePrimes, deviceSquares, hostPrimes.Length);
-        }
-    }
+		kernels.SmallPrimeFactorsPrimes = devicePrimes;
+		kernels.SmallPrimeFactorsSquares = deviceSquares;
 
-    public static SmallPrimeFactorScratch EnsureSmallPrimeFactorScratch(Accelerator accelerator, int slotCount)
-    {
-        var kernels = GetKernels(accelerator);
-        lock (kernels)
-        {
-            MemoryBuffer1D<ulong, Stride1D.Dense>? primeSlots = kernels.SmallPrimeFactorPrimeSlots;
-            if (primeSlots is null || primeSlots.Length < slotCount)
-            {
-                primeSlots?.Dispose();
-                primeSlots = accelerator.Allocate1D<ulong>(slotCount);
-                kernels.SmallPrimeFactorPrimeSlots = primeSlots;
-            }
+		return new SmallPrimeFactorTables(devicePrimes, deviceSquares, hostPrimes.Length);
+	}
 
-            MemoryBuffer1D<int, Stride1D.Dense>? exponentSlots = kernels.SmallPrimeFactorExponentSlots;
-            if (exponentSlots is null || exponentSlots.Length < slotCount)
-            {
-                exponentSlots?.Dispose();
-                exponentSlots = accelerator.Allocate1D<int>(slotCount);
-                kernels.SmallPrimeFactorExponentSlots = exponentSlots;
-            }
+	public static SmallPrimeFactorScratch EnsureSmallPrimeFactorScratch(KernelContainer kernels, Accelerator accelerator, int slotCount)
+	{
+		MemoryBuffer1D<ulong, Stride1D.Dense>? primeSlots = kernels.SmallPrimeFactorPrimeSlots;
+		if (primeSlots is null || primeSlots.Length < slotCount)
+		{
+			primeSlots?.Dispose();
+			primeSlots = accelerator.Allocate1D<ulong>(slotCount);
+			kernels.SmallPrimeFactorPrimeSlots = primeSlots;
+		}
 
-            MemoryBuffer1D<int, Stride1D.Dense>? countSlot = kernels.SmallPrimeFactorCountSlot;
-            if (countSlot is null)
-            {
-                countSlot = accelerator.Allocate1D<int>(1);
-                kernels.SmallPrimeFactorCountSlot = countSlot;
-            }
+		MemoryBuffer1D<int, Stride1D.Dense>? exponentSlots = kernels.SmallPrimeFactorExponentSlots;
+		if (exponentSlots is null || exponentSlots.Length < slotCount)
+		{
+			exponentSlots?.Dispose();
+			exponentSlots = accelerator.Allocate1D<int>(slotCount);
+			kernels.SmallPrimeFactorExponentSlots = exponentSlots;
+		}
 
-            MemoryBuffer1D<ulong, Stride1D.Dense>? remainingSlot = kernels.SmallPrimeFactorRemainingSlot;
-            if (remainingSlot is null)
-            {
-                remainingSlot = accelerator.Allocate1D<ulong>(1);
-                kernels.SmallPrimeFactorRemainingSlot = remainingSlot;
-            }
+		MemoryBuffer1D<int, Stride1D.Dense>? countSlot = kernels.SmallPrimeFactorCountSlot;
+		if (countSlot is null)
+		{
+			countSlot = accelerator.Allocate1D<int>(1);
+			kernels.SmallPrimeFactorCountSlot = countSlot;
+		}
 
-            return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!);
-        }
-    }
+		MemoryBuffer1D<ulong, Stride1D.Dense>? remainingSlot = kernels.SmallPrimeFactorRemainingSlot;
+		if (remainingSlot is null)
+		{
+			remainingSlot = accelerator.Allocate1D<ulong>(1);
+			kernels.SmallPrimeFactorRemainingSlot = remainingSlot;
+		}
 
-    public static GpuKernelLease GetKernel(bool useGpuOrder)
-    {
-        var limiter = GpuPrimeWorkLimiter.Acquire();
-        var gpu = Rent();
-        var accelerator = gpu.Accelerator;
-        var kernels = GetKernels(accelerator);
-        return GpuKernelLease.Rent(limiter, gpu, kernels);
-    }
+		return new SmallPrimeFactorScratch(primeSlots!, exponentSlots!, countSlot!, remainingSlot!);
+	}
 
-    /// <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)
-    {
-        var lease = GetKernel(useGpuOrder: true);
-        var accelerator = lease.Accelerator;
-        var stream = accelerator.CreateStream();
-        action(accelerator, stream);
-        stream.Dispose();
-        lease.Dispose();
+	public static SpecialMaxScratch EnsureSpecialMaxScratch(KernelContainer kernels, Accelerator accelerator, int factorCapacity)
+	{
+		MemoryBuffer1D<ulong, Stride1D.Dense>? factorValues = kernels.SpecialMaxFactors;
+		if (factorValues is null || factorValues.Length < factorCapacity)
+		{
+			factorValues?.Dispose();
+			factorValues = accelerator.Allocate1D<ulong>(factorCapacity);
+			kernels.SpecialMaxFactors = factorValues;
+		}
 
-    }
+		MemoryBuffer1D<ulong, Stride1D.Dense>? candidateValues = kernels.SpecialMaxCandidates;
+		if (candidateValues is null || candidateValues.Length < factorCapacity)
+		{
+			candidateValues?.Dispose();
+			candidateValues = accelerator.Allocate1D<ulong>(factorCapacity);
+			kernels.SpecialMaxCandidates = candidateValues;
+		}
+
+		MemoryBuffer1D<ushort, Stride1D.Dense>? resultSlot = kernels.SpecialMaxResult;
+		if (resultSlot is null)
+		{
+			resultSlot = accelerator.Allocate1D<ushort>(1);
+			kernels.SpecialMaxResult = resultSlot;
+		}
+
+		return new SpecialMaxScratch(factorValues!, candidateValues!, resultSlot!);
+	}
+
+	public static GpuKernelLease GetKernel(bool useGpuOrder)
+	{
+		GpuPrimeWorkLimiter.Acquire();
+		var gpu = Rent();
+		var kernels = GetKernels();
+		return GpuKernelLease.Rent(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)
+	{
+		var lease = GetKernel(useGpuOrder: true);
+		var accelerator = lease.Accelerator;
+		var stream = accelerator.CreateStream();
+		action(accelerator, stream);
+		stream.Dispose();
+		lease.Dispose();
+
+	}
 }

PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs

@@ -1,3 +1,6 @@
+using System;
+using System.Threading;
+
 namespace PerfectNumbers.Core.Gpu;
 
 public static class GpuPrimeWorkLimiter
@@ -5,41 +8,26 @@ public static class GpuPrimeWorkLimiter
     private static SemaphoreSlim _semaphore = new(1, int.MaxValue);
     private static int _currentLimit = 1;
 
-    public static Lease Acquire()
+	public static void Acquire() => _semaphore.Wait();
+	public static void Release() => _semaphore.Release();
+
+	public static void SetLimit(int value)
     {
-        // TODO: Inline the pooled limiter guard from the GpuLimiterThroughputBenchmarks so prime-sieve
-        // GPU jobs stop allocating Lease instances on every acquisition.
-        var sem = _semaphore;
-        // sem.Wait();
-        return new Lease(sem);
-    }
-
-    public static void SetLimit(int value)
-    {
-		if (value == _currentLimit)
-		{
-			return;
-		}
-
-		// TODO: Switch to a shared limiter implementation with GpuWorkLimiter so we can coordinate CPU/GPU prime work using
-		// the same semaphore pool without reallocating per adjustment.
-		_semaphore = new SemaphoreSlim(value, value);
-		_currentLimit = value;
-    }
-
-    public sealed class Lease
-    {
-        private readonly SemaphoreSlim _sem;
-
-        internal Lease(SemaphoreSlim sem)
+        if (value <= 0)
         {
-            _sem = sem;
+            value = 1;
         }
 
-        public void Dispose()
+        if (value == _currentLimit)
         {
-            // _sem.Release();
+            return;
         }
+
+        // TODO: Switch to a shared limiter implementation with GpuWorkLimiter so we can coordinate CPU/GPU prime work using
+        // the same semaphore pool without reallocating per adjustment.
+        var previous = _semaphore;
+        _semaphore = new SemaphoreSlim(value, value);
+        _currentLimit = value;
+        previous.Dispose();
     }
 }
-

PerfectNumbers.Core/Gpu/GpuStaticTableInitializer.cs

@@ -0,0 +1,39 @@
+using System.Collections.Concurrent;
+using System.Collections.Generic;
+using System.Runtime.CompilerServices;
+using ILGPU;
+using ILGPU.Runtime;
+using PerfectNumbers.Core;
+
+namespace PerfectNumbers.Core.Gpu;
+
+internal static class GpuStaticTableInitializer
+{
+    private static readonly ConcurrentDictionary<Accelerator, byte> Initialized = new(AcceleratorReferenceComparer.Instance);
+
+    internal static void EnsureStaticTables(KernelContainer kernels, Accelerator accelerator)
+    {
+        Initialized.GetOrAdd(accelerator, acc =>
+        {
+            PrimeTester.PrimeTesterGpuContextPool.EnsureStaticTables(kernels, acc);
+            GpuKernelPool.PreloadStaticTables(kernels, acc);
+            PrimeOrderGpuHeuristics.PreloadStaticTables( acc);
+            return 0;
+        });
+    }
+
+    private sealed class AcceleratorReferenceComparer : IEqualityComparer<Accelerator>
+    {
+        internal static AcceleratorReferenceComparer Instance { get; } = new();
+
+        public bool Equals(Accelerator? x, Accelerator? y)
+        {
+            return ReferenceEquals(x, y);
+        }
+
+        public int GetHashCode(Accelerator obj)
+        {
+            return RuntimeHelpers.GetHashCode(obj);
+        }
+    }
+}

PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs

@@ -0,0 +1,80 @@
+using ILGPU;
+using ILGPU.Runtime;
+using PerfectNumbers.Core;
+
+namespace PerfectNumbers.Core.Gpu;
+
+internal static partial class PrimeOrderGpuHeuristics
+{
+    internal static void EvaluateSpecialMaxCandidatesKernel(
+        Index1D index,
+        ulong phi,
+        ArrayView1D<ulong, Stride1D.Dense> factors,
+        int factorCount,
+        MontgomeryDivisorData divisor,
+        ArrayView1D<ulong, Stride1D.Dense> candidates,
+        ArrayView1D<ushort, Stride1D.Dense> resultOut)
+    {
+        if (index != 0)
+        {
+            return;
+        }
+
+        int actual = 0;
+        for (int i = 0; i < factorCount; i++)
+        {
+            ulong factor = factors[i];
+            if (factor <= 1UL)
+            {
+                continue;
+            }
+
+            ulong reduced = phi / factor;
+            if (reduced == 0UL)
+            {
+                continue;
+            }
+
+            candidates[actual] = reduced;
+            actual++;
+        }
+
+        if (actual == 0)
+        {
+            resultOut[0] = 1;
+            return;
+        }
+
+        for (int i = 1; i < actual; i++)
+        {
+            ulong current = candidates[i];
+            int insert = i - 1;
+            while (insert >= 0 && candidates[insert] > current)
+            {
+                candidates[insert + 1] = candidates[insert];
+                insert--;
+            }
+
+            candidates[insert + 1] = current;
+        }
+
+        var divisorPartial = new GpuDivisorPartialData(divisor.Modulus);
+        var stepper = new ExponentRemainderStepperGpu(divisorPartial);
+        if (stepper.InitializeIsUnityGpu(candidates[0]))
+        {
+            resultOut[0] = 0;
+            return;
+        }
+
+        for (int i = 1; i < actual; i++)
+        {
+            if (stepper.ComputeNextIsUnityGpu(candidates[i]))
+            {
+                resultOut[0] = 0;
+                return;
+            }
+        }
+
+        resultOut[0] = 1;
+    }
+}

PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs

@@ -92,14 +92,20 @@ internal static class PrimeTesterKernels
 
         int divisorLength = (int)divisorSquares.Length;
         int threadIndex = index;
-        bool skipThread = threadIndex >= divisorLength;
-        ulong divisorSquare = skipThread ? 0UL : divisorSquares[threadIndex];
-        skipThread = skipThread | (divisorSquare > maxDivisorSquare);
+        if (threadIndex >= divisorLength)
+        {
+            return;
+        }
 
-        ulong divisor = skipThread ? 1UL : divisors[threadIndex];
-        ulong remainder = skipThread ? 1UL : n % divisor;
-        bool hasFactor = !skipThread & (remainder == 0UL);
-        if (hasFactor)
+        ulong divisorSquare = divisorSquares[threadIndex];
+        if (divisorSquare > maxDivisorSquare)
+        {
+            return;
+        }
+
+        ulong divisor = divisors[threadIndex];
+        divisor = n % divisor;
+        if (divisor == 0UL)
         {
             resultFlag[0] = 1;
         }

PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs

@@ -31,11 +31,11 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
         ulong step3 = ((exponent % 3UL) << 1) % 3UL;
         ulong step5 = ((exponent % 5UL) << 1) % 5UL;
         GpuUInt128 twoPGpu = (GpuUInt128)twoP;
-        var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
+        var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(gpuLease.Kernels, accelerator);
         ResiduePrimeViews primeViews = default;
         if (_kernelType == GpuKernelType.Pow2Mod)
         {
-            primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
+            primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(gpuLease.Kernels, accelerator);
         }
 
         var orderBuffer = accelerator.Allocate1D<ulong>(batchSize);
@@ -121,7 +121,7 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u
         }
 
         ArrayPool<ulong>.Shared.Return(orderArray);
-        orderBuffer.Dispose();
+		orderBuffer.Dispose();
         gpuLease.Dispose();
     }
 }

PerfectNumbers.Core/Gpu/MersenneNumberLucasLehmerGpuTester.cs

@@ -54,13 +54,13 @@ public class MersenneNumberLucasLehmerGpuTester
             return false;
         }
 
-        var limiter = GpuPrimeWorkLimiter.Acquire();
+        GpuPrimeWorkLimiter.Acquire();
         bool result;
         if (exponent >= 128UL)
         {
             if (!TryGetNttParameters(exponent, out var nttMod, out var primitiveRoot))
             {
-                limiter.Dispose();
+				GpuPrimeWorkLimiter.Release();
                 throw new NotSupportedException("NTT parameters for the given exponent are not available.");
             }
 
@@ -80,7 +80,7 @@ public class MersenneNumberLucasLehmerGpuTester
             gpu.Dispose();
         }
 
-        limiter.Dispose();
+		GpuPrimeWorkLimiter.Release();
         return result;
     }
 

PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs

@@ -4,62 +4,63 @@ namespace PerfectNumbers.Core.Gpu;
 
 public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuOrder)
 {
-        private readonly GpuKernelType _kernelType = kernelType;
-        private readonly bool _useGpuOrder = useGpuOrder;
+	private readonly GpuKernelType _kernelType = kernelType;
+	private readonly bool _useGpuOrder = useGpuOrder;
 
-        public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
-        {
-                var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
-                var accelerator = gpuLease.Accelerator;
-                var stream = gpuLease.Stream;
-                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 = lastDigit == LastDigit.Seven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
+	public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
+	{
+		var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder);
+		var accelerator = gpuLease.Accelerator;
+		var stream = gpuLease.Stream;
+		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 = lastDigit == LastDigit.Seven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters
 
-                var kernel = _kernelType switch
-                {
-                        GpuKernelType.Pow2Mod => gpuLease.Pow2ModOrderKernel,
-                        _ => gpuLease.IncrementalOrderKernel,
-                };
+		var kernel = _kernelType switch
+		{
+			GpuKernelType.Pow2Mod => gpuLease.Pow2ModOrderKernel,
+			_ => gpuLease.IncrementalOrderKernel,
+		};
 
-                var foundBuffer = accelerator.Allocate1D<int>(1); // TODO: Replace this allocation with the pooled device buffer from GpuOrderKernelBenchmarks so repeated scans reuse the pinned staging memory instead of allocating per run.
-                exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
+		var foundBuffer = accelerator.Allocate1D<int>(1); // TODO: Replace this allocation with the pooled device buffer from GpuOrderKernelBenchmarks so repeated scans reuse the pinned staging memory instead of allocating per run.
+		exponent.Mod10_8_5_3Steps(out ulong step10, out ulong step8, out ulong step5, out ulong step3);
 
-                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;
-                        // This clear is required because the result may not always be set. Remove it after the kernels are modified to always set the result.
-                        foundBuffer.MemSetToZero();
-                        // Precompute residue automaton bases for this batch
-                        UInt128 q0 = twoP * kStart + 1UL;
-                        // TODO: Replace these modulo recomputations with the divisor-cycle stepping deltas captured in
-                        // MersenneDivisorCycleLengthGpuBenchmarks so each batch advances via cached residues from the
-                        // single snapshot block instead of recalculating mod 10/8/5/3 in every iteration.
-                        q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
-                        var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
+		UInt128 remaining;
+		int currentSize;
+		int found = 0;
 
-                        // Ensure device has small cycles table for early rejections in order kernels
-                        var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
-                        kernel(stream, currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, divMul,
-                                        kernelArgs, foundBuffer.View, smallCyclesView);
+		var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(gpuLease.Kernels, accelerator);
+		while (kStart <= maxK && Volatile.Read(ref isPrime))
+		{
+			remaining = maxK - kStart + 1UL;
+			currentSize = remaining > (UInt128)batchSize ? batchSize : (int)remaining;
+			// This clear is required because the result may not always be set. Remove it after the kernels are modified to always set the result.
+			foundBuffer.MemSetToZero();
+			// Precompute residue automaton bases for this batch
+			UInt128 q0 = twoP * kStart + 1UL;
+			// TODO: Replace these modulo recomputations with the divisor-cycle stepping deltas captured in
+			// MersenneDivisorCycleLengthGpuBenchmarks so each batch advances via cached residues from the
+			// single snapshot block instead of recalculating mod 10/8/5/3 in every iteration.
+			q0.Mod10_8_5_3(out ulong q0m10, out ulong q0m8, out ulong q0m5, out ulong q0m3);
+			var kernelArgs = new ResidueAutomatonArgs(q0m10, step10, q0m8, step8, q0m3, step3, q0m5, step5);
 
-                        stream.Synchronize();
-                        foundBuffer.View.CopyToCPU(ref found, 1);
-                        if (found != 0)
-                        {
-                                Volatile.Write(ref isPrime, false);
-                                break;
-                        }
+			// Ensure device has small cycles table for early rejections in order kernels
+			kernel(stream, currentSize, exponent, (GpuUInt128)twoP, (GpuUInt128)kStart, last, divMul,
+							kernelArgs, foundBuffer.View, smallCyclesView);
 
-                        kStart += (UInt128)currentSize;
-                }
+			stream.Synchronize();
+			foundBuffer.View.CopyToCPU(ref found, 1);
+			if (found != 0)
+			{
+				Volatile.Write(ref isPrime, false);
+				break;
+			}
 
-                foundBuffer.Dispose();
-                gpuLease.Dispose();
-        }
+			kStart += (UInt128)currentSize;
+		}
+
+		foundBuffer.Dispose();
+		gpuLease.Dispose();
+	}
 }

PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs

@@ -26,7 +26,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
     // GPU residue variant: check 2^p % q == 1 for q = 2*p*k + 1.
     public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime)
     {
-        var limiter = GpuPrimeWorkLimiter.Acquire();
+        GpuPrimeWorkLimiter.Acquire();
         var gpuLease = RentAccelerator();
         var accelerator = gpuLease.Accelerator;
 
@@ -34,15 +34,16 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
 
         var stream = accelerator.CreateStream();
         var kernel = GetPow2ModKernel(accelerator);
-        var resources = RentResources(accelerator, GpuConstants.ScanBatchSize);
+		var resources = RentResources(accelerator, GpuConstants.ScanBatchSize);
+		var gpuKernels = GpuKernelPool.GetKernels();
 
         var orderBuffer = resources.OrderBuffer;
         ulong[] orderArray = resources.OrderArray;
         int batchSize = resources.Capacity;
 
         // Ensure device has small cycles and primes tables for in-kernel lookup
-        var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(accelerator);
-        ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(accelerator);
+        var smallCyclesView = GpuKernelPool.EnsureSmallCyclesOnDevice(gpuKernels, accelerator);
+        ResiduePrimeViews primeViews = GpuKernelPool.EnsureSmallPrimesOnDevice(gpuKernels, accelerator);
 
         UInt128 kStart = 1UL;
         UInt128 limit = maxK + UInt128.One;
@@ -108,9 +109,10 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder)
 
         ReturnResources(accelerator, resources);
         stream.Dispose();
-        // Monitor.Exit(gpuLease.ExecutionLock);
-        ReturnAccelerator(gpuLease);
-        limiter.Dispose();
+		// Monitor.Exit(gpuLease.ExecutionLock);
+		gpuKernels.Dispose();
+		ReturnAccelerator(gpuLease);
+		GpuPrimeWorkLimiter.Release();
     }
 
     private Action<AcceleratorStream, Index1D, ulong, GpuUInt128, GpuUInt128, byte, ulong, ResidueAutomatonArgs, ArrayView<ulong>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<uint, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<ulong, Stride1D.Dense>> GetPow2ModKernel(Accelerator accelerator) =>

PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs

@@ -18,6 +18,11 @@ internal enum GpuPow2ModStatus
 
 internal static partial class PrimeOrderGpuHeuristics
 {
+    internal static void PreloadStaticTables(Accelerator accelerator)
+    {
+        _ = GetSmallPrimeDeviceCache(accelerator);
+    }
+
     private static readonly ConcurrentDictionary<ulong, byte> OverflowedPrimes = new();
     private static readonly ConcurrentDictionary<UInt128, byte> OverflowedPrimesWide = new();
     private static readonly ConcurrentDictionary<Accelerator, Action<AcceleratorStream, Index1D, ArrayView1D<ulong, Stride1D.Dense>, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>>> Pow2ModKernelCache = new();

PerfectNumbers.Core/Gpu/SharedGpuContext.cs

@@ -0,0 +1,19 @@
+using System;
+using System.Threading;
+using ILGPU;
+using ILGPU.Runtime;
+
+namespace PerfectNumbers.Core.Gpu;
+
+internal static class SharedGpuContext
+{
+    private static readonly Lazy<Context> s_context = new(() => ILGPU.Context.CreateDefault(), LazyThreadSafetyMode.ExecutionAndPublication);
+    private static readonly Lazy<Device> s_device = new(() => s_context.Value.GetPreferredDevice(false), LazyThreadSafetyMode.ExecutionAndPublication);
+    private static readonly Lazy<Accelerator> s_accelerator = new(() => s_device.Value.CreateAccelerator(s_context.Value), LazyThreadSafetyMode.ExecutionAndPublication);
+
+    internal static Context Context => s_context.Value;
+
+    internal static Device Device => s_device.Value;
+
+    internal static Accelerator Accelerator => s_accelerator.Value;
+}

PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs

@@ -1,5 +1,4 @@
 using Open.Numeric.Primes;
-using System;
 using ILGPU.Runtime;
 using System.Runtime.CompilerServices;
 using PerfectNumbers.Core.Gpu;
@@ -42,20 +41,20 @@ public sealed class HeuristicCombinedPrimeTester
 	}
 
 
-        public static void EnsureInitialized()
-        {
-        // Intentionally left blank. Accessing this method forces the static constructor to run.
-        }
-        
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding1Span => CombinedDivisorsEnding1;
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding3Span => CombinedDivisorsEnding3;
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding7Span => CombinedDivisorsEnding7;
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding9Span => CombinedDivisorsEnding9;
-        
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding1SquaresSpan => CombinedDivisorsEnding1Squares;
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding3SquaresSpan => CombinedDivisorsEnding3Squares;
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding7SquaresSpan => CombinedDivisorsEnding7Squares;
-        internal static ReadOnlySpan<ulong> CombinedDivisorsEnding9SquaresSpan => CombinedDivisorsEnding9Squares;
+	public static void EnsureInitialized()
+	{
+		// Intentionally left blank. Accessing this method forces the static constructor to run.
+	}
+
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding1Span => CombinedDivisorsEnding1;
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding3Span => CombinedDivisorsEnding3;
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding7Span => CombinedDivisorsEnding7;
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding9Span => CombinedDivisorsEnding9;
+
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding1SquaresSpan => CombinedDivisorsEnding1Squares;
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding3SquaresSpan => CombinedDivisorsEnding3Squares;
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding7SquaresSpan => CombinedDivisorsEnding7Squares;
+	internal static ReadOnlySpan<ulong> CombinedDivisorsEnding9SquaresSpan => CombinedDivisorsEnding9Squares;
 
 	private const ulong Wheel210 = 210UL;
 
@@ -223,39 +222,39 @@ public sealed class HeuristicCombinedPrimeTester
 	};
 
 
-        private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
-        {
-                var limiter = GpuPrimeWorkLimiter.Acquire();
-                var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(1);
-                var accelerator = gpu.Accelerator;
-                var kernel = gpu.HeuristicTrialDivisionKernel;
-                var flagView1D = gpu.HeuristicFlag.View;
-                var flagView = flagView1D.AsContiguous();
+	private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
+	{
+		GpuPrimeWorkLimiter.Acquire();
+		var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(1);
+		var accelerator = gpu.Accelerator;
+		var kernel = PrimeTester.PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.HeuristicTrialDivisionKernel;
+		var flagView1D = gpu.HeuristicFlag.View;
+		var flagView = flagView1D.AsContiguous();
 
-                bool compositeDetected = false;
-                int divisorLength = GetCombinedDivisors(nMod10).Length;
+		bool compositeDetected = false;
+		int divisorLength = GetCombinedDivisors(nMod10).Length;
 
-                int compositeFlag = 0;
-                flagView1D.CopyFromCPU(ref compositeFlag, 1);
-                kernel(
-                        divisorLength,
-                        flagView,
-                        n,
-                        maxDivisorSquare,
-                        HeuristicGpuDivisorTableKind.Combined,
-                        gpu.HeuristicGpuTables);
-                accelerator.Synchronize();
-                flagView1D.CopyToCPU(ref compositeFlag, 1);
-                compositeDetected = compositeFlag != 0;
+		int compositeFlag = 0;
+		flagView1D.CopyFromCPU(ref compositeFlag, 1);
+		kernel(
+				divisorLength,
+				flagView,
+				n,
+				maxDivisorSquare,
+				HeuristicGpuDivisorTableKind.Combined,
+				gpu.HeuristicGpuTables);
+		accelerator.Synchronize();
+		flagView1D.CopyToCPU(ref compositeFlag, 1);
+		compositeDetected = compositeFlag != 0;
 
-                gpu.Dispose();
-                limiter.Dispose();
-                return compositeDetected;
-        }
+		gpu.Dispose();
+		GpuPrimeWorkLimiter.Release();
+		return compositeDetected;
+	}
 
 
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        private static bool EvaluateWithOpenNumericFallback(ulong n)
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static bool EvaluateWithOpenNumericFallback(ulong n)
 	{
 		return Prime.Numbers.IsPrime(n);
 	}

PerfectNumbers.Core/HeuristicPrimeTester.cs

@@ -213,7 +213,7 @@ public sealed class HeuristicPrimeTester
 				int processed = 0;
 				while (processed < interleaveBatchSize)
 				{
-                    if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor))
+					if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor))
 					{
 						groupBHasMore = false;
 						break;
@@ -238,77 +238,77 @@ public sealed class HeuristicPrimeTester
 	}
 
 
-        private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
-        {
-                var limiter = GpuPrimeWorkLimiter.Acquire();
-                var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(1);
-                var accelerator = gpu.Accelerator;
-                var kernel = gpu.HeuristicTrialDivisionKernel;
-                var flagView1D = gpu.HeuristicFlag.View;
-                var flagView = flagView1D.AsContiguous();
+	private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10)
+	{
+		GpuPrimeWorkLimiter.Acquire();
+		var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(1);
+		var accelerator = gpu.Accelerator;
+		var kernel = PrimeTester.PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.HeuristicTrialDivisionKernel;
+		var flagView1D = gpu.HeuristicFlag.View;
+		var flagView = flagView1D.AsContiguous();
 
-                bool compositeDetected = false;
-                int compositeFlag = 0;
+		bool compositeDetected = false;
+		int compositeFlag = 0;
 
-                int groupALength = HeuristicPrimeSieves.GroupADivisors.Length;
+		int groupALength = HeuristicPrimeSieves.GroupADivisors.Length;
 
-                flagView1D.CopyFromCPU(ref compositeFlag, 1);
-                kernel(
-                        groupALength,
-                        flagView,
-                        n,
-                        maxDivisorSquare,
-                        HeuristicGpuDivisorTableKind.GroupA,
-                        gpu.HeuristicGpuTables);
-                accelerator.Synchronize();
-                flagView1D.CopyToCPU(ref compositeFlag, 1);
-                compositeDetected = compositeFlag != 0;
+		flagView1D.CopyFromCPU(ref compositeFlag, 1);
+		kernel(
+				groupALength,
+				flagView,
+				n,
+				maxDivisorSquare,
+				HeuristicGpuDivisorTableKind.GroupA,
+				gpu.HeuristicGpuTables);
+		accelerator.Synchronize();
+		flagView1D.CopyToCPU(ref compositeFlag, 1);
+		compositeDetected = compositeFlag != 0;
 
-                if (!compositeDetected)
-                {
-                        ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
-                        for (int i = 0; i < endingOrder.Length && !compositeDetected; i++)
-                        {
-                                byte ending = endingOrder[i];
-                                HeuristicGpuDivisorTableKind tableKind = ending switch
-                                {
-                                        1 => HeuristicGpuDivisorTableKind.GroupBEnding1,
-                                        7 => HeuristicGpuDivisorTableKind.GroupBEnding7,
-                                        9 => HeuristicGpuDivisorTableKind.GroupBEnding9,
-                                        _ => HeuristicGpuDivisorTableKind.GroupA,
-                                };
+		if (!compositeDetected)
+		{
+			ReadOnlySpan<byte> endingOrder = GetGroupBEndingOrder(nMod10);
+			for (int i = 0; i < endingOrder.Length && !compositeDetected; i++)
+			{
+				byte ending = endingOrder[i];
+				HeuristicGpuDivisorTableKind tableKind = ending switch
+				{
+					1 => HeuristicGpuDivisorTableKind.GroupBEnding1,
+					7 => HeuristicGpuDivisorTableKind.GroupBEnding7,
+					9 => HeuristicGpuDivisorTableKind.GroupBEnding9,
+					_ => HeuristicGpuDivisorTableKind.GroupA,
+				};
 
-                                if (tableKind == HeuristicGpuDivisorTableKind.GroupA)
-                                {
-                                        continue;
-                                }
+				if (tableKind == HeuristicGpuDivisorTableKind.GroupA)
+				{
+					continue;
+				}
 
-                                int divisorLength = GetGroupBDivisors(ending).Length;
+				int divisorLength = GetGroupBDivisors(ending).Length;
 
-                                compositeFlag = 0;
-                                flagView1D.CopyFromCPU(ref compositeFlag, 1);
-                                kernel(
-                                        divisorLength,
-                                        flagView,
-                                        n,
-                                        maxDivisorSquare,
-                                        tableKind,
-                                        gpu.HeuristicGpuTables);
-                                accelerator.Synchronize();
-                                flagView1D.CopyToCPU(ref compositeFlag, 1);
-                                compositeDetected = compositeFlag != 0;
-                        }
-                }
+				compositeFlag = 0;
+				flagView1D.CopyFromCPU(ref compositeFlag, 1);
+				kernel(
+						divisorLength,
+						flagView,
+						n,
+						maxDivisorSquare,
+						tableKind,
+						gpu.HeuristicGpuTables);
+				accelerator.Synchronize();
+				flagView1D.CopyToCPU(ref compositeFlag, 1);
+				compositeDetected = compositeFlag != 0;
+			}
+		}
 
-                gpu.Dispose();
-                limiter.Dispose();
-                return compositeDetected;
-        }
+		gpu.Dispose();
+		GpuPrimeWorkLimiter.Release();
+		return compositeDetected;
+	}
 
 
 
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        private static bool EvaluateWithOpenNumericFallback(ulong n)
+	[MethodImpl(MethodImplOptions.AggressiveInlining)]
+	private static bool EvaluateWithOpenNumericFallback(ulong n)
 	{
 		return Prime.Numbers.IsPrime(n);
 	}
@@ -327,7 +327,7 @@ public sealed class HeuristicPrimeTester
 	};
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-        private static bool TrySelectNextGroupBDivisor(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong maxDivisorSquare, out ulong divisor)
+	private static bool TrySelectNextGroupBDivisor(ReadOnlySpan<byte> endingOrder, Span<int> indices, ulong maxDivisorSquare, out ulong divisor)
 	{
 		ulong bestCandidate = ulong.MaxValue;
 		int bestEndingIndex = -1;

PerfectNumbers.Core/MersenneDivisorCycles.cs

@@ -412,7 +412,7 @@ public class MersenneDivisorCycles
 		// Console.WriteLine($"MersenneDivisorCycles.TryFactorIntoCountsInternal hits {Volatile.Read(ref _tryFactorIntoCountsInternalHits)}");
 
 		// HeuristicPrimeTester primeTester = _primeTester ??= new();
-		if (PrimeTester.IsPrime(remaining))
+		if (PrimeTester.IsPrimeGpu(remaining))
 		// if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None))
 		{
 			AddFactor(counts, remaining);
@@ -943,7 +943,7 @@ public class MersenneDivisorCycles
 		// HeuristicPrimeTester primeTester = _primeTester ??= new();
 		// if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None))
 		// if (!skipPrimeOrderHeuristic && primeTester.IsPrimeGpu(divisor))
-		if (!skipPrimeOrderHeuristic && PrimeTester.IsPrime(divisor))
+		if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeGpu(divisor))
 		{
 			ulong computedOrder = PrimeOrderCalculator.Calculate(
 					divisor,

PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs

@@ -1,4 +1,3 @@
-using System;
 using System.Buffers;
 using System.Diagnostics;
 using System.Numerics;
@@ -108,7 +107,7 @@ internal static partial class PrimeOrderCalculator
 		return candidateOrder;
 	}
 
-
+	// private static ulong specialMaxHits;
 	private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
 	{
 		int length = factors.Count;
@@ -121,16 +120,16 @@ internal static partial class PrimeOrderCalculator
 
 		ReadOnlySpan<FactorEntry> factorSpan = new(factors.Factors, 0, length);
 
-		if (length <= 32)
+		Span<ulong> stackBuffer = stackalloc ulong[length];
+
+		if (length <= 7)
 		{
-			Span<ulong> stackBuffer = stackalloc ulong[length];
 			return EvaluateSpecialMaxCandidates(stackBuffer, factorSpan, phi, prime, divisorData);
 		}
 
-		ulong[] rented = ThreadStaticPools.UlongPool.Rent(length);
-		bool result = EvaluateSpecialMaxCandidates(rented.AsSpan(0, length), factorSpan, phi, prime, divisorData);
-		ThreadStaticPools.UlongPool.Return(rented, clearArray: false);
-		return result;
+		// Atomic.Add(ref specialMaxHits, 1);
+		// Console.WriteLine($"Special max GPU hit {specialMaxHits} ({length})");
+		return EvaluateSpecialMaxCandidatesGpu(stackBuffer, factorSpan, phi, prime, divisorData);
 	}
 
 	private static bool EvaluateSpecialMaxCandidates(Span<ulong> buffer, ReadOnlySpan<FactorEntry> factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData)
@@ -1194,7 +1193,7 @@ internal static partial class PrimeOrderCalculator
 				// Console.WriteLine($"Partial factor pending hits {Volatile.Read(ref _partialFactorPendingHits)}");
 
 				// bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None);
-				bool isPrime = PrimeTester.IsPrime(composite);
+				bool isPrime = PrimeTester.IsPrimeGpu(composite);
 				// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite);
 
 				entry = entry.WithPrimality(isPrime);
@@ -1231,7 +1230,7 @@ internal static partial class PrimeOrderCalculator
 			// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
 			// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
 
-			cofactorIsPrime = PrimeTester.IsPrime(cofactor);
+			cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor);
 			// cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor);
 			// cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor);
 		}
@@ -1701,7 +1700,7 @@ internal static partial class PrimeOrderCalculator
 		// HeuristicPrimeTester tester = _tester ??= new();
 
 		// bool isPrime = primeTester.HeuristicIsPrimeGpu(value);
-		bool isPrime = PrimeTester.IsPrime(value);
+		bool isPrime = PrimeTester.IsPrimeGpu(value);
 		// bool isPrime = PrimeTester.IsPrimeGpu(value);
 		// bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value);
 		if (!knownComposite && isPrime)
@@ -1717,7 +1716,7 @@ internal static partial class PrimeOrderCalculator
 		// Console.WriteLine($"Factor completely after PollardRho hits {Volatile.Read(ref _factorCompletelyPollardRhoHits)}");
 
 		// isPrime = primeTester.HeuristicIsPrimeGpu(factor);
-		isPrime = PrimeTester.IsPrime(factor);
+		isPrime = PrimeTester.IsPrimeGpu(factor);
 		// isPrime = PrimeTester.IsPrimeGpu(factor);
 		// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor);
 		if (isPrime)
@@ -1754,7 +1753,7 @@ internal static partial class PrimeOrderCalculator
 			// Atomic.Add(ref _partialFactorCofactorHits, 1UL);
 			// Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}");
 
-			isPrime = PrimeTester.IsPrime(quotient);
+			isPrime = PrimeTester.IsPrimeGpu(quotient);
 			// isPrime = PrimeTester.IsPrimeGpu(quotient);
 			// isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient);
 		}

PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs

@@ -9,63 +9,111 @@ namespace PerfectNumbers.Core;
 
 internal static partial class PrimeOrderCalculator
 {
-        private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
+	private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2;
 
-        private const int GpuSmallPrimeFactorSlots = 64;
+	private const int GpuSmallPrimeFactorSlots = 64;
 
-        private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary<ulong, int> counts, out int factorCount, out ulong remaining)
-        {
-                var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots);
-                var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots);
-                Span<ulong> primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
-                Span<int> exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
-                remaining = value;
+	private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary<ulong, int> counts, out int factorCount, out ulong remaining)
+	{
+		var primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots);
+		var exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots);
+		Span<ulong> primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
+		Span<int> exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots);
+		remaining = value;
 
-                var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
-                Accelerator accelerator = lease.Accelerator;
-                AcceleratorStream stream = lease.Stream;
+		var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
+		Accelerator accelerator = lease.Accelerator;
+		AcceleratorStream stream = lease.Stream;
 
-                SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(accelerator);
-                SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(accelerator, GpuSmallPrimeFactorSlots);
+		SmallPrimeFactorTables tables = GpuKernelPool.EnsureSmallPrimeFactorTables(lease.Kernels, accelerator);
+		SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(lease.Kernels, accelerator, GpuSmallPrimeFactorSlots);
 
-                var kernel = lease.SmallPrimeFactorKernel;
-                kernel(
-                        stream,
-                        1,
-                        value,
-                        limit,
-                        tables.PrimesView,
-                        tables.SquaresView,
-                        tables.Count,
-                        scratch.PrimeSlots.View,
-                        scratch.ExponentSlots.View,
-                        scratch.CountSlot.View,
-                        scratch.RemainingSlot.View);
+		var kernel = lease.SmallPrimeFactorKernel;
+		kernel(
+				stream,
+				1,
+				value,
+				limit,
+				tables.PrimesView,
+				tables.SquaresView,
+				tables.Count,
+				scratch.PrimeSlots.View,
+				scratch.ExponentSlots.View,
+				scratch.CountSlot.View,
+				scratch.RemainingSlot.View);
 
-                stream.Synchronize();
+		stream.Synchronize();
 
-                factorCount = 0;
-                scratch.CountSlot.View.CopyToCPU(ref factorCount, 1);
-                factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots);
+		factorCount = 0;
+		scratch.CountSlot.View.CopyToCPU(ref factorCount, 1);
+		factorCount = Math.Min(factorCount, GpuSmallPrimeFactorSlots);
 
-                if (factorCount > 0)
-                {
-                        scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount);
-                        scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount);
-                }
+		if (factorCount > 0)
+		{
+			scratch.PrimeSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(primeBuffer), factorCount);
+			scratch.ExponentSlots.View.CopyToCPU(ref MemoryMarshal.GetReference(exponentBuffer), factorCount);
+		}
 
-                scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1);
+		scratch.RemainingSlot.View.CopyToCPU(ref remaining, 1);
 
-                for (int i = 0; i < factorCount; i++)
-                {
-                        ulong primeValue = primeBuffer[i];
-                        int exponent = exponentBuffer[i];
-                        counts.Add(primeValue, exponent);
-                }
+		for (int i = 0; i < factorCount; i++)
+		{
+			ulong primeValue = primeBuffer[i];
+			int exponent = exponentBuffer[i];
+			counts.Add(primeValue, exponent);
+		}
 
-                lease.Dispose();
-                ThreadStaticPools.UlongPool.Return(primeBufferArray);
-                ThreadStaticPools.IntPool.Return(exponentBufferArray);
-                return true;
-        }
+		lease.Dispose();
+		ThreadStaticPools.UlongPool.Return(primeBufferArray);
+		ThreadStaticPools.IntPool.Return(exponentBufferArray);
+		return true;
+	}
+
+	private static bool EvaluateSpecialMaxCandidatesGpu(
+			Span<ulong> buffer,
+			ReadOnlySpan<FactorEntry> factors,
+			ulong phi,
+			ulong prime,
+			in MontgomeryDivisorData divisorData)
+	{
+		if (factors.Length == 0)
+		{
+			return true;
+		}
+
+		int factorCount = factors.Length;
+		Span<ulong> factorSpan = buffer[..factorCount];
+		for (int i = 0; i < factorCount; i++)
+		{
+			factorSpan[i] = factors[i].Value;
+		}
+
+		var lease = GpuKernelPool.GetKernel(useGpuOrder: true);
+		Accelerator accelerator = lease.Accelerator;
+		AcceleratorStream stream = lease.Stream;
+
+		SpecialMaxScratch scratch = GpuKernelPool.EnsureSpecialMaxScratch(lease.Kernels, accelerator, factorCapacity: factorCount);
+
+		scratch.FactorsView.SubView(0, factorCount).CopyFromCPU(ref MemoryMarshal.GetReference(factorSpan), factorCount);
+
+		var kernel = lease.SpecialMaxKernel;
+		kernel(
+				stream,
+				1,
+				phi,
+				scratch.FactorsView,
+				factorCount,
+				divisorData,
+				scratch.CandidatesView,
+				scratch.ResultView);
+
+		stream.Synchronize();
+
+		ushort result = 0;
+		scratch.ResultView.CopyToCPU(ref result, 1);
+		lease.Dispose();
+
+		return result != 0;
+	}
 }
+

PerfectNumbers.Core/PrimeTester.cs

@@ -93,50 +93,43 @@ public sealed class PrimeTester
 	}
 
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
-	public bool IsPrimeGpu(ulong n)
+	public static bool IsPrimeGpu(ulong n)
 	{
-		// bool belowGpuRange = n < 31UL;
+		GpuPrimeWorkLimiter.Acquire();
+		var gpu = PrimeTesterGpuContextPool.Rent(1);
+		var accelerator = gpu.Accelerator;
+		var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.Kernel;
 
-		// if (belowGpuRange)
-		// {
-			return IsPrime(n);
-		// }
+		ulong value = n;
+		byte flag = 0;
 
-		// var limiter = GpuPrimeWorkLimiter.Acquire();
-		// var gpu = PrimeTesterGpuContextPool.Rent(1);
-		// var accelerator = gpu.Accelerator;
-		// var kernel = gpu.Kernel;
+		var input = gpu.Input;
+		var output = gpu.Output;
+		var inputView = input.View;
+		var outputView = output.View;
 
-		// ulong value = n;
-		// byte flag = 0;
+		inputView.CopyFromCPU(ref value, 1);
+		kernel(
+						1,
+						inputView,
+						gpu.DevicePrimesDefault.View,
+						gpu.DevicePrimesLastOne.View,
+						gpu.DevicePrimesLastSeven.View,
+						gpu.DevicePrimesLastThree.View,
+						gpu.DevicePrimesLastNine.View,
+						gpu.DevicePrimesPow2Default.View,
+						gpu.DevicePrimesPow2LastOne.View,
+						gpu.DevicePrimesPow2LastSeven.View,
+						gpu.DevicePrimesPow2LastThree.View,
+						gpu.DevicePrimesPow2LastNine.View,
+						outputView);
+		accelerator.Synchronize();
+		outputView.CopyToCPU(ref flag, 1);
 
-		// var input = gpu.Input;
-		// var output = gpu.Output;
-		// var inputView = input.View;
-		// var outputView = output.View;
+		gpu.Dispose();
+		GpuPrimeWorkLimiter.Release();
 
-		// inputView.CopyFromCPU(ref value, 1);
-		// kernel(
-		// 				1,
-		// 				inputView,
-		// 				gpu.DevicePrimesDefault.View,
-		// 				gpu.DevicePrimesLastOne.View,
-		// 				gpu.DevicePrimesLastSeven.View,
-		// 				gpu.DevicePrimesLastThree.View,
-		// 				gpu.DevicePrimesLastNine.View,
-		// 				gpu.DevicePrimesPow2Default.View,
-		// 				gpu.DevicePrimesPow2LastOne.View,
-		// 				gpu.DevicePrimesPow2LastSeven.View,
-		// 				gpu.DevicePrimesPow2LastThree.View,
-		// 				gpu.DevicePrimesPow2LastNine.View,
-		// 				outputView);
-		// accelerator.Synchronize();
-		// outputView.CopyToCPU(ref flag, 1);
-
-		// gpu.Dispose();
-		// limiter.Dispose();
-
-		// return flag != 0;
+		return flag != 0;
 	}
 
 	public static int GpuBatchSize { get; set; } = 262_144;
@@ -159,18 +152,7 @@ public sealed class PrimeTester
 				return;
 			}
 
-			int toWarm = target - WarmedGpuLeaseCount;
-			var leases = new PrimeTesterGpuContextPool.PrimeTesterGpuContextLease[toWarm];
-			for (int i = 0; i < toWarm; i++)
-			{
-				leases[i] = PrimeTesterGpuContextPool.Rent(GpuBatchSize);
-			}
-
-			for (int i = 0; i < toWarm; i++)
-			{
-				leases[i].Dispose();
-			}
-
+			PrimeTesterGpuContextPool.WarmUp(target, GpuBatchSize);
 			WarmedGpuLeaseCount = target;
 		}
 	}
@@ -178,10 +160,10 @@ public sealed class PrimeTester
 	[MethodImpl(MethodImplOptions.AggressiveInlining)]
 	public static void IsPrimeBatchGpu(ReadOnlySpan<ulong> values, Span<byte> results)
 	{
-		var limiter = GpuPrimeWorkLimiter.Acquire();
+		GpuPrimeWorkLimiter.Acquire();
 		var gpu = PrimeTesterGpuContextPool.Rent(GpuBatchSize);
 		var accelerator = gpu.Accelerator;
-		var kernel = gpu.Kernel;
+		var kernel = PrimeTesterGpuContextPool.PrimeTesterGpuContextLease.Kernel;
 		int totalLength = values.Length;
 		int batchSize = GpuBatchSize;
 
@@ -223,30 +205,25 @@ public sealed class PrimeTester
 		}
 
 		gpu.Dispose();
-		limiter.Dispose();
+		GpuPrimeWorkLimiter.Release();
 	}
 
 	internal static class PrimeTesterGpuContextPool
 	{
 		private static readonly ConcurrentQueue<PrimeTesterGpuContextLease> Pool = new();
+		private static readonly object WarmUpLock = new();
+		private static int WarmedLeaseCount;
 		private static readonly object SharedTablesLock = new();
-		private static readonly Dictionary<string, SharedHeuristicGpuTables> SharedTables = new();
-
-		private static string CreateAcceleratorKey(Accelerator accelerator)
-		{
-			return accelerator.AcceleratorType + ":" + accelerator.Name;
-		}
+		private static readonly Dictionary<Accelerator, SharedHeuristicGpuTables> SharedTables = new(AcceleratorReferenceComparer.Instance);
 
 		private static SharedHeuristicGpuTables RentSharedTables(Accelerator accelerator)
 		{
-			string key = CreateAcceleratorKey(accelerator);
-
 			lock (SharedTablesLock)
 			{
-				if (!SharedTables.TryGetValue(key, out var tables))
+				if (!SharedTables.TryGetValue(accelerator, out var tables))
 				{
-					tables = new SharedHeuristicGpuTables(accelerator, key);
-					SharedTables.Add(key, tables);
+					tables = new SharedHeuristicGpuTables(accelerator);
+					SharedTables.Add(accelerator, tables);
 				}
 
 				tables.AddReference();
@@ -267,21 +244,58 @@ public sealed class PrimeTester
 			}
 		}
 
-		internal static PrimeTesterGpuContextLease Rent(int minBufferCapacity = 1)
+		internal static void EnsureStaticTables(KernelContainer kernels, Accelerator accelerator)
 		{
-			if (!Pool.TryDequeue(out var lease))
+			var tables = RentSharedTables(accelerator);
+			var heuristicTables = tables.CreateHeuristicDivisorTables();
+			HeuristicGpuDivisorTables.InitializeShared(in heuristicTables);
+			ReleaseSharedTables(tables);
+		}
+
+		internal static void WarmUp(int target, int minBufferCapacity)
+		{
+			if (target <= 0)
 			{
-				lease = new PrimeTesterGpuContextLease(minBufferCapacity);
-			}
-			else
-			{
-				lease.EnsureCapacity(minBufferCapacity);
+				return;
 			}
 
-			lease.ResetReturnFlag();
-			return lease;
+			lock (WarmUpLock)
+			{
+				if (target <= WarmedLeaseCount)
+				{
+					return;
+				}
+
+				KernelContainer kernels = GpuKernelPool.GetKernels();
+				GpuStaticTableInitializer.EnsureStaticTables(kernels, SharedGpuContext.Accelerator);
+
+				int toCreate = target - WarmedLeaseCount;
+				for (int i = 0; i < toCreate; i++)
+				{
+					Pool.Enqueue(new PrimeTesterGpuContextLease(minBufferCapacity));
+				}
+
+				kernels.Dispose();
+				WarmedLeaseCount = target;
+			}
 		}
 
+                internal static PrimeTesterGpuContextLease Rent(int minBufferCapacity = 1)
+                {
+                        if (!Pool.TryDequeue(out var lease))
+                        {
+                                lease = new PrimeTesterGpuContextLease(minBufferCapacity);
+                        }
+                        else
+                        {
+                                lease.EnsureCapacity(minBufferCapacity);
+                        }
+
+                        lease.ResetReturnFlag();
+                        return lease;
+                }
+
+
 		internal static void Return(PrimeTesterGpuContextLease lease)
 		{
 			lease.Accelerator.Synchronize();
@@ -303,7 +317,11 @@ public sealed class PrimeTester
 				}
 
 				SharedTables.Clear();
+				var emptyTables = default(HeuristicGpuDivisorTables);
+				HeuristicGpuDivisorTables.InitializeShared(in emptyTables);
 			}
+
+			WarmedLeaseCount = 0;
 		}
 
 		internal static void Clear(Accelerator accelerator)
@@ -313,23 +331,42 @@ public sealed class PrimeTester
 				return;
 			}
 
-			var retained = new List<PrimeTesterGpuContextLease>();
-			while (Pool.TryDequeue(out var lease))
+			lock (WarmUpLock)
 			{
-				if (lease.Accelerator == accelerator)
+				var retained = new List<PrimeTesterGpuContextLease>();
+				while (Pool.TryDequeue(out var lease))
 				{
-					lease.DisposeResources();
+					if (lease.Accelerator == accelerator)
+					{
+						lease.DisposeResources();
+					}
+					else
+					{
+						retained.Add(lease);
+					}
 				}
-				else
+
+				int retainedCount = retained.Count;
+				for (int i = 0; i < retainedCount; i++)
 				{
-					retained.Add(lease);
+					Pool.Enqueue(retained[i]);
 				}
+
+				WarmedLeaseCount = Pool.Count;
 			}
 
-			int retainedCount = retained.Count;
-			for (int i = 0; i < retainedCount; i++)
+			lock (SharedTablesLock)
 			{
-				Pool.Enqueue(retained[i]);
+				if (SharedTables.Remove(accelerator, out var tables))
+				{
+					tables.Dispose();
+				}
+
+				if (SharedTables.Count == 0)
+				{
+					var emptyTables = default(HeuristicGpuDivisorTables);
+					HeuristicGpuDivisorTables.InitializeShared(in emptyTables);
+				}
 			}
 		}
 
@@ -342,8 +379,11 @@ public sealed class PrimeTester
 
 			internal readonly Context AcceleratorContext;
 			public readonly Accelerator Accelerator;
-			public readonly Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> Kernel;
-			public readonly Action<Index1D, ArrayView<int>, ulong, ulong, HeuristicGpuDivisorTableKind, HeuristicGpuDivisorTables> HeuristicTrialDivisionKernel;
+			
+			public static readonly Action<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>> Kernel = SharedGpuContext.Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterKernels.SmallPrimeSieveKernel);
+
+			public static readonly Action<Index1D, ArrayView<int>, ulong, ulong, HeuristicGpuDivisorTableKind, HeuristicGpuDivisorTables> HeuristicTrialDivisionKernel = SharedGpuContext.Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<int>, ulong, ulong, HeuristicGpuDivisorTableKind, HeuristicGpuDivisorTables>(PrimeTesterKernels.HeuristicTrialDivisionKernel);
+			
 			public MemoryBuffer1D<ulong, Stride1D.Dense> Input = null!;
 			public MemoryBuffer1D<byte, Stride1D.Dense> Output = null!;
 			public MemoryBuffer1D<int, Stride1D.Dense> HeuristicFlag = null!;
@@ -380,30 +420,15 @@ public sealed class PrimeTester
 
 			internal PrimeTesterGpuContextLease(int minBufferCapacity)
 			{
-				AcceleratorContext = Context.CreateDefault();
-				Accelerator = AcceleratorContext.GetPreferredDevice(false).CreateAccelerator(AcceleratorContext);
-				Kernel = Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterKernels.SmallPrimeSieveKernel);
-                                HeuristicTrialDivisionKernel = Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<int>, ulong, ulong, HeuristicGpuDivisorTableKind, HeuristicGpuDivisorTables>(PrimeTesterKernels.HeuristicTrialDivisionKernel);
+				AcceleratorContext = SharedGpuContext.Context;
+				Accelerator = SharedGpuContext.Accelerator;
+				// GpuStaticTableInitializer.EnsureStaticTables(Accelerator);
+				// Kernel = Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<ulong>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<uint>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<ulong>, ArrayView<byte>>(PrimeTesterKernels.SmallPrimeSieveKernel);
+				// HeuristicTrialDivisionKernel = Accelerator.LoadAutoGroupedStreamKernel<Index1D, ArrayView<int>, ulong, ulong, HeuristicGpuDivisorTableKind, HeuristicGpuDivisorTables>(PrimeTesterKernels.HeuristicTrialDivisionKernel);
 
 				_sharedTables = RentSharedTables(Accelerator);
 
-				_heuristicDivisorTables = new HeuristicGpuDivisorTables(
-					_sharedTables.HeuristicCombinedDivisorsEnding1.View,
-					_sharedTables.HeuristicCombinedDivisorsEnding3.View,
-					_sharedTables.HeuristicCombinedDivisorsEnding7.View,
-					_sharedTables.HeuristicCombinedDivisorsEnding9.View,
-					_sharedTables.HeuristicCombinedDivisorSquaresEnding1.View,
-					_sharedTables.HeuristicCombinedDivisorSquaresEnding3.View,
-					_sharedTables.HeuristicCombinedDivisorSquaresEnding7.View,
-					_sharedTables.HeuristicCombinedDivisorSquaresEnding9.View,
-					_sharedTables.HeuristicGroupADivisors.View,
-					_sharedTables.HeuristicGroupADivisorSquares.View,
-					_sharedTables.HeuristicGroupBDivisorsEnding1.View,
-					_sharedTables.HeuristicGroupBDivisorSquaresEnding1.View,
-					_sharedTables.HeuristicGroupBDivisorsEnding7.View,
-					_sharedTables.HeuristicGroupBDivisorSquaresEnding7.View,
-					_sharedTables.HeuristicGroupBDivisorsEnding9.View,
-					_sharedTables.HeuristicGroupBDivisorSquaresEnding9.View);
+				_heuristicDivisorTables = _sharedTables.CreateHeuristicDivisorTables();
 				HeuristicGpuDivisorTables.InitializeShared(in _heuristicDivisorTables);
 
 				EnsureCapacity(minBufferCapacity);
@@ -439,7 +464,7 @@ public sealed class PrimeTester
 				}
 
 				_returned = true;
-				PrimeTesterGpuContextPool.Return(this);
+				Return(this);
 			}
 
 			internal void DisposeResources()
@@ -448,8 +473,24 @@ public sealed class PrimeTester
 				Output?.Dispose();
 				HeuristicFlag?.Dispose();
 				ReleaseSharedTables(_sharedTables);
-				Accelerator.Dispose();
-				AcceleratorContext.Dispose();
+				// These resources are shared between GPU leases
+				// Accelerator.Dispose();
+				// AcceleratorContext.Dispose();
+			}
+		}
+
+		private sealed class AcceleratorReferenceComparer : IEqualityComparer<Accelerator>
+		{
+			internal static AcceleratorReferenceComparer Instance { get; } = new();
+
+			public bool Equals(Accelerator? x, Accelerator? y)
+			{
+				return ReferenceEquals(x, y);
+			}
+
+			public int GetHashCode(Accelerator obj)
+			{
+				return RuntimeHelpers.GetHashCode(obj);
 			}
 		}
 
@@ -457,9 +498,8 @@ public sealed class PrimeTester
 		{
 			private int _referenceCount;
 
-			internal SharedHeuristicGpuTables(Accelerator accelerator, string key)
+			internal SharedHeuristicGpuTables(Accelerator accelerator)
 			{
-				Key = key;
 
 				var primesDefault = DivisorGenerator.SmallPrimes;
 				DevicePrimesDefault = accelerator.Allocate1D<uint>(primesDefault.Length);
@@ -540,7 +580,6 @@ public sealed class PrimeTester
 				HeuristicCombinedDivisorSquaresEnding9 = CopySpanToDevice(accelerator, combinedSquaresEnding9);
 			}
 
-			internal string Key { get; }
 			internal MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesDefault { get; }
 			internal MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastOne { get; }
 			internal MemoryBuffer1D<uint, Stride1D.Dense> DevicePrimesLastSeven { get; }
@@ -568,6 +607,27 @@ public sealed class PrimeTester
 			internal MemoryBuffer1D<ulong, Stride1D.Dense> HeuristicCombinedDivisorsEnding9 { get; }
 			internal MemoryBuffer1D<ulong, Stride1D.Dense> HeuristicCombinedDivisorSquaresEnding9 { get; }
 
+			internal HeuristicGpuDivisorTables CreateHeuristicDivisorTables()
+			{
+				return new HeuristicGpuDivisorTables(
+					HeuristicCombinedDivisorsEnding1.View,
+					HeuristicCombinedDivisorsEnding3.View,
+					HeuristicCombinedDivisorsEnding7.View,
+					HeuristicCombinedDivisorsEnding9.View,
+					HeuristicCombinedDivisorSquaresEnding1.View,
+					HeuristicCombinedDivisorSquaresEnding3.View,
+					HeuristicCombinedDivisorSquaresEnding7.View,
+					HeuristicCombinedDivisorSquaresEnding9.View,
+					HeuristicGroupADivisors.View,
+					HeuristicGroupADivisorSquares.View,
+					HeuristicGroupBDivisorsEnding1.View,
+					HeuristicGroupBDivisorSquaresEnding1.View,
+					HeuristicGroupBDivisorsEnding7.View,
+					HeuristicGroupBDivisorSquaresEnding7.View,
+					HeuristicGroupBDivisorsEnding9.View,
+					HeuristicGroupBDivisorSquaresEnding9.View);
+			}
+
 			internal void AddReference()
 			{
 				_referenceCount++;