Merge branch 'codex/optimize-evenperfectbitscanner-for-static-arrays' from github

PerfectNumbers.Core/Gpu/GpuContextPool.cs

@@ -1,4 +1,3 @@
-using System;
 using System.Collections.Concurrent;
 using ILGPU;
 using ILGPU.Runtime;
@@ -15,14 +14,11 @@ public static class GpuContextPool
 		public readonly Context Context;
 		public readonly Accelerator Accelerator;
 
-		private static readonly Context _sharedContext = Context.CreateDefault();
-		private static readonly Device _sharedDevice = _sharedContext.GetPreferredDevice(false);
-		private static readonly Accelerator _sharedAccelerator = _sharedDevice.CreateAccelerator(_sharedContext);
-
 		public PooledContext()
 		{
-			Context = _sharedContext;
-			Accelerator = _sharedAccelerator;
+			Context = SharedGpuContext.Context;
+			Accelerator = SharedGpuContext.Accelerator;
+			GpuStaticTableInitializer.EnsureStaticTables(Accelerator);
 			// 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
@@ -48,7 +44,7 @@ public static class GpuContextPool
 	private static int WarmedGpuContextCount;
 
 	// Allows callers to choose CPU/GPU per use-case, decoupled from ForceCpu.
-        	public static GpuContextLease Rent()
+	public static GpuContextLease Rent()
 	{
 		if (GpuPool.TryDequeue(out var gpu))
 		{

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;
@@ -121,6 +122,20 @@ public sealed class GpuKernelLease
         }
     }
 
+    public Action<AcceleratorStream, Index1D, ulong, ArrayView1D<ulong, Stride1D.Dense>, int, MontgomeryDivisorData, ArrayView1D<ulong, Stride1D.Dense>, ArrayView1D<byte, 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<byte, 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<byte, 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

PerfectNumbers.Core/Gpu/GpuKernelPool.cs

@@ -1,3 +1,4 @@
+using PerfectNumbers.Core;
 using System.Collections.Concurrent;
 using ILGPU;
 using ILGPU.Runtime;
@@ -37,6 +38,7 @@ public sealed class KernelContainer
     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<byte, Stride1D.Dense>>? SpecialMax;
 
     // Optional device buffer with small divisor cycles (<= 4M). Index = divisor, value = cycle length.
     public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallCycles;
@@ -46,6 +48,9 @@ public sealed class KernelContainer
     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<byte, Stride1D.Dense>? SpecialMaxResult;
     public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastOne;
     public MemoryBuffer1D<ulong, Stride1D.Dense>? SmallPrimesPow2LastOne;
     public MemoryBuffer1D<uint, Stride1D.Dense>? SmallPrimesLastSeven;
@@ -129,6 +134,29 @@ public readonly struct SmallPrimeFactorScratch
     }
 }
 
+public readonly struct SpecialMaxScratch
+{
+    public readonly MemoryBuffer1D<ulong, Stride1D.Dense> FactorValues;
+    public readonly MemoryBuffer1D<ulong, Stride1D.Dense> CandidateValues;
+    public readonly MemoryBuffer1D<byte, Stride1D.Dense> ResultSlot;
+
+    public SpecialMaxScratch(
+        MemoryBuffer1D<ulong, Stride1D.Dense> factorValues,
+        MemoryBuffer1D<ulong, Stride1D.Dense> candidateValues,
+        MemoryBuffer1D<byte, Stride1D.Dense> resultSlot)
+    {
+        FactorValues = factorValues;
+        CandidateValues = candidateValues;
+        ResultSlot = resultSlot;
+    }
+
+    public ArrayView1D<ulong, Stride1D.Dense> FactorsView => FactorValues.View;
+
+    public ArrayView1D<ulong, Stride1D.Dense> CandidatesView => CandidateValues.View;
+
+    public ArrayView1D<byte, 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.
@@ -138,6 +166,13 @@ public class GpuKernelPool
         return KernelCache.GetOrAdd(accelerator, _ => new KernelContainer());
     }
 
+    internal static void PreloadStaticTables(Accelerator accelerator)
+    {
+        EnsureSmallCyclesOnDevice(accelerator);
+        _ = EnsureSmallPrimesOnDevice(accelerator);
+        _ = EnsureSmallPrimeFactorTables(accelerator);
+    }
+
     // 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)
@@ -277,6 +312,43 @@ public class GpuKernelPool
         }
     }
 
+    public static SpecialMaxScratch EnsureSpecialMaxScratch(Accelerator accelerator, int factorCapacity)
+    {
+        if (factorCapacity <= 0)
+        {
+            factorCapacity = 1;
+        }
+
+        var kernels = GetKernels(accelerator);
+        lock (kernels)
+        {
+            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<byte, Stride1D.Dense>? resultSlot = kernels.SpecialMaxResult;
+            if (resultSlot is null)
+            {
+                resultSlot = accelerator.Allocate1D<byte>(1);
+                kernels.SpecialMaxResult = resultSlot;
+            }
+
+            return new SpecialMaxScratch(factorValues!, candidateValues!, resultSlot!);
+        }
+    }
+
     public static GpuKernelLease GetKernel(bool useGpuOrder)
     {
         var limiter = GpuPrimeWorkLimiter.Acquire();

PerfectNumbers.Core/Gpu/GpuStaticTableInitializer.cs

@@ -0,0 +1,44 @@
+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(Accelerator accelerator)
+    {
+        if (accelerator is null)
+        {
+            return;
+        }
+
+        Initialized.GetOrAdd(accelerator, static acc =>
+        {
+            PrimeTester.PrimeTesterGpuContextPool.EnsureStaticTables(acc);
+            GpuKernelPool.PreloadStaticTables(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<byte, 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/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/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,10 +1,8 @@
-using System;
 using System.Buffers;
 using System.Diagnostics;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
-using ILGPU;
 using PerfectNumbers.Core.Gpu;
 
 namespace PerfectNumbers.Core;
@@ -108,7 +106,6 @@ internal static partial class PrimeOrderCalculator
 		return candidateOrder;
 	}
 
-
 	private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
 	{
 		int length = factors.Count;
@@ -121,16 +118,14 @@ 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;
+		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 +1189,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 +1226,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 +1696,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 +1712,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 +1749,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,112 @@ 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(accelerator);
+		SmallPrimeFactorScratch scratch = GpuKernelPool.EnsureSmallPrimeFactorScratch(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(accelerator, 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();
+
+		byte 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;
+		var limiter = GpuPrimeWorkLimiter.Acquire();
+		var gpu = PrimeTesterGpuContextPool.Rent(1);
+		var accelerator = gpu.Accelerator;
+		var kernel = gpu.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();
+		limiter.Dispose();
 
-		// 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;
@@ -221,23 +214,16 @@ public sealed class PrimeTester
 		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();
@@ -258,6 +244,19 @@ public sealed class PrimeTester
 			}
 		}
 
+		internal static void EnsureStaticTables(Accelerator accelerator)
+		{
+			if (accelerator is null)
+			{
+				return;
+			}
+
+			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)
@@ -272,6 +271,9 @@ public sealed class PrimeTester
 					return;
 				}
 
+				Accelerator accelerator = SharedGpuContext.Accelerator;
+				GpuStaticTableInitializer.EnsureStaticTables(accelerator);
+
 				int toCreate = target - WarmedLeaseCount;
 				for (int i = 0; i < toCreate; i++)
 				{
@@ -319,6 +321,8 @@ public sealed class PrimeTester
 				}
 
 				SharedTables.Clear();
+				var emptyTables = default(HeuristicGpuDivisorTables);
+				HeuristicGpuDivisorTables.InitializeShared(in emptyTables);
 			}
 
 			WarmedLeaseCount = 0;
@@ -354,6 +358,20 @@ public sealed class PrimeTester
 
 				WarmedLeaseCount = Pool.Count;
 			}
+
+			lock (SharedTablesLock)
+			{
+				if (SharedTables.Remove(accelerator, out var tables))
+				{
+					tables.Dispose();
+				}
+
+				if (SharedTables.Count == 0)
+				{
+					var emptyTables = default(HeuristicGpuDivisorTables);
+					HeuristicGpuDivisorTables.InitializeShared(in emptyTables);
+				}
+			}
 		}
 
 		internal sealed class PrimeTesterGpuContextLease : IDisposable
@@ -401,36 +419,17 @@ public sealed class PrimeTester
 
 			public HeuristicGpuDivisorTables HeuristicGpuTables => _heuristicDivisorTables;
 
-			private static readonly Context _sharedContext = Context.CreateDefault();
-			private static readonly Device _sharedDevice = _sharedContext.GetPreferredDevice(false);
-			private static readonly Accelerator _sharedAccelerator = _sharedDevice.CreateAccelerator(_sharedContext);
-
 			internal PrimeTesterGpuContextLease(int minBufferCapacity)
 			{
-				AcceleratorContext = _sharedContext;// Context.CreateDefault();
-				Accelerator = _sharedAccelerator; //AcceleratorContext.GetPreferredDevice(false).CreateAccelerator(AcceleratorContext);
+				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);
@@ -481,13 +480,27 @@ public sealed class PrimeTester
 			}
 		}
 
+		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);
+			}
+		}
+
 		private sealed class SharedHeuristicGpuTables
 		{
 			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);
@@ -568,7 +581,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; }
@@ -596,6 +608,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++;