From 4c79c4087e1acb0d8dc363f8784de14c0d87da52 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Wed, 29 Oct 2025 14:00:38 +0100 Subject: [PATCH] Various optimizations --- EvenPerfectBitScanner.Tests/ProgramTests.cs | 2 +- EvenPerfectBitScanner/Program.cs | 2 +- .../PrimeTester/PrimeTesterIsPrimeTests.cs | 2 +- .../Gpu/GpuPrimeWorkLimiter.cs | 29 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 8 +- .../Gpu/MersenneNumberResidueGpuTester.cs | 4 +- .../HeuristicCombinedPrimeTester.cs | 1334 ++++++++--------- PerfectNumbers.Core/HeuristicPrimeTester.cs | 14 +- PerfectNumbers.Core/PrimeTester.cs | 1013 ++++++------- 9 files changed, 1179 insertions(+), 1229 deletions(-) diff --git a/EvenPerfectBitScanner.Tests/ProgramTests.cs b/EvenPerfectBitScanner.Tests/ProgramTests.cs index 8b1a8fa..f08a676 100644 --- a/EvenPerfectBitScanner.Tests/ProgramTests.cs +++ b/EvenPerfectBitScanner.Tests/ProgramTests.cs @@ -55,7 +55,7 @@ public class ProgramTests [InlineData(81UL, false)] public void IsPrime_identifies_primes_correctly(ulong n, bool expected) { - new PrimeTester().IsPrime(n, CancellationToken.None).Should().Be(expected); + PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected); } [Fact] diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 48ff15b..86b4c20 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -622,7 +622,7 @@ internal static class Program if (!_cliArguments.UseByDivisor) { if (!(GpuContextPool.ForceCpu - ? PrimeTesters.Value!.IsPrime(p, CancellationToken.None) + ? PrimeTester.IsPrime(p, CancellationToken.None) : PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None))) { _lastCompositeP = true; diff --git a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs index c99d670..e86e1db 100644 --- a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs @@ -26,7 +26,7 @@ public class PrimeTesterIsPrimeTests { var tester = new PrimeTester(); - tester.IsPrime(n, CancellationToken.None).Should().Be(expected); + PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected); } [Theory] diff --git a/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs b/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs index 97c08ba..1d39b7d 100644 --- a/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs +++ b/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs @@ -3,7 +3,6 @@ namespace PerfectNumbers.Core.Gpu; public static class GpuPrimeWorkLimiter { private static SemaphoreSlim _semaphore = new(1, int.MaxValue); - private static readonly object _lock = new(); private static int _currentLimit = 1; public static Lease Acquire() @@ -11,29 +10,21 @@ public static class GpuPrimeWorkLimiter // 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(); + // sem.Wait(); return new Lease(sem); } public static void SetLimit(int value) { - if (value < 1) - { - value = 1; - } + if (value == _currentLimit) + { + return; + } - lock (_lock) - { - 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; - } + // 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 @@ -47,7 +38,7 @@ public static class GpuPrimeWorkLimiter public void Dispose() { - _sem.Release(); + // _sem.Release(); } } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index f480afd..2c9ca8b 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -94,7 +94,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv var gpuLease = RentAccelerator(); var accelerator = gpuLease.Accelerator; - Monitor.Enter(gpuLease.ExecutionLock); + // Monitor.Enter(gpuLease.ExecutionLock); var kernel = GetKernel(accelerator); var resources = RentBatchResources(accelerator, batchCapacity); @@ -121,7 +121,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv out processedCount); ReturnBatchResources(resources); - Monitor.Exit(gpuLease.ExecutionLock); + // Monitor.Exit(gpuLease.ExecutionLock); ReturnAccelerator(gpuLease); if (composite) @@ -560,7 +560,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } - Monitor.Enter(_lease.ExecutionLock); + // Monitor.Enter(_lease.ExecutionLock); EnsureExecutionResourcesLocked(length); @@ -601,7 +601,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ref byte hitRef = ref MemoryMarshal.GetReference(hitSlice); hitView.CopyToCPU(ref hitRef, length); - Monitor.Exit(_lease.ExecutionLock); + // Monitor.Exit(_lease.ExecutionLock); } public void Dispose() diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index ca42ddd..759d575 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -30,7 +30,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) var gpuLease = RentAccelerator(); var accelerator = gpuLease.Accelerator; - Monitor.Enter(gpuLease.ExecutionLock); + // Monitor.Enter(gpuLease.ExecutionLock); var stream = accelerator.CreateStream(); var kernel = GetPow2ModKernel(accelerator); @@ -108,7 +108,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) ReturnResources(accelerator, resources); stream.Dispose(); - Monitor.Exit(gpuLease.ExecutionLock); + // Monitor.Exit(gpuLease.ExecutionLock); ReturnAccelerator(gpuLease); limiter.Dispose(); } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 1ae2880..b9ec773 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -13,682 +13,660 @@ namespace PerfectNumbers.Core; public sealed class HeuristicCombinedPrimeTester { - static HeuristicCombinedPrimeTester() - { - HeuristicPrimeSieves.EnsureInitialized(); - - CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB); - - CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB); - - CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB); - - CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB; - CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB; - CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB; - CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB; - } - - private const ulong Wheel210 = 210UL; - - private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. - - private static readonly uint[] CombinedDivisorsEnding1TwoAOneB; - private static readonly uint[] CombinedDivisorsEnding3TwoAOneB; - private static readonly uint[] CombinedDivisorsEnding7TwoAOneB; - private static readonly uint[] CombinedDivisorsEnding9TwoAOneB; - - private static readonly uint[] CombinedDivisorsEnding1OneAOneB; - private static readonly uint[] CombinedDivisorsEnding3OneAOneB; - private static readonly uint[] CombinedDivisorsEnding7OneAOneB; - private static readonly uint[] CombinedDivisorsEnding9OneAOneB; - - private static readonly uint[] CombinedDivisorsEnding1ThreeAOneB; - private static readonly uint[] CombinedDivisorsEnding3ThreeAOneB; - private static readonly uint[] CombinedDivisorsEnding7ThreeAOneB; - private static readonly uint[] CombinedDivisorsEnding9ThreeAOneB; - - private static readonly uint[] CombinedDivisorsEnding1ThreeATwoB; - private static readonly uint[] CombinedDivisorsEnding3ThreeATwoB; - private static readonly uint[] CombinedDivisorsEnding7ThreeATwoB; - private static readonly uint[] CombinedDivisorsEnding9ThreeATwoB; - - private static readonly uint[] CombinedDivisorsEnding1; - private static readonly uint[] CombinedDivisorsEnding3; - private static readonly uint[] CombinedDivisorsEnding7; - private static readonly uint[] CombinedDivisorsEnding9; - - private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - - internal enum CombinedDivisorPattern : byte - { - TwoAOneB = 0, - OneAOneB = 1, - ThreeAOneB = 2, - ThreeATwoB = 3, - } - - internal enum HeuristicDivisorGroup : byte - { - None = 0, - GroupAConstant = 1, - GroupAWheel = 2, - GroupB = 3, - } - - internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue) - { - public readonly ulong Value = value; - public readonly HeuristicDivisorGroup Group = group; - public readonly byte Ending = ending; - public readonly byte PriorityIndex = priorityIndex; - public readonly ushort WheelResidue = wheelResidue; - } - - internal readonly struct HeuristicDivisorPreparation( - in HeuristicDivisorCandidate candidate, - in MontgomeryDivisorData divisorData, - ulong cycleLengthHint, - bool hasCycleLengthHint) - { - public readonly HeuristicDivisorCandidate Candidate = candidate; - public readonly MontgomeryDivisorData DivisorData = divisorData; - public readonly ulong CycleLengthHint = cycleLengthHint; - public readonly bool HasCycleLengthHint = hasCycleLengthHint; - public readonly bool RequiresCycleComputation = !hasCycleLengthHint; - } - - - public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096; - - public static int HeuristicDivisorInterleaveBatchSize { get; set; } = 64; - - [ThreadStatic] - private static HeuristicCombinedPrimeTester? _tester; - - public static HeuristicCombinedPrimeTester Exclusive => _tester ??= new(); - - public bool IsPrimeCpu(ulong n, CancellationToken ct) - { - ct.ThrowIfCancellationRequested(); - - byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); - - if (sqrtLimit < 3UL) - { - return EvaluateWithOpenNumericFallback(n); - } - - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); - } - - public bool IsPrimeCpu(ulong n) - { - return IsPrimeCpu(n, CancellationToken.None); - } - - public bool IsPrimeGpu(ulong n) - { - byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); - return IsPrimeGpu(n, sqrtLimit, nMod10); - } - - public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) - { - return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10); - } - - private bool HeuristicIsPrimeGpuCore(ulong n, ulong sqrtLimit, byte nMod10) - { - // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? - if (sqrtLimit < 3UL) - { - return EvaluateWithOpenNumericFallback(n); - } - - if (!UseHeuristicGroupBTrialDivision) - { - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); - } - - bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); - return !compositeDetected; - } - - private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10) - { - ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); - if (combinedDivisors.IsEmpty) - { - return EvaluateWithOpenNumericFallback(n); - } - - int length = combinedDivisors.Length; - for (int i = 0; i < length; i++) - { - uint entry = combinedDivisors[i]; - if (entry > sqrtLimit) - { - break; - } - - if (n % entry == 0UL) - { - return false; - } - } - - return EvaluateWithOpenNumericFallback(n); - } - - private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) - { - int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); - var divisorPool = ThreadStaticPools.UlongPool; - var hitPool = ThreadStaticPools.BytePool; - - ulong[] divisorArray = divisorPool.Rent(batchCapacity); - byte[] hitFlags = hitPool.Rent(batchCapacity); - - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(); - var accelerator = gpu.Accelerator; - var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); - var scratch = state.RentScratch(batchCapacity, accelerator); - - bool compositeDetected = false; - int count = 0; - - lock (gpu.ExecutionLock) - { - var input = scratch.Input; - var output = scratch.Output; - - bool ProcessBatch(int length) - { - if (length <= 0) - { - return false; - } - - input.View.CopyFromCPU(ref divisorArray[0], length); - state.HeuristicTrialDivisionKernel(length, input.View, n, output.View); - accelerator.Synchronize(); - output.View.CopyToCPU(ref hitFlags[0], length); - - for (int i = 0; i < length; i++) - { - if (hitFlags[i] == 0) - { - continue; - } - - ulong divisor = divisorArray[i]; - if (divisor > 1UL && n % divisor == 0UL) - { - return true; - } - } - - return false; - } - - ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); - if (combinedDivisors.IsEmpty) - { - goto Cleanup; - } - - int length = combinedDivisors.Length; - for (int i = 0; i < length; i++) - { - uint entry = combinedDivisors[i]; - if (entry > sqrtLimit) - { - break; - } - - divisorArray[count] = entry; - count++; - - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; - } - } - - if (!compositeDetected && count > 0) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; - } - } - -Cleanup: - state.ReturnScratch(scratch); - divisorPool.Return(divisorArray); - hitPool.Return(hitFlags); - gpu.Dispose(); - limiter.Dispose(); - return compositeDetected; - } - - private static bool EvaluateWithOpenNumericFallback(ulong n) - { - return Prime.Numbers.IsPrime(n); - } - - internal static ulong ComputeHeuristicSqrt(ulong n) - { - ulong sqrt = (ulong)Math.Sqrt(n); - UInt128 square = (UInt128)sqrt * sqrt; - - while (square > n) - { - sqrt--; - square = (UInt128)sqrt * sqrt; - } - - ulong next = sqrt + 1UL; - UInt128 nextSquare = (UInt128)next * next; - while (nextSquare <= n) - { - sqrt = next; - next++; - nextSquare = (UInt128)next * next; - } - - return sqrt; - } - - private static ReadOnlySpan GetCombinedDivisors(byte nMod10) - { - return nMod10 switch - { - 1 => CombinedDivisorsEnding1, - 3 => CombinedDivisorsEnding3, - 7 => CombinedDivisorsEnding7, - 9 => CombinedDivisorsEnding9, - _ => ReadOnlySpan.Empty, - }; - } - - internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) - { - return pattern switch - { - CombinedDivisorPattern.TwoAOneB => nMod10 switch - { - 1 => CombinedDivisorsEnding1TwoAOneB, - 3 => CombinedDivisorsEnding3TwoAOneB, - 7 => CombinedDivisorsEnding7TwoAOneB, - 9 => CombinedDivisorsEnding9TwoAOneB, - _ => ReadOnlySpan.Empty, - }, - CombinedDivisorPattern.OneAOneB => nMod10 switch - { - 1 => CombinedDivisorsEnding1OneAOneB, - 3 => CombinedDivisorsEnding3OneAOneB, - 7 => CombinedDivisorsEnding7OneAOneB, - 9 => CombinedDivisorsEnding9OneAOneB, - _ => ReadOnlySpan.Empty, - }, - CombinedDivisorPattern.ThreeAOneB => nMod10 switch - { - 1 => CombinedDivisorsEnding1ThreeAOneB, - 3 => CombinedDivisorsEnding3ThreeAOneB, - 7 => CombinedDivisorsEnding7ThreeAOneB, - 9 => CombinedDivisorsEnding9ThreeAOneB, - _ => ReadOnlySpan.Empty, - }, - CombinedDivisorPattern.ThreeATwoB => nMod10 switch - { - 1 => CombinedDivisorsEnding1ThreeATwoB, - 3 => CombinedDivisorsEnding3ThreeATwoB, - 7 => CombinedDivisorsEnding7ThreeATwoB, - 9 => CombinedDivisorsEnding9ThreeATwoB, - _ => ReadOnlySpan.Empty, - }, - _ => ReadOnlySpan.Empty, - }; - } - - private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) - { - ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; - ReadOnlySpan groupB = GetGroupBDivisors(nMod10); - - var combined = new List(groupA.Length + groupB.Length); - - int groupAIndex = 0; - while (groupAIndex < groupA.Length && groupA[groupAIndex] <= 19) - { - combined.Add((uint)groupA[groupAIndex]); - groupAIndex++; - } - - int groupBIndex = GetGroupBStartIndex(groupB); - int groupAInterleaveCount; - int groupBInterleaveCount; - switch (pattern) - { - case CombinedDivisorPattern.OneAOneB: - groupAInterleaveCount = 1; - groupBInterleaveCount = 1; - break; - case CombinedDivisorPattern.TwoAOneB: - groupAInterleaveCount = 2; - groupBInterleaveCount = 1; - break; - case CombinedDivisorPattern.ThreeAOneB: - groupAInterleaveCount = 3; - groupBInterleaveCount = 1; - break; - case CombinedDivisorPattern.ThreeATwoB: - groupAInterleaveCount = 3; - groupBInterleaveCount = 2; - break; - default: - groupAInterleaveCount = 2; - groupBInterleaveCount = 1; - break; - } - - while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) - { - int addedA = 0; - while (addedA < groupAInterleaveCount && groupAIndex < groupA.Length) - { - combined.Add((uint)groupA[groupAIndex]); - groupAIndex++; - addedA++; - } - - int addedB = 0; - while (addedB < groupBInterleaveCount && groupBIndex < groupB.Length) - { - combined.Add(groupB[groupBIndex]); - groupBIndex++; - addedB++; - } - } - - return combined.ToArray(); - } - - - - private static int GetGroupBStartIndex(ReadOnlySpan divisors) - { - int index = 0; - while (index < divisors.Length && IsGroupAPrefixValue(divisors[index])) - { - index++; - } - - return index; - } - - private static bool IsGroupAPrefixValue(uint value) => value is 3U or 7U or 11U or 13U; - - private static HeuristicDivisorGroup ResolveGroup(uint value) => value switch - { - 3U or 7U or 11U or 13U => HeuristicDivisorGroup.GroupAConstant, - _ when value % 10U == 3U => HeuristicDivisorGroup.GroupAWheel, - _ => HeuristicDivisorGroup.GroupB, - }; - - private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch - { - 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, - 3 => ReadOnlySpan.Empty, - 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, - 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, - _ => ReadOnlySpan.Empty, - }; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) - { - return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); - } - - internal struct HeuristicGroupBSequenceState - { - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor) - { - return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) - { - ulong divisor = candidate.Value; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); - return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong ResolveHeuristicCycleLength( - ulong exponent, - in HeuristicDivisorPreparation preparation, - out bool cycleFromHint, - out bool cycleComputed, - out bool primeOrderFailed) - { - if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL) - { - cycleFromHint = true; - cycleComputed = true; - primeOrderFailed = false; - return preparation.CycleLengthHint; - } - - ulong divisor = preparation.Candidate.Value; - MontgomeryDivisorData divisorData = preparation.DivisorData; - - bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( - divisor, - exponent, - divisorData, - out ulong computedCycle, - out bool primeOrderFailedLocal); - - if (trySuccess && computedCycle != 0UL) - { - cycleFromHint = false; - cycleComputed = true; - primeOrderFailed = primeOrderFailedLocal; - return computedCycle; - } - - primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL; - ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); - cycleFromHint = false; - cycleComputed = resolvedCycle != 0UL; - return resolvedCycle; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength) - { - ulong[] snapshot = HeuristicSmallCycleSnapshot; - if (divisor < (ulong)snapshot.Length) - { - cycleLength = snapshot[(int)divisor]; - if (cycleLength != 0UL) - { - return true; - } - } - - cycleLength = 0UL; - return false; - } - - internal ref struct HeuristicDivisorEnumerator - { - private readonly ulong sqrtLimit; - private readonly ReadOnlySpan combinedDivisors; - private int index; - - public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) - { - this.sqrtLimit = sqrtLimit; - _ = groupBBuffer; - combinedDivisors = GetCombinedDivisors(nMod10); - index = 0; - } - - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - while (index < combinedDivisors.Length) - { - uint entry = combinedDivisors[index++]; - if (entry > sqrtLimit) - { - index = combinedDivisors.Length; - break; - } - - HeuristicDivisorGroup group = ResolveGroup(entry); - ulong value = entry; - byte ending = (byte)(value % 10UL); - ushort residue = (ushort)(value % Wheel210); - candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); - return true; - } - - candidate = default; - return false; - } - } - - internal ref struct MersenneHeuristicDivisorEnumerator - { - private readonly GpuUInt128 step; - private readonly GpuUInt128 limit; - private GpuUInt128 current; - private MersenneDivisorResidueStepper residueStepper; - private bool active; - private ulong processedCount; - private ulong lastDivisor; - - public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor) - { - var stepLocal = new GpuUInt128(exponent); - stepLocal.ShiftLeft(1); - step = stepLocal; - - var limitLocal = new GpuUInt128(maxDivisor); - limit = limitLocal; - - var firstDivisor = stepLocal; - firstDivisor.Add(1UL); - - bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0; - - current = hasCandidates ? firstDivisor : GpuUInt128.Zero; - residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default; - active = hasCandidates; - processedCount = 0UL; - lastDivisor = 0UL; - } - - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - while (active) - { - ulong value = current.Low; - processedCount++; - lastDivisor = value; - - bool admissible = residueStepper.IsAdmissible(); - - Advance(); - - if (admissible) - { - candidate = CreateCandidate(value); - return true; - } - } - - candidate = default; - return false; - } - - private void Advance() - { - GpuUInt128 next = current; - next.Add(step); - - if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0) - { - active = false; - return; - } - - current = next; - residueStepper.Advance(); - } - - public ulong ProcessedCount => processedCount; - - public ulong LastDivisor => lastDivisor; - - public bool Exhausted => !active; - - private static HeuristicDivisorCandidate CreateCandidate(ulong value) - { - byte ending = (byte)(value % 10UL); - HeuristicDivisorGroup group; - - if (value == 3UL || value == 7UL || value == 11UL || value == 13UL) - { - group = HeuristicDivisorGroup.GroupAConstant; - } - else if (ending == 3) - { - group = HeuristicDivisorGroup.GroupAWheel; - } - else - { - group = HeuristicDivisorGroup.GroupB; - } - - ushort residue = (ushort)(value % Wheel210); - return new HeuristicDivisorCandidate(value, group, ending, 0, residue); - } - } + static HeuristicCombinedPrimeTester() + { + HeuristicPrimeSieves.EnsureInitialized(); + + CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB); + + CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB); + + CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB); + + CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB; + CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB; + CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB; + CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB; + } + + private const ulong Wheel210 = 210UL; + + private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. + + private static readonly uint[] CombinedDivisorsEnding1TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding3TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding7TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding9TwoAOneB; + + private static readonly uint[] CombinedDivisorsEnding1OneAOneB; + private static readonly uint[] CombinedDivisorsEnding3OneAOneB; + private static readonly uint[] CombinedDivisorsEnding7OneAOneB; + private static readonly uint[] CombinedDivisorsEnding9OneAOneB; + + private static readonly uint[] CombinedDivisorsEnding1ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding3ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding7ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding9ThreeAOneB; + + private static readonly uint[] CombinedDivisorsEnding1ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding3ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding7ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding9ThreeATwoB; + + private static readonly uint[] CombinedDivisorsEnding1; + private static readonly uint[] CombinedDivisorsEnding3; + private static readonly uint[] CombinedDivisorsEnding7; + private static readonly uint[] CombinedDivisorsEnding9; + + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + + internal enum CombinedDivisorPattern : byte + { + TwoAOneB = 0, + OneAOneB = 1, + ThreeAOneB = 2, + ThreeATwoB = 3, + } + + internal enum HeuristicDivisorGroup : byte + { + None = 0, + GroupAConstant = 1, + GroupAWheel = 2, + GroupB = 3, + } + + internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue) + { + public readonly ulong Value = value; + public readonly HeuristicDivisorGroup Group = group; + public readonly byte Ending = ending; + public readonly byte PriorityIndex = priorityIndex; + public readonly ushort WheelResidue = wheelResidue; + } + + internal readonly struct HeuristicDivisorPreparation( + in HeuristicDivisorCandidate candidate, + in MontgomeryDivisorData divisorData, + ulong cycleLengthHint, + bool hasCycleLengthHint) + { + public readonly HeuristicDivisorCandidate Candidate = candidate; + public readonly MontgomeryDivisorData DivisorData = divisorData; + public readonly ulong CycleLengthHint = cycleLengthHint; + public readonly bool HasCycleLengthHint = hasCycleLengthHint; + public readonly bool RequiresCycleComputation = !hasCycleLengthHint; + } + + + internal static int HeuristicGpuDivisorBatchSize = 4_096; + internal static int HeuristicDivisorInterleaveBatchSize = 64; + + [ThreadStatic] + private static HeuristicCombinedPrimeTester? _tester; + + public static HeuristicCombinedPrimeTester Exclusive => _tester ??= new(); + + public static bool IsPrimeCpu(ulong n) + { + byte nMod10 = (byte)n.Mod10(); + ulong sqrtLimit = ComputeHeuristicSqrt(n); + + if (sqrtLimit < 3UL) + { + return EvaluateWithOpenNumericFallback(n); + } + + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); + } + + public static bool IsPrimeGpu(ulong n) + { + byte nMod10 = (byte)n.Mod10(); + ulong sqrtLimit = ComputeHeuristicSqrt(n); + return IsPrimeGpu(n, sqrtLimit, nMod10); + } + + public static bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) + { + // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? + if (sqrtLimit < 3UL) + { + return EvaluateWithOpenNumericFallback(n); + } + + bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); + return !compositeDetected; + } + + private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10) + { + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + if (combinedDivisors.IsEmpty) + { + return EvaluateWithOpenNumericFallback(n); + } + + int length = combinedDivisors.Length; + for (int i = 0; i < length; i++) + { + uint entry = combinedDivisors[i]; + if (entry > sqrtLimit) + { + break; + } + + if (n % entry == 0UL) + { + return false; + } + } + + return EvaluateWithOpenNumericFallback(n); + } + + private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) + { + int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); + var divisorPool = ThreadStaticPools.UlongPool; + var hitPool = ThreadStaticPools.BytePool; + + ulong[] divisorArray = divisorPool.Rent(batchCapacity); + byte[] hitFlags = hitPool.Rent(batchCapacity); + + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(); + var accelerator = gpu.Accelerator; + var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); + var scratch = state.RentScratch(batchCapacity, accelerator); + + bool compositeDetected = false; + int count = 0; + + var input = scratch.Input; + var output = scratch.Output; + + bool ProcessBatch(int length) + { + if (length <= 0) + { + return false; + } + + input.View.CopyFromCPU(ref divisorArray[0], length); + state.HeuristicTrialDivisionKernel(length, input.View, n, output.View); + accelerator.Synchronize(); + output.View.CopyToCPU(ref hitFlags[0], length); + + for (int i = 0; i < length; i++) + { + if (hitFlags[i] == 0) + { + continue; + } + + return true; + // ulong divisor = divisorArray[i]; + // if (divisor > 1UL && n % divisor == 0UL) + // { + // return true; + // } + } + + return false; + } + + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + if (combinedDivisors.IsEmpty) + { + goto Cleanup; + } + + int length = combinedDivisors.Length; + for (int i = 0; i < length; i++) + { + uint entry = combinedDivisors[i]; + if (entry > sqrtLimit) + { + break; + } + + divisorArray[count] = entry; + count++; + + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + + count = 0; + } + } + + if (!compositeDetected && count > 0) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + } + + Cleanup: + state.ReturnScratch(scratch); + divisorPool.Return(divisorArray); + hitPool.Return(hitFlags); + gpu.Dispose(); + limiter.Dispose(); + return compositeDetected; + } + + private static bool EvaluateWithOpenNumericFallback(ulong n) + { + return Prime.Numbers.IsPrime(n); + } + + internal static ulong ComputeHeuristicSqrt(ulong n) + { + ulong sqrt = (ulong)Math.Sqrt(n); + UInt128 square = (UInt128)sqrt * sqrt; + + while (square > n) + { + sqrt--; + square = (UInt128)sqrt * sqrt; + } + + ulong next = sqrt + 1UL; + UInt128 nextSquare = (UInt128)next * next; + while (nextSquare <= n) + { + sqrt = next; + next++; + nextSquare = (UInt128)next * next; + } + + return sqrt; + } + + private static ReadOnlySpan GetCombinedDivisors(byte nMod10) + { + return nMod10 switch + { + 1 => CombinedDivisorsEnding1, + 3 => CombinedDivisorsEnding3, + 7 => CombinedDivisorsEnding7, + 9 => CombinedDivisorsEnding9, + _ => ReadOnlySpan.Empty, + }; + } + + internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) + { + return pattern switch + { + CombinedDivisorPattern.TwoAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1TwoAOneB, + 3 => CombinedDivisorsEnding3TwoAOneB, + 7 => CombinedDivisorsEnding7TwoAOneB, + 9 => CombinedDivisorsEnding9TwoAOneB, + _ => ReadOnlySpan.Empty, + }, + CombinedDivisorPattern.OneAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1OneAOneB, + 3 => CombinedDivisorsEnding3OneAOneB, + 7 => CombinedDivisorsEnding7OneAOneB, + 9 => CombinedDivisorsEnding9OneAOneB, + _ => ReadOnlySpan.Empty, + }, + CombinedDivisorPattern.ThreeAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeAOneB, + 3 => CombinedDivisorsEnding3ThreeAOneB, + 7 => CombinedDivisorsEnding7ThreeAOneB, + 9 => CombinedDivisorsEnding9ThreeAOneB, + _ => ReadOnlySpan.Empty, + }, + CombinedDivisorPattern.ThreeATwoB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeATwoB, + 3 => CombinedDivisorsEnding3ThreeATwoB, + 7 => CombinedDivisorsEnding7ThreeATwoB, + 9 => CombinedDivisorsEnding9ThreeATwoB, + _ => ReadOnlySpan.Empty, + }, + _ => ReadOnlySpan.Empty, + }; + } + + private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) + { + ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; + ReadOnlySpan groupB = GetGroupBDivisors(nMod10); + + var combined = new List(groupA.Length + groupB.Length); + + int groupAIndex = 0; + while (groupAIndex < groupA.Length && groupA[groupAIndex] <= 19) + { + combined.Add((uint)groupA[groupAIndex]); + groupAIndex++; + } + + int groupBIndex = GetGroupBStartIndex(groupB); + int groupAInterleaveCount; + int groupBInterleaveCount; + switch (pattern) + { + case CombinedDivisorPattern.OneAOneB: + groupAInterleaveCount = 1; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.TwoAOneB: + groupAInterleaveCount = 2; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.ThreeAOneB: + groupAInterleaveCount = 3; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.ThreeATwoB: + groupAInterleaveCount = 3; + groupBInterleaveCount = 2; + break; + default: + groupAInterleaveCount = 2; + groupBInterleaveCount = 1; + break; + } + + while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) + { + int addedA = 0; + while (addedA < groupAInterleaveCount && groupAIndex < groupA.Length) + { + combined.Add((uint)groupA[groupAIndex]); + groupAIndex++; + addedA++; + } + + int addedB = 0; + while (addedB < groupBInterleaveCount && groupBIndex < groupB.Length) + { + combined.Add(groupB[groupBIndex]); + groupBIndex++; + addedB++; + } + } + + return combined.ToArray(); + } + + + + private static int GetGroupBStartIndex(ReadOnlySpan divisors) + { + int index = 0; + while (index < divisors.Length && IsGroupAPrefixValue(divisors[index])) + { + index++; + } + + return index; + } + + private static bool IsGroupAPrefixValue(uint value) => value is 3U or 7U or 11U or 13U; + + private static HeuristicDivisorGroup ResolveGroup(uint value) => value switch + { + 3U or 7U or 11U or 13U => HeuristicDivisorGroup.GroupAConstant, + _ when value % 10U == 3U => HeuristicDivisorGroup.GroupAWheel, + _ => HeuristicDivisorGroup.GroupB, + }; + + private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch + { + 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, + 3 => ReadOnlySpan.Empty, + 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, + 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, + _ => ReadOnlySpan.Empty, + }; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); + } + + internal struct HeuristicGroupBSequenceState + { + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) + { + ulong divisor = candidate.Value; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); + return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong ResolveHeuristicCycleLength( + ulong exponent, + in HeuristicDivisorPreparation preparation, + out bool cycleFromHint, + out bool cycleComputed, + out bool primeOrderFailed) + { + if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL) + { + cycleFromHint = true; + cycleComputed = true; + primeOrderFailed = false; + return preparation.CycleLengthHint; + } + + ulong divisor = preparation.Candidate.Value; + MontgomeryDivisorData divisorData = preparation.DivisorData; + + bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + divisor, + exponent, + divisorData, + out ulong computedCycle, + out bool primeOrderFailedLocal); + + if (trySuccess && computedCycle != 0UL) + { + cycleFromHint = false; + cycleComputed = true; + primeOrderFailed = primeOrderFailedLocal; + return computedCycle; + } + + primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL; + ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); + cycleFromHint = false; + cycleComputed = resolvedCycle != 0UL; + return resolvedCycle; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength) + { + ulong[] snapshot = HeuristicSmallCycleSnapshot; + if (divisor < (ulong)snapshot.Length) + { + cycleLength = snapshot[(int)divisor]; + if (cycleLength != 0UL) + { + return true; + } + } + + cycleLength = 0UL; + return false; + } + + internal ref struct HeuristicDivisorEnumerator + { + private readonly ulong sqrtLimit; + private readonly ReadOnlySpan combinedDivisors; + private int index; + + public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + this.sqrtLimit = sqrtLimit; + _ = groupBBuffer; + combinedDivisors = GetCombinedDivisors(nMod10); + index = 0; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (index < combinedDivisors.Length) + { + uint entry = combinedDivisors[index++]; + if (entry > sqrtLimit) + { + index = combinedDivisors.Length; + break; + } + + HeuristicDivisorGroup group = ResolveGroup(entry); + ulong value = entry; + byte ending = (byte)(value % 10UL); + ushort residue = (ushort)(value % Wheel210); + candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); + return true; + } + + candidate = default; + return false; + } + } + + internal ref struct MersenneHeuristicDivisorEnumerator + { + private readonly GpuUInt128 step; + private readonly GpuUInt128 limit; + private GpuUInt128 current; + private MersenneDivisorResidueStepper residueStepper; + private bool active; + private ulong processedCount; + private ulong lastDivisor; + + public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + var stepLocal = new GpuUInt128(exponent); + stepLocal.ShiftLeft(1); + step = stepLocal; + + var limitLocal = new GpuUInt128(maxDivisor); + limit = limitLocal; + + var firstDivisor = stepLocal; + firstDivisor.Add(1UL); + + bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0; + + current = hasCandidates ? firstDivisor : GpuUInt128.Zero; + residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default; + active = hasCandidates; + processedCount = 0UL; + lastDivisor = 0UL; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (active) + { + ulong value = current.Low; + processedCount++; + lastDivisor = value; + + bool admissible = residueStepper.IsAdmissible(); + + Advance(); + + if (admissible) + { + candidate = CreateCandidate(value); + return true; + } + } + + candidate = default; + return false; + } + + private void Advance() + { + GpuUInt128 next = current; + next.Add(step); + + if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0) + { + active = false; + return; + } + + current = next; + residueStepper.Advance(); + } + + public ulong ProcessedCount => processedCount; + + public ulong LastDivisor => lastDivisor; + + public bool Exhausted => !active; + + private static HeuristicDivisorCandidate CreateCandidate(ulong value) + { + byte ending = (byte)(value % 10UL); + HeuristicDivisorGroup group; + + if (value == 3UL || value == 7UL || value == 11UL || value == 13UL) + { + group = HeuristicDivisorGroup.GroupAConstant; + } + else if (ending == 3) + { + group = HeuristicDivisorGroup.GroupAWheel; + } + else + { + group = HeuristicDivisorGroup.GroupB; + } + + ushort residue = (ushort)(value % Wheel210); + return new HeuristicDivisorCandidate(value, group, ending, 0, residue); + } + } } diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index c3983a1..ce7e227 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -122,10 +122,10 @@ public sealed class HeuristicPrimeTester return EvaluateWithOpenNumericFallback(n); } - if (!UseHeuristicGroupBTrialDivision) - { - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); - } + // if (!UseHeuristicGroupBTrialDivision) + // { + // return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); + // } bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); return !compositeDetected; @@ -260,8 +260,8 @@ public sealed class HeuristicPrimeTester private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) { - int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); - int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); + int batchCapacity = HeuristicGpuDivisorBatchSize; + int interleaveBatchSize = HeuristicDivisorInterleaveBatchSize; var divisorPool = ThreadStaticPools.UlongPool; var hitPool = ThreadStaticPools.BytePool; @@ -277,7 +277,7 @@ public sealed class HeuristicPrimeTester bool compositeDetected = false; int count = 0; - lock (gpu.ExecutionLock) + // lock (gpu.ExecutionLock) { var input = scratch.Input; var output = scratch.Output; diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 476eea9..7c4e3af 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -1,11 +1,8 @@ -using System; using System.Buffers; using System.Collections.Concurrent; using System.Numerics; using System.Runtime.CompilerServices; -using System.Threading; using ILGPU; -using ILGPU.Algorithms; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; @@ -13,516 +10,500 @@ namespace PerfectNumbers.Core; public sealed class PrimeTester { - public PrimeTester() - { - } - - [ThreadStatic] - private static PrimeTester? _tester; - - public static PrimeTester Exclusive => _tester ??= new(); - - public bool IsPrime(ulong n, CancellationToken ct) - { - if (n <= 1UL) - { - return false; - } - - if (n == 2UL) - { - throw new InvalidOperationException("PrimeTester.IsPrime encountered the sentinel input 2."); - } - - bool isOdd = (n & 1UL) != 0UL; - bool result = isOdd; - - bool requiresTrialDivision = result && n >= 7UL; - - if (requiresTrialDivision) - { - // EvenPerfectBitScanner streams exponents starting at 136,279,841, so the Mod10/GCD guard never fires on the - // production path. Leave the logic commented out as instrumentation for diagnostic builds. - // bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); - // result &= !sharesMaxExponentFactor; - - if (result) - { - uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; - ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; - - ulong nMod10 = n.Mod10(); - switch (nMod10) - { - case 1UL: - smallPrimeDivisors = PrimesGenerator.SmallPrimesLastOne; - smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastOne; - break; - case 3UL: - smallPrimeDivisors = DivisorGenerator.SmallPrimesLastThree; - smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastThree; - break; - case 7UL: - smallPrimeDivisors = PrimesGenerator.SmallPrimesLastSeven; - smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastSeven; - break; - case 9UL: - smallPrimeDivisors = DivisorGenerator.SmallPrimesLastNine; - smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastNine; - break; - } - - int smallPrimeDivisorsLength = smallPrimeDivisors.Length; - for (int i = 0; i < smallPrimeDivisorsLength; i++) - { - if (smallPrimeDivisorsMul[i] > n) - { - break; - } - - if (n % smallPrimeDivisors[i] == 0) - { - result = false; - break; - } - } - } - } - - return result; - } - - public static bool IsPrimeGpu(ulong n) - { - return Exclusive.IsPrimeGpu(n, CancellationToken.None); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool IsPrimeGpu(ulong n, CancellationToken ct) - { - bool forceCpu = GpuContextPool.ForceCpu; - bool belowGpuRange = n < 31UL; - - if (forceCpu || belowGpuRange) - { - return IsPrimeCpu(n, ct); - } - - Span one = stackalloc ulong[1]; - Span outFlags = stackalloc byte[1]; - one[0] = n; - outFlags[0] = 0; - - IsPrimeBatchGpu(one, outFlags); - - return outFlags[0] != 0; - } - - public static bool IsPrimeCpu(ulong n, CancellationToken ct) - { - // Preserve the legacy entry point for callers that bypass the thread-local caches. - return Exclusive.IsPrime(n, ct); - } - - - public static int GpuBatchSize { get; set; } = 262_144; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static void IsPrimeBatchGpu(ReadOnlySpan values, Span results) - { - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = PrimeTesterGpuContextPool.Rent(); - var state = gpu.State; - var accelerator = gpu.Accelerator; - int totalLength = values.Length; - int batchSize = Math.Max(1, GpuBatchSize); - - lock (gpu.ExecutionLock) - { - var scratch = state.RentScratch(batchSize, accelerator); - var input = scratch.Input; - var output = scratch.Output; - ArrayPool pool = ThreadStaticPools.UlongPool; - ulong[] temp = pool.Rent(batchSize); - - int pos = 0; - while (pos < totalLength) - { - int remaining = totalLength - pos; - int count = remaining > batchSize ? batchSize : remaining; - - values.Slice(pos, count).CopyTo(temp); - input.View.CopyFromCPU(ref temp[0], count); - - state.Kernel( - count, - input.View, - state.DevicePrimesDefault.View, - state.DevicePrimesLastOne.View, - state.DevicePrimesLastSeven.View, - state.DevicePrimesLastThree.View, - state.DevicePrimesLastNine.View, - state.DevicePrimesPow2Default.View, - state.DevicePrimesPow2LastOne.View, - state.DevicePrimesPow2LastSeven.View, - state.DevicePrimesPow2LastThree.View, - state.DevicePrimesPow2LastNine.View, - output.View); - accelerator.Synchronize(); - output.View.CopyToCPU(ref results[pos], count); - - pos += count; - } - - pool.Return(temp, clearArray: false); - state.ReturnScratch(scratch); - } - - gpu.Dispose(); - limiter.Dispose(); - } - - internal static class PrimeTesterGpuContextPool - { - internal sealed class PooledContext - { - public Context Context { get; } - - public Accelerator Accelerator { get; } - - public object ExecutionLock { get; } = new(); - - private KernelState? _kernelState; - - public PooledContext() - { - Context = Context.CreateDefault(); - Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context); - _kernelState = GpuKernelState.GetOrCreate(Accelerator); - } - - public KernelState KernelState - { - get - { - var state = _kernelState; - if (state is null || state.IsDisposed) - { - state = GpuKernelState.GetOrCreate(Accelerator); - _kernelState = state; - } - - return state; - } - } - - public void Dispose() - { - ClearGpuCaches(Accelerator); - _kernelState = null; - Accelerator.Dispose(); - Context.Dispose(); - } - } - - private static readonly ConcurrentQueue Pool = new(); - private static readonly object CreationLock = new(); - - internal static PrimeTesterGpuContextLease Rent() - { - if (Pool.TryDequeue(out var ctx)) - { - return new PrimeTesterGpuContextLease(ctx); - } - - lock (CreationLock) - { - return new PrimeTesterGpuContextLease(new PooledContext()); - } - } - - internal static void DisposeAll() - { - while (Pool.TryDequeue(out var ctx)) - { - ctx.Dispose(); - } - } - - private static void Return(PooledContext ctx) - { - ctx.Accelerator.Synchronize(); - Pool.Enqueue(ctx); - } - - internal struct PrimeTesterGpuContextLease - { - private PooledContext? _ctx; - - internal PrimeTesterGpuContextLease(PooledContext ctx) - { - _ctx = ctx; - } - - public Accelerator Accelerator => _ctx!.Accelerator; - - public object ExecutionLock => _ctx!.ExecutionLock; - - public KernelState State => _ctx!.KernelState; - - public void Dispose() - { - Return(_ctx ?? throw new InvalidOperationException("GPU context lease is not initialized.")); - } - } - } - - // Per-accelerator GPU state for prime sieve (kernel + uploaded primes). - internal sealed class KernelState - { - public Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> Kernel { get; } - public Action, ulong, ArrayView> HeuristicTrialDivisionKernel { get; } - public MemoryBuffer1D DevicePrimesDefault { get; } - public MemoryBuffer1D DevicePrimesLastOne { get; } - public MemoryBuffer1D DevicePrimesLastSeven { get; } - public MemoryBuffer1D DevicePrimesLastThree { get; } - public MemoryBuffer1D DevicePrimesLastNine { get; } - public MemoryBuffer1D DevicePrimesPow2Default { get; } - public MemoryBuffer1D DevicePrimesPow2LastOne { get; } - public MemoryBuffer1D DevicePrimesPow2LastSeven { get; } - public MemoryBuffer1D DevicePrimesPow2LastThree { get; } - public MemoryBuffer1D DevicePrimesPow2LastNine { get; } - public bool IsDisposed { get; private set; } - private readonly System.Collections.Concurrent.ConcurrentBag _scratchPool = []; - // TODO: Replace this ConcurrentBag with the lock-free ring buffer variant validated in - // GpuModularArithmeticBenchmarks so renting scratch buffers stops contending on the bag's internal locks when - // thousands of GPU batches execute per second. - - public KernelState(Accelerator accelerator) - { - IsDisposed = false; - // Compile once per accelerator and upload primes once. - Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); - HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel, ulong, ArrayView>(PrimeTesterKernels.HeuristicTrialDivisionKernel); - - var primesDefault = DivisorGenerator.SmallPrimes; - DevicePrimesDefault = accelerator.Allocate1D(primesDefault.Length); - DevicePrimesDefault.View.CopyFromCPU(primesDefault); - - var primesLastOne = DivisorGenerator.SmallPrimesLastOne; - DevicePrimesLastOne = accelerator.Allocate1D(primesLastOne.Length); - DevicePrimesLastOne.View.CopyFromCPU(primesLastOne); - - var primesLastSeven = DivisorGenerator.SmallPrimesLastSeven; - DevicePrimesLastSeven = accelerator.Allocate1D(primesLastSeven.Length); - DevicePrimesLastSeven.View.CopyFromCPU(primesLastSeven); - - var primesLastThree = DivisorGenerator.SmallPrimesLastThree; - DevicePrimesLastThree = accelerator.Allocate1D(primesLastThree.Length); - DevicePrimesLastThree.View.CopyFromCPU(primesLastThree); - - var primesLastNine = DivisorGenerator.SmallPrimesLastNine; - DevicePrimesLastNine = accelerator.Allocate1D(primesLastNine.Length); - DevicePrimesLastNine.View.CopyFromCPU(primesLastNine); - - var primesPow2Default = DivisorGenerator.SmallPrimesPow2; - DevicePrimesPow2Default = accelerator.Allocate1D(primesPow2Default.Length); - DevicePrimesPow2Default.View.CopyFromCPU(primesPow2Default); - - var primesPow2LastOne = DivisorGenerator.SmallPrimesPow2LastOne; - DevicePrimesPow2LastOne = accelerator.Allocate1D(primesPow2LastOne.Length); - DevicePrimesPow2LastOne.View.CopyFromCPU(primesPow2LastOne); - - var primesPow2LastSeven = DivisorGenerator.SmallPrimesPow2LastSeven; - DevicePrimesPow2LastSeven = accelerator.Allocate1D(primesPow2LastSeven.Length); - DevicePrimesPow2LastSeven.View.CopyFromCPU(primesPow2LastSeven); - - var primesPow2LastThree = DivisorGenerator.SmallPrimesPow2LastThree; - DevicePrimesPow2LastThree = accelerator.Allocate1D(primesPow2LastThree.Length); - DevicePrimesPow2LastThree.View.CopyFromCPU(primesPow2LastThree); - - var primesPow2LastNine = DivisorGenerator.SmallPrimesPow2LastNine; - DevicePrimesPow2LastNine = accelerator.Allocate1D(primesPow2LastNine.Length); - DevicePrimesPow2LastNine.View.CopyFromCPU(primesPow2LastNine); - } - - internal sealed class ScratchBuffers - { - public MemoryBuffer1D Input { get; private set; } - public MemoryBuffer1D Output { get; private set; } - public int Capacity { get; private set; } - - public ScratchBuffers(Accelerator accel, int capacity) - { - Capacity = Math.Max(1, capacity); - Input = accel.Allocate1D(Capacity); - Output = accel.Allocate1D(Capacity); - } - - public void Dispose() - { - Output.Dispose(); - Input.Dispose(); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ScratchBuffers RentScratch(int minCapacity, Accelerator accel) - { - minCapacity = Math.Max(1, minCapacity); - while (_scratchPool.TryTake(out var sb)) - { - if (sb.Capacity >= minCapacity) - { - return sb; - } - - sb.Dispose(); - } - - return new ScratchBuffers(accel, minCapacity); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void ReturnScratch(ScratchBuffers scratch) - { - _scratchPool.Add(scratch); - } - - public void Clear() - { - if (IsDisposed) - { - return; - } - - while (_scratchPool.TryTake(out var sb)) - { - sb.Dispose(); - } - - DevicePrimesDefault.Dispose(); - DevicePrimesLastOne.Dispose(); - DevicePrimesLastSeven.Dispose(); - DevicePrimesLastThree.Dispose(); - DevicePrimesLastNine.Dispose(); - DevicePrimesPow2Default.Dispose(); - DevicePrimesPow2LastOne.Dispose(); - DevicePrimesPow2LastSeven.Dispose(); - DevicePrimesPow2LastThree.Dispose(); - DevicePrimesPow2LastNine.Dispose(); - IsDisposed = true; - } - } - - internal static class GpuKernelState - { - // Map accelerator to cached state; use Lazy to serialize kernel creation - private static readonly System.Collections.Concurrent.ConcurrentDictionary> States = new(); - // TODO: Prewarm this per-accelerator cache during startup (and reuse a simple array keyed by accelerator index) - // once the kernel pool exposes deterministic ordering; the Lazy wrappers showed measurable overhead in the - // GpuModularArithmeticBenchmarks hot path. - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static KernelState GetOrCreate(Accelerator accelerator) - { - var lazy = States.GetOrAdd( - accelerator, - acc => new Lazy(() => new KernelState(acc), System.Threading.LazyThreadSafetyMode.ExecutionAndPublication)); - return lazy.Value; - } - - public static void Clear(Accelerator accelerator) - { - if (States.TryRemove(accelerator, out var lazy)) - { - if (lazy.IsValueCreated) - { - var state = lazy.Value; - state.Clear(); - } - } - } - } - - // Expose cache clearing for accelerator disposal coordination - public static void ClearGpuCaches(Accelerator accelerator) - { - GpuKernelState.Clear(accelerator); - } - - internal static void DisposeGpuContexts() - { - PrimeTesterGpuContextPool.DisposeAll(); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool SharesFactorWithMaxExponent(ulong n) - { - // TODO: Replace this on-the-fly GCD probe with the cached factor table derived from - // ResidueComputationBenchmarks so divisor-cycle metadata can short-circuit the test - // instead of recomputing binary GCD for every candidate. - ulong m = (ulong)BitOperations.Log2(n); - return BinaryGcd(n, m) != 1UL; - } - - internal static void SharesFactorWithMaxExponentBatch(ReadOnlySpan values, Span results) - { - // TODO: Route this batch helper through the shared GPU kernel pool from - // GpuUInt128BinaryGcdBenchmarks so we reuse cached kernels, pinned host buffers, - // and divisor-cycle staging instead of allocating new device buffers per call. - var gpu = PrimeTesterGpuContextPool.Rent(); - var accelerator = gpu.Accelerator; - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel); - - int length = values.Length; - var inputBuffer = accelerator.Allocate1D(length); - var resultBuffer = accelerator.Allocate1D(length); - - ArrayPool pool = ThreadStaticPools.UlongPool; - ulong[] temp = pool.Rent(length); - values.CopyTo(temp); - inputBuffer.View.CopyFromCPU(ref temp[0], length); - kernel(length, inputBuffer.View, resultBuffer.View); - accelerator.Synchronize(); - resultBuffer.View.CopyToCPU(ref results[0], length); - pool.Return(temp, clearArray: false); - resultBuffer.Dispose(); - inputBuffer.Dispose(); - gpu.Dispose(); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong BinaryGcd(ulong u, ulong v) - { - // TODO: Swap this handwritten binary GCD for the optimized helper measured in - // GpuUInt128BinaryGcdBenchmarks so CPU callers share the faster subtract-less - // ladder once the common implementation is promoted into PerfectNumbers.Core. - if (u == 0UL) - { - return v; - } - - if (v == 0UL) - { - return u; - } - - int shift = BitOperations.TrailingZeroCount(u | v); - u >>= BitOperations.TrailingZeroCount(u); - - do - { - v >>= BitOperations.TrailingZeroCount(v); - if (u > v) - { - (u, v) = (v, u); - } - - v -= u; - } - while (v != 0UL); - - return u << shift; - } - } + public PrimeTester() + { + } + + [ThreadStatic] + private static PrimeTester? _tester; + + public static PrimeTester Exclusive => _tester ??= new(); + + public static bool IsPrime(ulong n, CancellationToken ct) + { + if (n <= 1UL) + { + return false; + } + + if (n == 2UL) + { + throw new InvalidOperationException("PrimeTester.IsPrime encountered the sentinel input 2."); + } + + bool isOdd = (n & 1UL) != 0UL; + bool result = isOdd; + + bool requiresTrialDivision = result && n >= 7UL; + + if (requiresTrialDivision) + { + // EvenPerfectBitScanner streams exponents starting at 136,279,841, so the Mod10/GCD guard never fires on the + // production path. Leave the logic commented out as instrumentation for diagnostic builds. + // bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); + // result &= !sharesMaxExponentFactor; + + if (result) + { + uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; + ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; + + ulong nMod10 = n.Mod10(); + switch (nMod10) + { + case 1UL: + smallPrimeDivisors = PrimesGenerator.SmallPrimesLastOne; + smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastOne; + break; + case 3UL: + smallPrimeDivisors = DivisorGenerator.SmallPrimesLastThree; + smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastThree; + break; + case 7UL: + smallPrimeDivisors = PrimesGenerator.SmallPrimesLastSeven; + smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastSeven; + break; + case 9UL: + smallPrimeDivisors = DivisorGenerator.SmallPrimesLastNine; + smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastNine; + break; + } + + int smallPrimeDivisorsLength = smallPrimeDivisors.Length; + for (int i = 0; i < smallPrimeDivisorsLength; i++) + { + if (smallPrimeDivisorsMul[i] > n) + { + break; + } + + if (n % smallPrimeDivisors[i] == 0) + { + result = false; + break; + } + } + } + } + + return result; + } + + public static bool IsPrimeGpu(ulong n) + { + return Exclusive.IsPrimeGpu(n, CancellationToken.None); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool IsPrimeGpu(ulong n, CancellationToken ct) + { + bool forceCpu = GpuContextPool.ForceCpu; + bool belowGpuRange = n < 31UL; + + if (forceCpu || belowGpuRange) + { + return IsPrimeCpu(n, ct); + } + + Span one = stackalloc ulong[1]; + Span outFlags = stackalloc byte[1]; + one[0] = n; + outFlags[0] = 0; + + IsPrimeBatchGpu(one, outFlags); + + return outFlags[0] != 0; + } + + public static bool IsPrimeCpu(ulong n, CancellationToken ct) + { + // Preserve the legacy entry point for callers that bypass the thread-local caches. + return IsPrime(n, ct); + } + + + public static int GpuBatchSize { get; set; } = 262_144; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void IsPrimeBatchGpu(ReadOnlySpan values, Span results) + { + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTesterGpuContextPool.Rent(); + var state = gpu.State; + var accelerator = gpu.Accelerator; + int totalLength = values.Length; + int batchSize = Math.Max(1, GpuBatchSize); + + var scratch = state.RentScratch(batchSize, accelerator); + var input = scratch.Input; + var output = scratch.Output; + ArrayPool pool = ThreadStaticPools.UlongPool; + ulong[] temp = pool.Rent(batchSize); + + int pos = 0; + while (pos < totalLength) + { + int remaining = totalLength - pos; + int count = remaining > batchSize ? batchSize : remaining; + + values.Slice(pos, count).CopyTo(temp); + input.View.CopyFromCPU(ref temp[0], count); + + state.Kernel( + count, + input.View, + state.DevicePrimesDefault.View, + state.DevicePrimesLastOne.View, + state.DevicePrimesLastSeven.View, + state.DevicePrimesLastThree.View, + state.DevicePrimesLastNine.View, + state.DevicePrimesPow2Default.View, + state.DevicePrimesPow2LastOne.View, + state.DevicePrimesPow2LastSeven.View, + state.DevicePrimesPow2LastThree.View, + state.DevicePrimesPow2LastNine.View, + output.View); + accelerator.Synchronize(); + output.View.CopyToCPU(ref results[pos], count); + + pos += count; + } + + pool.Return(temp, clearArray: false); + state.ReturnScratch(scratch); + + gpu.Dispose(); + limiter.Dispose(); + } + + internal static class PrimeTesterGpuContextPool + { + internal sealed class PooledContext + { + public Context Context { get; } + + public Accelerator Accelerator { get; } + + private KernelState? _kernelState; + + public PooledContext() + { + Context = Context.CreateDefault(); + Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context); + _kernelState = GpuKernelState.GetOrCreate(Accelerator); + } + + public KernelState KernelState + { + get + { + var state = _kernelState; + if (state is null || state.IsDisposed) + { + state = GpuKernelState.GetOrCreate(Accelerator); + _kernelState = state; + } + + return state; + } + } + + public void Dispose() + { + ClearGpuCaches(Accelerator); + _kernelState = null; + Accelerator.Dispose(); + Context.Dispose(); + } + } + + private static readonly ConcurrentQueue Pool = new(); + + internal static PrimeTesterGpuContextLease Rent() + { + if (Pool.TryDequeue(out var ctx)) + { + return new PrimeTesterGpuContextLease(ctx); + } + + return new PrimeTesterGpuContextLease(new PooledContext()); + } + + internal static void DisposeAll() + { + while (Pool.TryDequeue(out var ctx)) + { + ctx.Dispose(); + } + } + + private static void Return(PooledContext ctx) + { + ctx.Accelerator.Synchronize(); + Pool.Enqueue(ctx); + } + + internal readonly struct PrimeTesterGpuContextLease(PrimeTesterGpuContextPool.PooledContext ctx) + { + private readonly PooledContext? _ctx = ctx; + + public readonly Accelerator Accelerator => _ctx!.Accelerator; + + + public readonly KernelState State => _ctx!.KernelState; + + public readonly void Dispose() + { + Return(_ctx ?? throw new InvalidOperationException("GPU context lease is not initialized.")); + } + } + } + + // Per-accelerator GPU state for prime sieve (kernel + uploaded primes). + internal sealed class KernelState + { + public Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> Kernel; + public Action, ulong, ArrayView> HeuristicTrialDivisionKernel; + public MemoryBuffer1D DevicePrimesDefault; + public MemoryBuffer1D DevicePrimesLastOne; + public MemoryBuffer1D DevicePrimesLastSeven; + public MemoryBuffer1D DevicePrimesLastThree; + public MemoryBuffer1D DevicePrimesLastNine; + public MemoryBuffer1D DevicePrimesPow2Default; + public MemoryBuffer1D DevicePrimesPow2LastOne; + public MemoryBuffer1D DevicePrimesPow2LastSeven; + public MemoryBuffer1D DevicePrimesPow2LastThree; + public MemoryBuffer1D DevicePrimesPow2LastNine; + public bool IsDisposed; + private readonly ConcurrentBag _scratchPool = []; + // TODO: Replace this ConcurrentBag with the lock-free ring buffer variant validated in + // GpuModularArithmeticBenchmarks so renting scratch buffers stops contending on the bag's internal locks when + // thousands of GPU batches execute per second. + + public KernelState(Accelerator accelerator) + { + IsDisposed = false; + // Compile once per accelerator and upload primes once. + Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); + HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel, ulong, ArrayView>(PrimeTesterKernels.HeuristicTrialDivisionKernel); + + var primesDefault = DivisorGenerator.SmallPrimes; + DevicePrimesDefault = accelerator.Allocate1D(primesDefault.Length); + DevicePrimesDefault.View.CopyFromCPU(primesDefault); + + var primesLastOne = DivisorGenerator.SmallPrimesLastOne; + DevicePrimesLastOne = accelerator.Allocate1D(primesLastOne.Length); + DevicePrimesLastOne.View.CopyFromCPU(primesLastOne); + + var primesLastSeven = DivisorGenerator.SmallPrimesLastSeven; + DevicePrimesLastSeven = accelerator.Allocate1D(primesLastSeven.Length); + DevicePrimesLastSeven.View.CopyFromCPU(primesLastSeven); + + var primesLastThree = DivisorGenerator.SmallPrimesLastThree; + DevicePrimesLastThree = accelerator.Allocate1D(primesLastThree.Length); + DevicePrimesLastThree.View.CopyFromCPU(primesLastThree); + + var primesLastNine = DivisorGenerator.SmallPrimesLastNine; + DevicePrimesLastNine = accelerator.Allocate1D(primesLastNine.Length); + DevicePrimesLastNine.View.CopyFromCPU(primesLastNine); + + var primesPow2Default = DivisorGenerator.SmallPrimesPow2; + DevicePrimesPow2Default = accelerator.Allocate1D(primesPow2Default.Length); + DevicePrimesPow2Default.View.CopyFromCPU(primesPow2Default); + + var primesPow2LastOne = DivisorGenerator.SmallPrimesPow2LastOne; + DevicePrimesPow2LastOne = accelerator.Allocate1D(primesPow2LastOne.Length); + DevicePrimesPow2LastOne.View.CopyFromCPU(primesPow2LastOne); + + var primesPow2LastSeven = DivisorGenerator.SmallPrimesPow2LastSeven; + DevicePrimesPow2LastSeven = accelerator.Allocate1D(primesPow2LastSeven.Length); + DevicePrimesPow2LastSeven.View.CopyFromCPU(primesPow2LastSeven); + + var primesPow2LastThree = DivisorGenerator.SmallPrimesPow2LastThree; + DevicePrimesPow2LastThree = accelerator.Allocate1D(primesPow2LastThree.Length); + DevicePrimesPow2LastThree.View.CopyFromCPU(primesPow2LastThree); + + var primesPow2LastNine = DivisorGenerator.SmallPrimesPow2LastNine; + DevicePrimesPow2LastNine = accelerator.Allocate1D(primesPow2LastNine.Length); + DevicePrimesPow2LastNine.View.CopyFromCPU(primesPow2LastNine); + } + + internal sealed class ScratchBuffers + { + public MemoryBuffer1D Input { get; private set; } + public MemoryBuffer1D Output { get; private set; } + public int Capacity { get; private set; } + + public ScratchBuffers(Accelerator accel, int capacity) + { + Capacity = Math.Max(1, capacity); + Input = accel.Allocate1D(Capacity); + Output = accel.Allocate1D(Capacity); + } + + public void Dispose() + { + Output.Dispose(); + Input.Dispose(); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ScratchBuffers RentScratch(int minCapacity, Accelerator accel) + { + while (_scratchPool.TryTake(out var sb)) + { + if (sb.Capacity >= minCapacity) + { + return sb; + } + + sb.Dispose(); + } + + return new ScratchBuffers(accel, minCapacity); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void ReturnScratch(ScratchBuffers scratch) + { + _scratchPool.Add(scratch); + } + + public void Clear() + { + if (IsDisposed) + { + return; + } + + while (_scratchPool.TryTake(out var sb)) + { + sb.Dispose(); + } + + DevicePrimesDefault.Dispose(); + DevicePrimesLastOne.Dispose(); + DevicePrimesLastSeven.Dispose(); + DevicePrimesLastThree.Dispose(); + DevicePrimesLastNine.Dispose(); + DevicePrimesPow2Default.Dispose(); + DevicePrimesPow2LastOne.Dispose(); + DevicePrimesPow2LastSeven.Dispose(); + DevicePrimesPow2LastThree.Dispose(); + DevicePrimesPow2LastNine.Dispose(); + IsDisposed = true; + } + } + + internal static class GpuKernelState + { + // Map accelerator to cached state; use Lazy to serialize kernel creation + private static readonly System.Collections.Concurrent.ConcurrentDictionary> States = new(); + // TODO: Prewarm this per-accelerator cache during startup (and reuse a simple array keyed by accelerator index) + // once the kernel pool exposes deterministic ordering; the Lazy wrappers showed measurable overhead in the + // GpuModularArithmeticBenchmarks hot path. + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static KernelState GetOrCreate(Accelerator accelerator) + { + var lazy = States.GetOrAdd( + accelerator, + acc => new Lazy(() => new KernelState(acc), LazyThreadSafetyMode.ExecutionAndPublication)); + return lazy.Value; + } + + public static void Clear(Accelerator accelerator) + { + if (States.TryRemove(accelerator, out var lazy)) + { + if (lazy.IsValueCreated) + { + var state = lazy.Value; + state.Clear(); + } + } + } + } + + // Expose cache clearing for accelerator disposal coordination + public static void ClearGpuCaches(Accelerator accelerator) + { + GpuKernelState.Clear(accelerator); + } + + internal static void DisposeGpuContexts() + { + PrimeTesterGpuContextPool.DisposeAll(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool SharesFactorWithMaxExponent(ulong n) + { + // TODO: Replace this on-the-fly GCD probe with the cached factor table derived from + // ResidueComputationBenchmarks so divisor-cycle metadata can short-circuit the test + // instead of recomputing binary GCD for every candidate. + ulong m = (ulong)BitOperations.Log2(n); + return BinaryGcd(n, m) != 1UL; + } + + internal static void SharesFactorWithMaxExponentBatch(ReadOnlySpan values, Span results) + { + // TODO: Route this batch helper through the shared GPU kernel pool from + // GpuUInt128BinaryGcdBenchmarks so we reuse cached kernels, pinned host buffers, + // and divisor-cycle staging instead of allocating new device buffers per call. + var gpu = PrimeTesterGpuContextPool.Rent(); + var accelerator = gpu.Accelerator; + var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel); + + int length = values.Length; + var inputBuffer = accelerator.Allocate1D(length); + var resultBuffer = accelerator.Allocate1D(length); + + ArrayPool pool = ThreadStaticPools.UlongPool; + ulong[] temp = pool.Rent(length); + values.CopyTo(temp); + inputBuffer.View.CopyFromCPU(ref temp[0], length); + kernel(length, inputBuffer.View, resultBuffer.View); + accelerator.Synchronize(); + resultBuffer.View.CopyToCPU(ref results[0], length); + pool.Return(temp, clearArray: false); + resultBuffer.Dispose(); + inputBuffer.Dispose(); + gpu.Dispose(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong BinaryGcd(ulong u, ulong v) + { + // TODO: Swap this handwritten binary GCD for the optimized helper measured in + // GpuUInt128BinaryGcdBenchmarks so CPU callers share the faster subtract-less + // ladder once the common implementation is promoted into PerfectNumbers.Core. + if (u == 0UL) + { + return v; + } + + if (v == 0UL) + { + return u; + } + + int shift = BitOperations.TrailingZeroCount(u | v); + u >>= BitOperations.TrailingZeroCount(u); + + do + { + v >>= BitOperations.TrailingZeroCount(v); + if (u > v) + { + (u, v) = (v, u); + } + + v -= u; + } + while (v != 0UL); + + return u << shift; + } +}