diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index f63e315..db1c2fb 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -29,8 +29,7 @@ public class PrimeTesterBenchmarks return; } - CancellationToken cancellationToken = CancellationToken.None; - _ = _heuristicTester.IsPrimeCpu(_candidates[0], cancellationToken); + _ = _heuristicTester.IsPrimeCpu(_candidates[0]); _ = _heuristicTester.IsPrimeGpu(_candidates[0]); _ = Prime.Numbers.IsPrime(_candidates[0]); } @@ -40,12 +39,11 @@ public class PrimeTesterBenchmarks { HeuristicPrimeTester tester = _heuristicTester; ulong[] values = _candidates; - CancellationToken cancellationToken = CancellationToken.None; int primeCount = 0; for (int i = 0; i < values.Length; i++) { - if (tester.IsPrimeCpu(values[i], cancellationToken)) + if (tester.IsPrimeCpu(values[i])) { primeCount++; } @@ -75,7 +73,6 @@ public class PrimeTesterBenchmarks [Benchmark] public int HeuristicCombinedGpu() { - HeuristicCombinedPrimeTester tester = _heuristicCombinedTester; ulong[] values = _candidates; int primeCount = 0; @@ -93,7 +90,6 @@ public class PrimeTesterBenchmarks [Benchmark] public int HeuristicCombinedCpu() { - HeuristicCombinedPrimeTester tester = _heuristicCombinedTester; ulong[] values = _candidates; int primeCount = 0; @@ -115,12 +111,11 @@ public class PrimeTesterBenchmarks { PrimeTester tester = _nonHeuristicTester; ulong[] values = _candidates; - CancellationToken cancellationToken = CancellationToken.None; int primeCount = 0; for (int i = 0; i < values.Length; i++) { - if (tester.IsPrimeGpu(values[i], cancellationToken)) + if (tester.IsPrimeGpu(values[i])) { primeCount++; } @@ -133,12 +128,11 @@ public class PrimeTesterBenchmarks public int NonHeuristicCpu() { ulong[] values = _candidates; - CancellationToken cancellationToken = CancellationToken.None; int primeCount = 0; for (int i = 0; i < values.Length; i++) { - if (PrimeTester.IsPrimeCpu(values[i], cancellationToken)) + if (PrimeTester.IsPrime(values[i])) { primeCount++; } diff --git a/EvenPerfectBitScanner.Tests/ProgramTests.cs b/EvenPerfectBitScanner.Tests/ProgramTests.cs index f08a676..cf16bdb 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) { - PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected); + PrimeTester.IsPrime(n).Should().Be(expected); } [Fact] diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 86b4c20..bf755d4 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -622,8 +622,8 @@ internal static class Program if (!_cliArguments.UseByDivisor) { if (!(GpuContextPool.ForceCpu - ? PrimeTester.IsPrime(p, CancellationToken.None) - : PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None))) + ? PrimeTester.IsPrime(p) + : PrimeTesters.Value!.IsPrimeGpu(p))) { _lastCompositeP = true; return false; diff --git a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs index 5e8faf9..6213b07 100644 --- a/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs +++ b/PerfectNumbers.Core.Tests/Gpu/PrimeTesterGpuTests.cs @@ -17,7 +17,7 @@ public class PrimeTesterGpuTests ulong[] primes = [31UL, 61UL, 89UL, 107UL, 127UL, 521UL]; foreach (ulong prime in primes) { - tester.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue(); + tester.IsPrimeGpu(prime).Should().BeTrue(); } } finally @@ -37,7 +37,7 @@ public class PrimeTesterGpuTests ulong[] composites = [33UL, 39UL, 51UL, 77UL, 91UL]; foreach (ulong composite in composites) { - tester.IsPrimeGpu(composite, CancellationToken.None).Should().BeFalse(); + tester.IsPrimeGpu(composite).Should().BeFalse(); } } finally diff --git a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs index 2d76983..d3d2051 100644 --- a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs +++ b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs @@ -19,7 +19,7 @@ public class HeuristicCombinedPrimeTesterTests Span groupBState = stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4]; - var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500, lastDigit, groupBState); + var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500UL * 500UL, lastDigit, groupBState); var results = new List(32); while (results.Count < 20 && enumerator.TryGetNext(out var candidate)) @@ -230,14 +230,10 @@ public class HeuristicCombinedPrimeTesterTests [InlineData(18_446_744_073_709_551_615UL)] [InlineData(18_446_744_073_709_551_614UL)] [InlineData(18_446_744_073_709_551_613UL)] - public void ComputeHeuristicSqrt_returns_floor_square_root(ulong value) + public void ComputeHeuristicDivisorSquareLimit_returns_input_value(ulong value) { - ulong sqrt = HeuristicCombinedPrimeTester.ComputeHeuristicSqrt(value); + ulong limit = HeuristicCombinedPrimeTester.ComputeHeuristicDivisorSquareLimit(value); - bool squareBelowOrEqual = (UInt128)sqrt * sqrt <= (UInt128)value; - bool nextSquareAbove = (UInt128)(sqrt + 1UL) * (UInt128)(sqrt + 1UL) > (UInt128)value; - - squareBelowOrEqual.Should().BeTrue(); - nextSquareAbove.Should().BeTrue(); + limit.Should().Be(value); } } diff --git a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs index e86e1db..f7dd5b7 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(); - PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected); + PrimeTester.IsPrime(n).Should().Be(expected); } [Theory] @@ -40,7 +40,7 @@ public class PrimeTesterIsPrimeTests GpuContextPool.ForceCpu = false; var tester = new PrimeTester(); - tester.IsPrimeGpu(n, CancellationToken.None).Should().Be(expected); + tester.IsPrimeGpu(n).Should().Be(expected); } [Fact] @@ -50,9 +50,9 @@ public class PrimeTesterIsPrimeTests var tester = new PrimeTester(); GpuContextPool.ForceCpu = true; - tester.IsPrimeGpu(137UL, CancellationToken.None).Should().BeTrue(); + tester.IsPrimeGpu(137UL).Should().BeTrue(); - tester.IsPrimeGpu(341UL, CancellationToken.None).Should().BeFalse(); + tester.IsPrimeGpu(341UL).Should().BeFalse(); GpuContextPool.ForceCpu = false; } @@ -71,7 +71,7 @@ public class PrimeTesterIsPrimeTests var testers = new ThreadLocal(() => new PrimeTester(), trackAllValues: true); Parallel.For(0, primes.Length, i => { - results[i] = testers.Value!.IsPrimeGpu(primes[i], CancellationToken.None); + results[i] = testers.Value!.IsPrimeGpu(primes[i]); }); results.Should().AllBeEquivalentTo(true); @@ -80,7 +80,7 @@ public class PrimeTesterIsPrimeTests { foreach (ulong prime in primes) { - testers.Value!.IsPrimeGpu(prime, CancellationToken.None).Should().BeTrue(); + testers.Value!.IsPrimeGpu(prime).Should().BeTrue(); } }); } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index f60596a..5b9bf5c 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -48,7 +48,7 @@ internal static class PrimeTesterKernels break; } - int length = XMath.Min((int)primes.Length, (int)primeSquares.Length); + int length = (int)primes.Length; for (int i = 0; i < length; i++) { ulong prime = primes[i]; @@ -68,16 +68,10 @@ internal static class PrimeTesterKernels results[index] = result; } - public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView divisors, ulong n, ArrayView results) + public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView divisors, ulong n, ArrayView resultFlag) { ulong divisor = divisors[index]; - if (divisor <= 1UL) - { - results[index] = 0; - return; - } - - results[index] = n % divisor == 0UL ? (byte)1 : (byte)0; + Atomic.Or(ref resultFlag[0], divisor <= 1UL ? 0 : (n % divisor == 0UL ? 1 : 0)); } public static void SharesFactorKernel(Index1D index, ArrayView numbers, ArrayView results) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 2c9ca8b..27a066a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -292,6 +292,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv coveredRange = true; return true; } + + divisorSpan[admissibleCount] = divisorValue; + divisorDataSpan[admissibleCount] = new GpuDivisorPartialData(divisorValue); + offsetSpan[admissibleCount] = admissibleCount; + countSpan[admissibleCount] = 1; + cycleSpan[admissibleCount] = divisorCycle; + exponentSpan[admissibleCount] = prime; + admissibleCount++; } nextDivisor128 += twoP128; @@ -444,13 +452,25 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv // { // return 0UL; // } + ulong cappedLimit = divisorLimit; + const ulong TestPrimeDivisorCap = 170_000_000UL; + if (prime < 10_000_000UL && cappedLimit > TestPrimeDivisorCap) + { + cappedLimit = TestPrimeDivisorCap; + } + if (prime - 1UL >= 64UL) { - return divisorLimit; + return cappedLimit; } ulong computedLimit = (1UL << (int)(prime - 1UL)) - 1UL; - return computedLimit < divisorLimit ? computedLimit : divisorLimit; + if (prime < 10_000_000UL && computedLimit > cappedLimit) + { + computedLimit = cappedLimit; + } + + return computedLimit < cappedLimit ? computedLimit : cappedLimit; } public sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 6175a28..72add2f 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -17,29 +17,33 @@ public sealed class HeuristicCombinedPrimeTester { HeuristicPrimeSieves.EnsureInitialized(); - CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding1TwoAOneBSquares); + CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding3TwoAOneBSquares); + CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding7TwoAOneBSquares); + CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB, out CombinedDivisorsEnding9TwoAOneBSquares); - CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding1OneAOneBSquares); + CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding3OneAOneBSquares); + CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding7OneAOneBSquares); + CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB, out CombinedDivisorsEnding9OneAOneBSquares); - CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding1ThreeAOneBSquares); + CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding3ThreeAOneBSquares); + CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding7ThreeAOneBSquares); + CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB, out CombinedDivisorsEnding9ThreeAOneBSquares); - CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding1ThreeATwoBSquares); + CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding3ThreeATwoBSquares); + CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding7ThreeATwoBSquares); + CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB, out CombinedDivisorsEnding9ThreeATwoBSquares); CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB; CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB; CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB; CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB; + CombinedDivisorsEnding1Squares = CombinedDivisorsEnding1OneAOneBSquares; + CombinedDivisorsEnding3Squares = CombinedDivisorsEnding3OneAOneBSquares; + CombinedDivisorsEnding7Squares = CombinedDivisorsEnding7OneAOneBSquares; + CombinedDivisorsEnding9Squares = CombinedDivisorsEnding9OneAOneBSquares; } private const ulong Wheel210 = 210UL; @@ -71,6 +75,31 @@ public sealed class HeuristicCombinedPrimeTester private static readonly uint[] CombinedDivisorsEnding7; private static readonly uint[] CombinedDivisorsEnding9; + private static readonly ulong[] CombinedDivisorsEnding1TwoAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding3TwoAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding7TwoAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding9TwoAOneBSquares; + + private static readonly ulong[] CombinedDivisorsEnding1OneAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding3OneAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding7OneAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding9OneAOneBSquares; + + private static readonly ulong[] CombinedDivisorsEnding1ThreeAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding3ThreeAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding7ThreeAOneBSquares; + private static readonly ulong[] CombinedDivisorsEnding9ThreeAOneBSquares; + + private static readonly ulong[] CombinedDivisorsEnding1ThreeATwoBSquares; + private static readonly ulong[] CombinedDivisorsEnding3ThreeATwoBSquares; + private static readonly ulong[] CombinedDivisorsEnding7ThreeATwoBSquares; + private static readonly ulong[] CombinedDivisorsEnding9ThreeATwoBSquares; + + private static readonly ulong[] CombinedDivisorsEnding1Squares; + private static readonly ulong[] CombinedDivisorsEnding3Squares; + private static readonly ulong[] CombinedDivisorsEnding7Squares; + private static readonly ulong[] CombinedDivisorsEnding9Squares; + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); internal enum CombinedDivisorPattern : byte @@ -123,53 +152,51 @@ public sealed class HeuristicCombinedPrimeTester public static bool IsPrimeCpu(ulong n) { byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); + ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); - if (sqrtLimit < 3UL) + if (maxDivisorSquare < 9UL) { return EvaluateWithOpenNumericFallback(n); } - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); + return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10); } public static bool IsPrimeGpu(ulong n) { byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); - return IsPrimeGpu(n, sqrtLimit, nMod10); + ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); + return IsPrimeGpu(n, maxDivisorSquare, nMod10); } - public static bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) + public static bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10) { // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? - if (sqrtLimit < 3UL) + if (maxDivisorSquare < 9UL) { return EvaluateWithOpenNumericFallback(n); } - bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); + bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10); return !compositeDetected; } - private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10) + private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10) { ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); - if (combinedDivisors.IsEmpty) - { - return EvaluateWithOpenNumericFallback(n); - } + ReadOnlySpan combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); int length = combinedDivisors.Length; + for (int i = 0; i < length; i++) { - uint entry = combinedDivisors[i]; - if (entry > sqrtLimit) + ulong divisorSquare = combinedDivisorSquares[i]; + if (divisorSquare > maxDivisorSquare) { break; } - if (n % entry == 0UL) + if (n % combinedDivisors[i] == 0UL) { return false; } @@ -178,127 +205,87 @@ public sealed class HeuristicCombinedPrimeTester return EvaluateWithOpenNumericFallback(n); } - private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) + private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10) { - int batchCapacity = HeuristicGpuDivisorBatchSize; - var divisorPool = ThreadStaticPools.UlongPool; - var hitPool = ThreadStaticPools.BytePool; + int batchCapacity = HeuristicGpuDivisorBatchSize; + var divisorPool = ThreadStaticPools.UlongPool; - ulong[] divisorArray = divisorPool.Rent(batchCapacity); - byte[] hitFlags = hitPool.Rent(batchCapacity); + ulong[] divisorArray = divisorPool.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); + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity); + var accelerator = gpu.Accelerator; + var kernel = gpu.HeuristicTrialDivisionKernel; - bool compositeDetected = false; - int count = 0; + bool compositeDetected = false; + int count = 0; - var input = scratch.Input; - var output = scratch.Output; + var inputView = gpu.Input.View; + var flagView = gpu.HeuristicFlag.View; - bool ProcessBatch(int length) - { - input.View.CopyFromCPU(ref divisorArray[0], length); - state.HeuristicTrialDivisionKernel(length, input.View, n, output.View); - accelerator.Synchronize(); - output.View.CopyToCPU(ref hitFlags[0], length); + bool ProcessBatch(int length) + { + inputView.CopyFromCPU(ref divisorArray[0], length); + int compositeFlag = 0; + flagView.CopyFromCPU(ref compositeFlag, 1); + kernel(length, inputView, n, flagView); + accelerator.Synchronize(); + flagView.CopyToCPU(ref compositeFlag, 1); - for (int i = 0; i < length; i++) - { - if (hitFlags[i] == 0) - { - continue; - } + return compositeFlag != 0; + } - return true; - // ulong divisor = divisorArray[i]; - // if (divisor > 1UL && n % divisor == 0UL) - // { - // return true; - // } - } + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + ReadOnlySpan combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); - return false; - } + int length = combinedDivisors.Length; - ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); - if (combinedDivisors.IsEmpty) - { - goto Cleanup; - } + for (int i = 0; i < length; i++) + { + ulong divisorSquare = combinedDivisorSquares[i]; + if (divisorSquare > maxDivisorSquare) + { + break; + } - int length = combinedDivisors.Length; - for (int i = 0; i < length; i++) - { - uint entry = combinedDivisors[i]; - if (entry > sqrtLimit) - { - break; - } + divisorArray[count] = combinedDivisors[i]; + count++; - divisorArray[count] = entry; - count++; + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } + count = 0; + } + } - count = 0; - } - } + if (count > 0 && ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } - 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; +Cleanup: + divisorPool.Return(divisorArray); + 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; - } + internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n) + { + return n; + } private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { @@ -308,7 +295,19 @@ public sealed class HeuristicCombinedPrimeTester 3 => CombinedDivisorsEnding3, 7 => CombinedDivisorsEnding7, 9 => CombinedDivisorsEnding9, - _ => ReadOnlySpan.Empty, + _ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."), + }; + } + + private static ReadOnlySpan GetCombinedDivisorSquares(byte nMod10) + { + return nMod10 switch + { + 1 => CombinedDivisorsEnding1Squares, + 3 => CombinedDivisorsEnding3Squares, + 7 => CombinedDivisorsEnding7Squares, + 9 => CombinedDivisorsEnding9Squares, + _ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."), }; } @@ -322,7 +321,7 @@ public sealed class HeuristicCombinedPrimeTester 3 => CombinedDivisorsEnding3TwoAOneB, 7 => CombinedDivisorsEnding7TwoAOneB, 9 => CombinedDivisorsEnding9TwoAOneB, - _ => ReadOnlySpan.Empty, + _ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."), }, CombinedDivisorPattern.OneAOneB => nMod10 switch { @@ -330,7 +329,7 @@ public sealed class HeuristicCombinedPrimeTester 3 => CombinedDivisorsEnding3OneAOneB, 7 => CombinedDivisorsEnding7OneAOneB, 9 => CombinedDivisorsEnding9OneAOneB, - _ => ReadOnlySpan.Empty, + _ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."), }, CombinedDivisorPattern.ThreeAOneB => nMod10 switch { @@ -338,7 +337,7 @@ public sealed class HeuristicCombinedPrimeTester 3 => CombinedDivisorsEnding3ThreeAOneB, 7 => CombinedDivisorsEnding7ThreeAOneB, 9 => CombinedDivisorsEnding9ThreeAOneB, - _ => ReadOnlySpan.Empty, + _ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."), }, CombinedDivisorPattern.ThreeATwoB => nMod10 switch { @@ -346,13 +345,53 @@ public sealed class HeuristicCombinedPrimeTester 3 => CombinedDivisorsEnding3ThreeATwoB, 7 => CombinedDivisorsEnding7ThreeATwoB, 9 => CombinedDivisorsEnding9ThreeATwoB, - _ => ReadOnlySpan.Empty, + _ => throw new InvalidOperationException($"Unsupported combined divisor selector for digit {nMod10}."), }, - _ => ReadOnlySpan.Empty, + _ => throw new InvalidOperationException($"Unsupported combined divisor pattern: {pattern}."), }; } - private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) + internal static ReadOnlySpan GetCombinedDivisorSquares(byte nMod10, CombinedDivisorPattern pattern) + { + return pattern switch + { + CombinedDivisorPattern.TwoAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1TwoAOneBSquares, + 3 => CombinedDivisorsEnding3TwoAOneBSquares, + 7 => CombinedDivisorsEnding7TwoAOneBSquares, + 9 => CombinedDivisorsEnding9TwoAOneBSquares, + _ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."), + }, + CombinedDivisorPattern.OneAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1OneAOneBSquares, + 3 => CombinedDivisorsEnding3OneAOneBSquares, + 7 => CombinedDivisorsEnding7OneAOneBSquares, + 9 => CombinedDivisorsEnding9OneAOneBSquares, + _ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."), + }, + CombinedDivisorPattern.ThreeAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeAOneBSquares, + 3 => CombinedDivisorsEnding3ThreeAOneBSquares, + 7 => CombinedDivisorsEnding7ThreeAOneBSquares, + 9 => CombinedDivisorsEnding9ThreeAOneBSquares, + _ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."), + }, + CombinedDivisorPattern.ThreeATwoB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeATwoBSquares, + 3 => CombinedDivisorsEnding3ThreeATwoBSquares, + 7 => CombinedDivisorsEnding7ThreeATwoBSquares, + 9 => CombinedDivisorsEnding9ThreeATwoBSquares, + _ => throw new InvalidOperationException($"Unsupported combined divisor square selector for digit {nMod10}."), + }, + _ => throw new InvalidOperationException($"Unsupported combined divisor pattern for squares: {pattern}."), + }; + } + + private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern, out ulong[] squares) { ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; ReadOnlySpan groupB = GetGroupBDivisors(nMod10); @@ -412,7 +451,16 @@ public sealed class HeuristicCombinedPrimeTester } } - return combined.ToArray(); + uint[] result = combined.ToArray(); + ulong[] resultSquares = new ulong[result.Length]; + for (int i = 0; i < result.Length; i++) + { + ulong divisor = result[i]; + resultSquares[i] = divisor * divisor; + } + + squares = resultSquares; + return result; } @@ -447,9 +495,9 @@ public sealed class HeuristicCombinedPrimeTester }; [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) { - return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); + return new HeuristicDivisorEnumerator(maxDivisorSquare, nMod10, groupBBuffer); } internal struct HeuristicGroupBSequenceState @@ -530,44 +578,50 @@ public sealed class HeuristicCombinedPrimeTester } internal ref struct HeuristicDivisorEnumerator - { - private readonly ulong sqrtLimit; - private readonly ReadOnlySpan combinedDivisors; - private int index; + { + private readonly ulong maxDivisorSquare; + private readonly ReadOnlySpan combinedDivisors; + private readonly ReadOnlySpan combinedDivisorSquares; + private int index; - public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) - { - this.sqrtLimit = sqrtLimit; - _ = groupBBuffer; - combinedDivisors = GetCombinedDivisors(nMod10); - index = 0; - } + public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) + { + this.maxDivisorSquare = maxDivisorSquare; + _ = groupBBuffer; + combinedDivisors = GetCombinedDivisors(nMod10); + combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); + index = 0; + } - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - while (index < combinedDivisors.Length) - { - uint entry = combinedDivisors[index++]; - if (entry > sqrtLimit) - { - index = combinedDivisors.Length; - break; - } + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (index < combinedDivisors.Length) + { + int currentIndex = index; + ulong divisorSquare = combinedDivisorSquares[currentIndex]; + if (divisorSquare > maxDivisorSquare) + { + 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; - } + uint entry = combinedDivisors[currentIndex]; + index = currentIndex + 1; - candidate = default; - return false; - } - } + 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; + } - internal ref struct MersenneHeuristicDivisorEnumerator + candidate = default; + return false; + } + } + + internal ref struct MersenneHeuristicDivisorEnumerator { private readonly GpuUInt128 step; private readonly GpuUInt128 limit; diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index ce7e227..2f491df 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -21,8 +21,6 @@ public sealed class HeuristicPrimeTester } // TODO: Don't use these tiny arrays. Hard-code the values / checks instead - private static readonly ulong[] GroupAConstantDivisors = [3UL, 7UL, 11UL, 13UL]; - private static readonly byte[] GroupAIncrementPattern = [20, 10]; private static readonly byte[] GroupBEndingOrderMod1 = [9, 1]; private static readonly byte[] GroupBEndingOrderMod3 = [9, 7]; private static readonly byte[] GroupBEndingOrderMod7 = [7, 1]; @@ -30,10 +28,8 @@ public sealed class HeuristicPrimeTester // TODO: Would it speed up things if we implement something like ExponentRemainderStepper, // MersenneDivisorResidueStepper, or CycleRemainderStepper work with these? - private static readonly ushort[] Wheel210ResiduesEnding1 = [1, 11, 31, 41, 61, 71, 101, 121, 131, 151, 181, 191]; - private static readonly ushort[] Wheel210ResiduesEnding7 = [17, 37, 47, 67, 97, 107, 127, 137, 157, 167, 187, 197]; - private static readonly ushort[] Wheel210ResiduesEnding9 = [19, 29, 59, 79, 89, 109, 139, 149, 169, 179, 199, 209]; private const ulong Wheel210 = 210UL; + private const int GroupAConstantCount = 4; internal const int MaxGroupBSequences = 36; private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. @@ -81,76 +77,73 @@ public sealed class HeuristicPrimeTester public static HeuristicPrimeTester Exclusive => _tester ??= new(); - public bool IsPrimeCpu(ulong n, CancellationToken ct) + public bool IsPrimeCpu(ulong n) { - ct.ThrowIfCancellationRequested(); - byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); + ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); - if (sqrtLimit < 3UL) + if (maxDivisorSquare < 9UL) { return EvaluateWithOpenNumericFallback(n); } bool includeGroupB = UseHeuristicGroupBTrialDivision; - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB); - } - - public bool IsPrimeCpu(ulong n) - { - return IsPrimeCpu(n, CancellationToken.None); + return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB); } public bool IsPrimeGpu(ulong n) { byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); - return IsPrimeGpu(n, sqrtLimit, nMod10); + ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); + return IsPrimeGpu(n, maxDivisorSquare, nMod10); } - public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) + public bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10) { - return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10); + return HeuristicIsPrimeGpuCore(n, maxDivisorSquare, nMod10); } - private bool HeuristicIsPrimeGpuCore(ulong n, ulong sqrtLimit, byte nMod10) + private bool HeuristicIsPrimeGpuCore(ulong n, ulong maxDivisorSquare, byte nMod10) { // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? - if (sqrtLimit < 3UL) + if (maxDivisorSquare < 9UL) { return EvaluateWithOpenNumericFallback(n); } // if (!UseHeuristicGroupBTrialDivision) // { - // return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); + // return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB: false); // } - bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); + bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10); return !compositeDetected; } - private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB) + + private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10, bool includeGroupB) { ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; + ReadOnlySpan groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); int groupAIndex = 0; if (!includeGroupB) { - while (groupAIndex < groupADivisors.Length) + int groupALength = groupADivisors.Length; + while (groupAIndex < groupALength) { int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) + while (processed < interleaveBatchSize && groupAIndex < groupALength) { - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (divisor > sqrtLimit) + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) { - groupAIndex = groupADivisors.Length; + groupAIndex = groupALength; break; } + ulong divisor = (ulong)groupADivisors[groupAIndex]; if (n % divisor == 0UL) { return false; @@ -160,14 +153,14 @@ public sealed class HeuristicPrimeTester processed++; } - if (groupAIndex >= groupADivisors.Length) + if (groupAIndex >= groupALength) { break; } - if ((ulong)groupADivisors[groupAIndex] > sqrtLimit) + if (groupADivisorSquares[groupAIndex] > maxDivisorSquare) { - groupAIndex = groupADivisors.Length; + groupAIndex = groupALength; break; } } @@ -185,8 +178,8 @@ public sealed class HeuristicPrimeTester ? stackalloc int[endingOrder.Length] : new int[endingOrder.Length]; - bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; - bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); + bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare; + bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); while (groupAHasMore || groupBHasMore) { @@ -195,14 +188,15 @@ public sealed class HeuristicPrimeTester int processed = 0; while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) { - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (divisor > sqrtLimit) + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) { groupAIndex = groupADivisors.Length; groupAHasMore = false; break; } + ulong divisor = (ulong)groupADivisors[groupAIndex]; if (n % divisor == 0UL) { return false; @@ -218,7 +212,7 @@ public sealed class HeuristicPrimeTester } else { - groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; + groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare; if (!groupAHasMore) { groupAIndex = groupADivisors.Length; @@ -231,7 +225,7 @@ public sealed class HeuristicPrimeTester int processed = 0; while (processed < interleaveBatchSize) { - if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) + if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor)) { groupBHasMore = false; break; @@ -247,7 +241,7 @@ public sealed class HeuristicPrimeTester if (groupBHasMore) { - groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); + groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); } } } @@ -258,198 +252,162 @@ public sealed class HeuristicPrimeTester - private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) + + + + private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10) { int batchCapacity = HeuristicGpuDivisorBatchSize; int interleaveBatchSize = HeuristicDivisorInterleaveBatchSize; 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 gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity); var accelerator = gpu.Accelerator; - var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); - var scratch = state.RentScratch(batchCapacity, accelerator); + var kernel = gpu.HeuristicTrialDivisionKernel; bool compositeDetected = false; int count = 0; - // lock (gpu.ExecutionLock) + var inputView = gpu.Input.View; + var flagView = gpu.HeuristicFlag.View; + + bool ProcessBatch(int length) { - var input = scratch.Input; - var output = scratch.Output; - - bool ProcessBatch(int length) + if (length <= 0) { - 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 groupADivisors = HeuristicPrimeSieves.GroupADivisors; - int groupAIndex = 0; - ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); - Span indices = endingOrder.Length <= 8 - ? stackalloc int[endingOrder.Length] - : new int[endingOrder.Length]; + inputView.CopyFromCPU(ref divisorArray[0], length); + int compositeFlag = 0; + flagView.CopyFromCPU(ref compositeFlag, 1); + kernel(length, inputView, n, flagView); + accelerator.Synchronize(); + flagView.CopyToCPU(ref compositeFlag, 1); + return compositeFlag != 0; + } - bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; - bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, sqrtLimit); + ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; + ReadOnlySpan groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; + int groupAIndex = 0; + ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); + Span indices = endingOrder.Length <= 8 + ? stackalloc int[endingOrder.Length] + : new int[endingOrder.Length]; - while (!compositeDetected && (groupAHasMore || groupBHasMore)) + bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare; + bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); + + while (!compositeDetected && (groupAHasMore || groupBHasMore)) + { + if (groupAHasMore) { - if (groupAHasMore) + int processed = 0; + while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) { - int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) { - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (divisor > sqrtLimit) + groupAIndex = groupADivisors.Length; + groupAHasMore = false; + break; + } + + divisorArray[count] = (ulong)groupADivisors[groupAIndex]; + count++; + + if (count == batchCapacity) + { + if (ProcessBatch(count)) { - groupAIndex = groupADivisors.Length; - groupAHasMore = false; + compositeDetected = true; break; } - divisorArray[count] = divisor; - count++; + count = 0; + } - if (count == batchCapacity) + groupAIndex++; + processed++; + } + + if (compositeDetected) + { + break; + } + + if (groupAIndex >= groupADivisors.Length) + { + groupAHasMore = false; + } + else + { + groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare; + if (!groupAHasMore) + { + groupAIndex = groupADivisors.Length; + } + } + } + + if (!compositeDetected && groupBHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize) + { + if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor)) + { + groupBHasMore = false; + break; + } + + divisorArray[count] = divisor; + count++; + + if (count == batchCapacity) + { + if (ProcessBatch(count)) { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; + compositeDetected = true; + break; } - groupAIndex++; - processed++; + count = 0; } - if (groupAIndex >= groupADivisors.Length) - { - groupAHasMore = false; - } - else - { - groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; - if (!groupAHasMore) - { - groupAIndex = groupADivisors.Length; - } - } + processed++; } if (!compositeDetected && groupBHasMore) { - int processed = 0; - while (processed < interleaveBatchSize) - { - if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) - { - groupBHasMore = false; - break; - } - - divisorArray[count] = divisor; - count++; - - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; - } - - processed++; - } - - if (groupBHasMore) - { - groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); - } + groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); } } - - if (!compositeDetected && count > 0) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; - } } -Cleanup: - state.ReturnScratch(scratch); + if (!compositeDetected && count > 0 && ProcessBatch(count)) + { + compositeDetected = true; + } + 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) + internal static ulong ComputeHeuristicDivisorSquareLimit(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; + return n; } private static ReadOnlySpan GetGroupBEndingOrder(byte nMod10) => nMod10 switch @@ -462,8 +420,9 @@ Cleanup: }; + [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool TrySelectNextGroupBDivisor(ReadOnlySpan endingOrder, Span indices, ulong sqrtLimit, out ulong divisor) + private static bool TrySelectNextGroupBDivisor(ReadOnlySpan endingOrder, Span indices, ulong maxDivisorSquare, out ulong divisor) { ulong bestCandidate = ulong.MaxValue; int bestEndingIndex = -1; @@ -471,19 +430,21 @@ Cleanup: for (int i = 0; i < endingOrder.Length; i++) { ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + ReadOnlySpan squares = GetGroupBDivisorSquares(endingOrder[i]); int index = indices[i]; if ((uint)index >= (uint)divisors.Length) { continue; } - ulong candidate = (ulong)divisors[index]; - if (candidate > sqrtLimit) + ulong candidateSquare = squares[index]; + if (candidateSquare > maxDivisorSquare) { indices[i] = divisors.Length; continue; } + ulong candidate = (ulong)divisors[index]; if (bestEndingIndex == -1 || candidate < bestCandidate) { bestCandidate = candidate; @@ -497,25 +458,29 @@ Cleanup: return false; } - divisor = bestCandidate; + ReadOnlySpan selectedDivisors = GetGroupBDivisors(endingOrder[bestEndingIndex]); + divisor = (ulong)selectedDivisors[indices[bestEndingIndex]]; indices[bestEndingIndex]++; return true; } + + + [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool HasGroupBCandidates(ReadOnlySpan endingOrder, Span indices, ulong sqrtLimit) + private static bool HasGroupBCandidates(ReadOnlySpan endingOrder, Span indices, ulong maxDivisorSquare) { for (int i = 0; i < endingOrder.Length; i++) { ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + ReadOnlySpan squares = GetGroupBDivisorSquares(endingOrder[i]); int index = indices[i]; if ((uint)index >= (uint)divisors.Length) { continue; } - ulong candidate = (ulong)divisors[index]; - if (candidate > sqrtLimit) + if (squares[index] > maxDivisorSquare) { indices[i] = divisors.Length; continue; @@ -527,6 +492,8 @@ Cleanup: return false; } + + private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch { 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, @@ -536,35 +503,32 @@ Cleanup: _ => ReadOnlySpan.Empty, }; - private static ReadOnlySpan GetWheel210ResiduesForEnding(byte ending) => ending switch + private static ReadOnlySpan GetGroupBDivisorSquares(byte ending) => ending switch { - 1 => Wheel210ResiduesEnding1, - 7 => Wheel210ResiduesEnding7, - 9 => Wheel210ResiduesEnding9, - _ => ReadOnlySpan.Empty, + 1 => DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree, + 3 => ReadOnlySpan.Empty, + 7 => DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree, + 9 => DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree, + _ => ReadOnlySpan.Empty, }; private static int InitializeGroupBStates(byte nMod10, Span buffer) { ReadOnlySpan endings = GetGroupBEndingOrder(nMod10); - int count = 0; + int count = endings.Length; - for (int priorityIndex = 0; priorityIndex < endings.Length; priorityIndex++) + for (int i = 0; i < count; i++) { - ReadOnlySpan residues = GetWheel210ResiduesForEnding(endings[priorityIndex]); - for (int residueIndex = 0; residueIndex < residues.Length; residueIndex++) - { - buffer[count++] = new HeuristicGroupBSequenceState(residues[residueIndex], endings[priorityIndex], (byte)priorityIndex); - } + buffer[i] = new HeuristicGroupBSequenceState(endings[i], (byte)i); } return count; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) { - return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); + return new HeuristicDivisorEnumerator(maxDivisorSquare, nMod10, groupBBuffer); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -640,158 +604,143 @@ Cleanup: return false; } - internal ref struct HeuristicDivisorEnumerator + +internal ref struct HeuristicDivisorEnumerator +{ + private readonly ulong maxDivisorSquare; + private readonly ReadOnlySpan groupADivisors; + private readonly ReadOnlySpan groupADivisorSquares; + private int groupAIndex; + private Span groupBStates; + + public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) { - private readonly ulong sqrtLimit; - private int groupAConstantIndex; - private ulong groupAWheelCandidate; - private int groupAIncrementIndex; - private bool groupAWheelActive; - private Span groupBStates; + this.maxDivisorSquare = maxDivisorSquare; + groupADivisors = HeuristicPrimeSieves.GroupADivisors; + groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; + groupAIndex = 0; - public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + int count = InitializeGroupBStates(nMod10, groupBBuffer); + groupBStates = count == 0 ? Span.Empty : groupBBuffer[..count]; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (groupAIndex < groupADivisors.Length) { - this.sqrtLimit = sqrtLimit; - groupAConstantIndex = 0; - groupAWheelCandidate = 23UL; - groupAIncrementIndex = 0; - groupAWheelActive = sqrtLimit >= 23UL; - - int count = InitializeGroupBStates(nMod10, groupBBuffer); - groupBStates = count == 0 ? Span.Empty : groupBBuffer[..count]; - - for (int i = 0; i < groupBStates.Length; i++) + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) { - groupBStates[i].TryNormalize(); - } - } - - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - // TODO: Hard-code constant A divisor checks here instead of using an array - while (groupAConstantIndex < GroupAConstantDivisors.Length) - { - ulong value = GroupAConstantDivisors[groupAConstantIndex++]; - if (value > sqrtLimit) - { - groupAConstantIndex = GroupAConstantDivisors.Length; - break; - } - - // TODO: We're constantly dividing constants mod 10 which is constant too. - candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAConstant, (byte)(value % 10UL), 0, (ushort)value); - return true; + groupAIndex = groupADivisors.Length; + break; } - if (groupAWheelActive) - { - ulong value = groupAWheelCandidate; - if (value <= sqrtLimit) - { - candidate = new HeuristicDivisorCandidate(value, HeuristicDivisorGroup.GroupAWheel, 3, 0, (ushort)(value % Wheel210)); - byte increment = GroupAIncrementPattern[groupAIncrementIndex]; - groupAIncrementIndex ^= 1; + int currentIndex = groupAIndex; + ulong divisor = (ulong)groupADivisors[currentIndex]; + groupAIndex++; - ulong next = value + increment; - if (next <= value) - { - groupAWheelActive = false; - } - else - { - groupAWheelCandidate = next; - } + HeuristicDivisorGroup group = currentIndex < GroupAConstantCount + ? HeuristicDivisorGroup.GroupAConstant + : HeuristicDivisorGroup.GroupAWheel; - return true; - } - - groupAWheelActive = false; - } - - ulong bestCandidate = ulong.MaxValue; - int bestIndex = -1; - - for (int i = 0; i < groupBStates.Length; i++) - { - ref HeuristicGroupBSequenceState state = ref groupBStates[i]; - ulong value = state.Candidate; - if (value == ulong.MaxValue) - { - continue; - } - - if (value > sqrtLimit) - { - state.Candidate = ulong.MaxValue; - continue; - } - - if (bestIndex == -1 || value < bestCandidate || - (value == bestCandidate && state.PriorityIndex < groupBStates[bestIndex].PriorityIndex)) - { - bestCandidate = value; - bestIndex = i; - } - } - - if (bestIndex == -1) - { - candidate = default; - return false; - } - - ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex]; - var result = new HeuristicDivisorCandidate(bestState.Candidate, HeuristicDivisorGroup.GroupB, bestState.Ending, bestState.PriorityIndex, bestState.Residue); - bestState.Advance(); - candidate = result; + candidate = new HeuristicDivisorCandidate( + divisor, + group, + (byte)(divisor % 10UL), + 0, + (ushort)(divisor % Wheel210)); return true; } - } - internal struct HeuristicGroupBSequenceState(ushort residue, byte ending, byte priorityIndex) - { - public ulong Candidate = residue; - public byte PriorityIndex = priorityIndex; - public byte Ending = ending; - public ushort Residue = residue; + ulong bestCandidate = ulong.MaxValue; + int bestIndex = -1; - public void TryNormalize() + for (int i = 0; i < groupBStates.Length; i++) { - ulong candidate = Candidate; - if (candidate == ulong.MaxValue) + ref HeuristicGroupBSequenceState state = ref groupBStates[i]; + if (!state.TryPeek(maxDivisorSquare, out ulong candidateValue)) { - return; + continue; } - while (candidate <= 13UL) + if (bestIndex == -1 || candidateValue < bestCandidate || + (candidateValue == bestCandidate && state.PriorityIndex < groupBStates[bestIndex].PriorityIndex)) { - ulong next = candidate + Wheel210; - if (next <= candidate) - { - Candidate = ulong.MaxValue; - return; - } - - candidate = next; + bestCandidate = candidateValue; + bestIndex = i; } - - Candidate = candidate; } - public void Advance() + if (bestIndex == -1) { - ulong candidate = Candidate; - if (candidate == ulong.MaxValue) - { - return; - } - - ulong next = candidate + Wheel210; - Candidate = next > candidate ? next : ulong.MaxValue; + candidate = default; + return false; } - + ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex]; + bestState.Advance(); + candidate = new HeuristicDivisorCandidate( + bestCandidate, + HeuristicDivisorGroup.GroupB, + bestState.Ending, + bestState.PriorityIndex, + (ushort)(bestCandidate % Wheel210)); + return true; } +} + + + +internal struct HeuristicGroupBSequenceState +{ + public byte Ending; + public byte PriorityIndex; + public int Index; + + public HeuristicGroupBSequenceState(byte ending, byte priorityIndex) + { + Ending = ending; + PriorityIndex = priorityIndex; + Index = 0; + } + + public bool TryPeek(ulong maxDivisorSquare, out ulong divisor) + { + ReadOnlySpan divisors = GetGroupBDivisors(Ending); + ReadOnlySpan squares = GetGroupBDivisorSquares(Ending); + int index = Index; + if ((uint)index >= (uint)divisors.Length) + { + divisor = 0UL; + return false; + } + + if (squares[index] > maxDivisorSquare) + { + Index = divisors.Length; + divisor = 0UL; + return false; + } + + divisor = (ulong)divisors[index]; + return true; + } + + public void Advance() + { + int next = Index + 1; + ReadOnlySpan divisors = GetGroupBDivisors(Ending); + if ((uint)next >= (uint)divisors.Length) + { + Index = divisors.Length; + } + else + { + Index = next; + } + } +} + internal ref struct MersenneHeuristicDivisorEnumerator { diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 851473c..149bc77 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -411,8 +411,8 @@ public class MersenneDivisorCycles // Atomic.Add(ref _tryFactorIntoCountsInternalHits, 1UL); // Console.WriteLine($"MersenneDivisorCycles.TryFactorIntoCountsInternal hits {Volatile.Read(ref _tryFactorIntoCountsInternalHits)}"); - HeuristicPrimeTester primeTester = _primeTester ??= new(); - if (primeTester.IsPrimeGpu(remaining)) + // HeuristicPrimeTester primeTester = _primeTester ??= new(); + if (PrimeTester.IsPrime(remaining)) // if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None)) { AddFactor(counts, remaining); @@ -917,8 +917,8 @@ public class MersenneDivisorCycles return CalculateCycleLengthFallback(divisor); } - [ThreadStatic] - private static HeuristicPrimeTester? _primeTester; + // [ThreadStatic] + // private static HeuristicPrimeTester? _primeTester; internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false) { @@ -940,10 +940,10 @@ public class MersenneDivisorCycles // Console.WriteLine($"MersenneDivisorCycles.CalculateCycleLengthFallback hits {Volatile.Read(ref _tryCalculateCycleLengthHeuristicHits)}"); // } - HeuristicPrimeTester primeTester = _primeTester ??= new(); + // HeuristicPrimeTester primeTester = _primeTester ??= new(); // if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None)) - if (!skipPrimeOrderHeuristic && primeTester.IsPrimeGpu(divisor)) - // if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeCpu(divisor, CancellationToken.None)) + // if (!skipPrimeOrderHeuristic && primeTester.IsPrimeGpu(divisor)) + if (!skipPrimeOrderHeuristic && PrimeTester.IsPrime(divisor)) { ulong computedOrder = PrimeOrderCalculator.Calculate( divisor, diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs index a160a3b..93dcb58 100644 --- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs +++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs @@ -34,7 +34,7 @@ public static class MersennePrimeFactorTester // Atomic.Add(ref _isPrimeFactorHits, 1UL); // Console.WriteLine($"MersennePrimeFactorTester.IsPrimeFactor hits {Volatile.Read(ref _isPrimeFactorHits)}"); - if (!PrimeTester.IsPrimeCpu(q, ct)) + if (!PrimeTester.IsPrime(q)) { return false; } @@ -247,7 +247,7 @@ public static class MersennePrimeFactorTester // Console.WriteLine($"MersennePrimeFactorTester.Factor64 hits {Volatile.Read(ref _factor64Hits)}"); // if (HeuristicPrimeTester.Exclusive.IsPrime(m, ct)) - if (PrimeTester.IsPrimeCpu(m, ct)) + if (PrimeTester.IsPrime(m)) { dict[m] = 1; continue; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 1be9074..3f5773e 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -1079,7 +1079,7 @@ internal static partial class PrimeOrderCalculator PartialFactorResult result; PrimeTester primeTester = PrimeTester.Exclusive; CancellationToken ct = CancellationToken.None; - HeuristicPrimeTester? tester = _tester ??= new(); + // HeuristicPrimeTester? tester = _tester ??= new(); if (!gpuFactored) { @@ -1152,7 +1152,7 @@ internal static partial class PrimeOrderCalculator // Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}"); // bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); - bool isPrime = tester.IsPrimeGpu(composite); + bool isPrime = PrimeTester.IsPrime(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) @@ -1220,7 +1220,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 = tester.IsPrimeGpu(composite); + bool isPrime = PrimeTester.IsPrime(composite); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); @@ -1257,7 +1257,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorCofactorHits, 1UL); // Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}"); - cofactorIsPrime = tester.IsPrimeGpu(cofactor); + cofactorIsPrime = PrimeTester.IsPrime(cofactor); // cofactorIsPrime = PrimeTester.IsPrimeGpu(cofactor); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1782,8 +1782,8 @@ internal static partial class PrimeOrderCalculator // private static ulong _factorCompletelyHits; // private static ulong _factorCompletelyPollardRhoHits; - [ThreadStatic] - private static HeuristicPrimeTester? _tester; + // [ThreadStatic] + // private static HeuristicPrimeTester? _tester; private static void FactorCompletely(ulong value, Dictionary counts, bool knownComposite) { @@ -1794,10 +1794,10 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _factorCompletelyHits, 1UL); // Console.WriteLine($"Factor completely hits {Volatile.Read(ref _factorCompletelyHits)}"); - HeuristicPrimeTester tester = _tester ??= new(); + // HeuristicPrimeTester tester = _tester ??= new(); // bool isPrime = primeTester.HeuristicIsPrimeGpu(value); - bool isPrime = tester.IsPrimeGpu(value); + bool isPrime = PrimeTester.IsPrime(value); // bool isPrime = PrimeTester.IsPrimeGpu(value); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value); if (!knownComposite && isPrime) @@ -1813,7 +1813,7 @@ internal static partial class PrimeOrderCalculator // Console.WriteLine($"Factor completely after PollardRho hits {Volatile.Read(ref _factorCompletelyPollardRhoHits)}"); // isPrime = primeTester.HeuristicIsPrimeGpu(factor); - isPrime = tester.IsPrimeGpu(factor); + isPrime = PrimeTester.IsPrime(factor); // isPrime = PrimeTester.IsPrimeGpu(factor); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor); if (isPrime) @@ -1850,7 +1850,7 @@ internal static partial class PrimeOrderCalculator // Atomic.Add(ref _partialFactorCofactorHits, 1UL); // Console.WriteLine($"Partial factor cofactor hits {Volatile.Read(ref _partialFactorCofactorHits)}"); - isPrime = tester.IsPrimeGpu(quotient); + isPrime = PrimeTester.IsPrime(quotient); // isPrime = PrimeTester.IsPrimeGpu(quotient); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient); } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs index 2f923c6..2fe2791 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -950,7 +950,7 @@ internal static partial class PrimeOrderCalculator if (value <= ulong.MaxValue) { // return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None); - return PrimeTester.IsPrimeCpu((ulong)value, CancellationToken.None); + return PrimeTester.IsPrime((ulong)value); } if ((value & UInt128.One) == UInt128.Zero) diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 7c4e3af..52bb5bf 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -2,6 +2,7 @@ using System.Buffers; using System.Collections.Concurrent; using System.Numerics; using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; @@ -19,7 +20,7 @@ public sealed class PrimeTester public static PrimeTester Exclusive => _tester ??= new(); - public static bool IsPrime(ulong n, CancellationToken ct) + public static bool IsPrime(ulong n) { if (n <= 1UL) { @@ -89,56 +90,70 @@ public sealed class PrimeTester return result; } - public static bool IsPrimeGpu(ulong n) - { - return Exclusive.IsPrimeGpu(n, CancellationToken.None); - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool IsPrimeGpu(ulong n, CancellationToken ct) + public bool IsPrimeGpu(ulong n) { bool forceCpu = GpuContextPool.ForceCpu; bool belowGpuRange = n < 31UL; if (forceCpu || belowGpuRange) { - return IsPrimeCpu(n, ct); + return IsPrime(n); } - Span one = stackalloc ulong[1]; - Span outFlags = stackalloc byte[1]; - one[0] = n; - outFlags[0] = 0; + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTesterGpuContextPool.Rent(1); + var accelerator = gpu.Accelerator; + var kernel = gpu.Kernel; - IsPrimeBatchGpu(one, outFlags); + ulong value = n; + byte flag = 0; - return outFlags[0] != 0; + var input = gpu.Input; + var output = gpu.Output; + var inputView = input.View; + var outputView = output.View; + + 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; } - 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 gpu = PrimeTesterGpuContextPool.Rent(GpuBatchSize); var accelerator = gpu.Accelerator; + var kernel = gpu.Kernel; int totalLength = values.Length; - int batchSize = Math.Max(1, GpuBatchSize); + int batchSize = GpuBatchSize; - var scratch = state.RentScratch(batchSize, accelerator); - var input = scratch.Input; - var output = scratch.Output; - ArrayPool pool = ThreadStaticPools.UlongPool; - ulong[] temp = pool.Rent(batchSize); + var input = gpu.Input; + var output = gpu.Output; + var inputView = input.View; + var outputView = output.View; int pos = 0; while (pos < totalLength) @@ -146,290 +161,230 @@ public sealed class PrimeTester int remaining = totalLength - pos; int count = remaining > batchSize ? batchSize : remaining; - values.Slice(pos, count).CopyTo(temp); - input.View.CopyFromCPU(ref temp[0], count); + var valueSlice = values.Slice(pos, count); + ref ulong valueRef = ref MemoryMarshal.GetReference(valueSlice); + inputView.CopyFromCPU(ref valueRef, 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); + kernel( + count, + 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(); - output.View.CopyToCPU(ref results[pos], count); + var resultSlice = results.Slice(pos, count); + ref byte resultRef = ref MemoryMarshal.GetReference(resultSlice); + outputView.CopyToCPU(ref resultRef, count); pos += count; } - pool.Return(temp, clearArray: false); - state.ReturnScratch(scratch); - gpu.Dispose(); limiter.Dispose(); } internal static class PrimeTesterGpuContextPool { - internal sealed class PooledContext + private static readonly ConcurrentQueue Pool = new(); + + internal static PrimeTesterGpuContextLease Rent(int minBufferCapacity = 1) { - public Context Context { get; } - - public Accelerator Accelerator { get; } - - private KernelState? _kernelState; - - public PooledContext() + if (!Pool.TryDequeue(out var lease)) { - Context = Context.CreateDefault(); - Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context); - _kernelState = GpuKernelState.GetOrCreate(Accelerator); + lease = new PrimeTesterGpuContextLease(minBufferCapacity); + } + else + { + lease.EnsureCapacity(minBufferCapacity); } - 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(); - } + lease.ResetReturnFlag(); + return lease; } - private static readonly ConcurrentQueue Pool = new(); - - internal static PrimeTesterGpuContextLease Rent() + internal static void Return(PrimeTesterGpuContextLease lease) { - if (Pool.TryDequeue(out var ctx)) - { - return new PrimeTesterGpuContextLease(ctx); - } - - return new PrimeTesterGpuContextLease(new PooledContext()); + lease.Accelerator.Synchronize(); + Pool.Enqueue(lease); } internal static void DisposeAll() { - while (Pool.TryDequeue(out var ctx)) + while (Pool.TryDequeue(out var lease)) { - ctx.Dispose(); + lease.DisposeResources(); } } - private static void Return(PooledContext ctx) + internal static void Clear(Accelerator accelerator) { - 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) + if (accelerator is null) { return; } - while (_scratchPool.TryTake(out var sb)) + var retained = new List(); + while (Pool.TryDequeue(out var lease)) { - 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) + if (lease.Accelerator == accelerator) { - var state = lazy.Value; - state.Clear(); + lease.DisposeResources(); + } + else + { + retained.Add(lease); } } + + int retainedCount = retained.Count; + for (int i = 0; i < retainedCount; i++) + { + Pool.Enqueue(retained[i]); + } + } + + internal sealed class PrimeTesterGpuContextLease : IDisposable + { + private bool _returned; + + internal readonly Context AcceleratorContext; + public readonly Accelerator Accelerator; + public readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> Kernel; + public readonly Action, ulong, ArrayView> HeuristicTrialDivisionKernel; + public readonly MemoryBuffer1D DevicePrimesDefault; + public readonly MemoryBuffer1D DevicePrimesLastOne; + public readonly MemoryBuffer1D DevicePrimesLastSeven; + public readonly MemoryBuffer1D DevicePrimesLastThree; + public readonly MemoryBuffer1D DevicePrimesLastNine; + public readonly MemoryBuffer1D DevicePrimesPow2Default; + public readonly MemoryBuffer1D DevicePrimesPow2LastOne; + public readonly MemoryBuffer1D DevicePrimesPow2LastSeven; + public readonly MemoryBuffer1D DevicePrimesPow2LastThree; + public readonly MemoryBuffer1D DevicePrimesPow2LastNine; + public MemoryBuffer1D Input; + public MemoryBuffer1D Output; + public MemoryBuffer1D HeuristicFlag; + public int BufferCapacity; + + internal PrimeTesterGpuContextLease(int minBufferCapacity) + { + AcceleratorContext = Context.CreateDefault(); + Accelerator = AcceleratorContext.GetPreferredDevice(false).CreateAccelerator(AcceleratorContext); + 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); + + BufferCapacity = minBufferCapacity < 1 ? 1 : minBufferCapacity; + Input = Accelerator.Allocate1D(BufferCapacity); + Output = Accelerator.Allocate1D(BufferCapacity); + HeuristicFlag = Accelerator.Allocate1D(BufferCapacity); + _returned = false; + } + + public void EnsureCapacity(int minCapacity) + { + int required = minCapacity < 1 ? 1 : minCapacity; + if (required <= BufferCapacity) + { + return; + } + + Input.Dispose(); + Output.Dispose(); + HeuristicFlag.Dispose(); + Input = Accelerator.Allocate1D(required); + Output = Accelerator.Allocate1D(required); + HeuristicFlag = Accelerator.Allocate1D(required); + BufferCapacity = required; + } + + internal void ResetReturnFlag() + { + _returned = false; + } + + public void Dispose() + { + if (_returned) + { + return; + } + + _returned = true; + PrimeTesterGpuContextPool.Return(this); + } + + internal void DisposeResources() + { + Input.Dispose(); + Output.Dispose(); + HeuristicFlag.Dispose(); + DevicePrimesDefault.Dispose(); + DevicePrimesLastOne.Dispose(); + DevicePrimesLastSeven.Dispose(); + DevicePrimesLastThree.Dispose(); + DevicePrimesLastNine.Dispose(); + DevicePrimesPow2Default.Dispose(); + DevicePrimesPow2LastOne.Dispose(); + DevicePrimesPow2LastSeven.Dispose(); + DevicePrimesPow2LastThree.Dispose(); + DevicePrimesPow2LastNine.Dispose(); + Accelerator.Dispose(); + AcceleratorContext.Dispose(); + } } } - // Expose cache clearing for accelerator disposal coordination public static void ClearGpuCaches(Accelerator accelerator) { - GpuKernelState.Clear(accelerator); + PrimeTesterGpuContextPool.Clear(accelerator); } internal static void DisposeGpuContexts() diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md index e2157c9..cc024cd 100644 --- a/docs/HeuristicIsPrimePlan.md +++ b/docs/HeuristicIsPrimePlan.md @@ -131,7 +131,7 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, * [done] The CPU heuristic now operates as a pure trial-division sweep without maintaining per-call summaries, relying on `PrepareHeuristicDivisor` only when downstream consumers request Montgomery data. * [done] Introduced the temporary `UseHeuristicGroupBTrialDivision` gate so current builds execute Group A locally and then fall back to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, keeping the full Group B implementation available for reactivation. 7. [done] **Implement `IsPrimeGpu`** - * [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `sqrtLimit`/`nMod10` inputs, and falling back to CPU when GPUs are disabled. + * [done] Mirror CPU structure while batching divisor checks on the accelerator via `HeuristicTrialDivisionGpu`, drawing divisors from `HeuristicPrimeSieves`, reusing caller-provided `maxDivisorSquare`/`nMod10` inputs, and falling back to CPU when GPUs are disabled. * [done] Added `PrimeTesterKernels.HeuristicTrialDivisionKernel` so GPU batches respect Group A/B ordering while emitting hit flags for early exits. * [done] Kept GPU divisibility in standard modular arithmetic while reserving Montgomery transforms for CPU confirmation logic. * [done] When the temporary Group B gate is disabled, delegate GPU calls to the CPU fallback so Group A coverage remains in place while the accelerator path stays available for future re-enablement.