diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index e1e9f9c..981de9b 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -12,8 +12,8 @@ namespace EvenPerfectBitScanner.Benchmarks; [MemoryDiagnoser] public class PrimeTesterBenchmarks { - private PrimeTester _cpuTester = null!; - private PrimeTester _gpuTester = null!; + private HeuristicPrimeTester _cpuTester = null!; + private HeuristicPrimeTester _gpuTester = null!; private ulong[] _candidates = []; [ParamsSource(nameof(GetBenchmarkCases))] @@ -22,9 +22,8 @@ public class PrimeTesterBenchmarks [GlobalSetup] public void GlobalSetup() { - PrimeTester.EnableHeuristicPrimeTesting = true; - _cpuTester = new PrimeTester(); - _gpuTester = new PrimeTester(); + _cpuTester = new HeuristicPrimeTester(); + _gpuTester = new HeuristicPrimeTester(); _candidates = BenchmarkCase.Candidates; if (_candidates.Length == 0) @@ -41,7 +40,7 @@ public class PrimeTesterBenchmarks [Benchmark(Baseline = true)] public int HeuristicCpu() { - PrimeTester tester = _cpuTester; + HeuristicPrimeTester tester = _cpuTester; ulong[] values = _candidates; CancellationToken cancellationToken = CancellationToken.None; int primeCount = 0; @@ -60,14 +59,14 @@ public class PrimeTesterBenchmarks [Benchmark] public int HeuristicGpu() { - PrimeTester tester = _gpuTester; + HeuristicPrimeTester tester = _gpuTester; ulong[] values = _candidates; CancellationToken cancellationToken = CancellationToken.None; int primeCount = 0; for (int i = 0; i < values.Length; i++) { - if (tester.HeuristicIsPrimeGpu(values[i])) + if (tester.IsPrimeGpu(values[i])) { primeCount++; } @@ -87,7 +86,7 @@ public class PrimeTesterBenchmarks for (int i = 0; i < values.Length; i++) { - if (PrimeTester.LegacyIsPrimeGpu(values[i], cancellationToken)) + if (PrimeTester.IsPrimeGpuFallback(values[i], cancellationToken)) { primeCount++; } @@ -106,7 +105,7 @@ public class PrimeTesterBenchmarks for (int i = 0; i < values.Length; i++) { - if (PrimeTester.LegacyIsPrimeInternal(values[i], cancellationToken)) + if (PrimeTester.IsPrimeInternal(values[i], cancellationToken)) { primeCount++; } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index e4f4924..c061da1 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -106,13 +106,17 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv processedAll = false; processedCount = 0UL; + // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); + // while (enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate)) while (enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate)) { processedCount = enumerator.ProcessedCount; + // HeuristicPrimeTester.HeuristicDivisorPreparation preparation = HeuristicPrimeTester.PrepareHeuristicDivisor(in candidate); PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate); + // ulong divisorCycle = HeuristicPrimeTester.ResolveHeuristicCycleLength( ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength( prime, in preparation, diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index 9e9b729..e1bc1de 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -189,6 +189,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv processedCount = 0UL; lastProcessed = 0UL; + // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); Span divisorSpan = divisors.AsSpan(); @@ -203,12 +204,15 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv int admissibleCount = 0; int exponentIndex = 0; + // while (admissibleCount < batchCapacity && enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate)) while (admissibleCount < batchCapacity && enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate)) { processedCount = enumerator.ProcessedCount; ulong divisorValue = candidate.Value; + // HeuristicPrimeTester.HeuristicDivisorPreparation preparation = HeuristicPrimeTester.PrepareHeuristicDivisor(in candidate); PrimeTester.HeuristicDivisorPreparation preparation = PrimeTester.PrepareHeuristicDivisor(in candidate); + // ulong divisorCycle = HeuristicPrimeTester.ResolveHeuristicCycleLength( ulong divisorCycle = PrimeTester.ResolveHeuristicCycleLength( prime, in preparation, diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs new file mode 100644 index 0000000..d5d70a6 --- /dev/null +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -0,0 +1,968 @@ +using System; +using System.Buffers; +using System.Collections.Concurrent; +using Open.Numeric.Primes; +using System.Collections.Generic; +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Threading; +using ILGPU; +using ILGPU.Algorithms; +using ILGPU.Runtime; +using PerfectNumbers.Core.Gpu; + +namespace PerfectNumbers.Core; + +public sealed class HeuristicPrimeTester +{ + static HeuristicPrimeTester() + { + HeuristicPrimeSieves.EnsureInitialized(); + } + + // 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]; + private static readonly byte[] GroupBEndingOrderMod9 = [9, 7, 1]; + + // 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; + internal const int MaxGroupBSequences = 36; + + private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. + + + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + + 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 bool EnableHeuristicPrimeTesting { get; set; } = true; + + public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096; + + public static int HeuristicDivisorInterleaveBatchSize { get; set; } = 64; + + [ThreadStatic] + private static HeuristicPrimeTester? _tester; + + public static HeuristicPrimeTester Exclusive => _tester ??= new(); + + public bool IsPrime(ulong n, CancellationToken ct) + { + if (!EnableHeuristicPrimeTesting) + { + return PrimeTester.IsPrimeInternal(n, ct); + } + + return HeuristicIsPrimeCpu(n); + } + + public bool IsPrime(ulong n) + { + return IsPrime(n, CancellationToken.None); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool IsPrimeGpu(ulong n, CancellationToken ct) + { + if (!EnableHeuristicPrimeTesting) + { + return PrimeTester.IsPrimeGpuFallback(n, ct); + } + + return IsPrimeGpu(n); + } + + 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) + { + if (!EnableHeuristicPrimeTesting) + { + return PrimeTester.IsPrimeGpuFallback(n, CancellationToken.None); + } + + return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10); + } + + public bool HeuristicIsPrimeGpu(ulong n) + { + return IsPrimeGpu(n); + } + + public static bool HeuristicIsPrimeCpu(ulong n) + { + return HeuristicIsPrimeCpu(n, 0UL, (byte)n.Mod10()); + } + + public static bool HeuristicIsPrimeCpu(ulong n, ulong sqrtLimit, byte nMod10) + { + if (!EnableHeuristicPrimeTesting) + { + return PrimeTester.IsPrimeInternal(n, CancellationToken.None); + } + + if (sqrtLimit == 0UL) + { + sqrtLimit = ComputeHeuristicSqrt(n); + } + + if (sqrtLimit < 3UL) + { + return EvaluateWithOpenNumericFallback(n); + } + + bool includeGroupB = UseHeuristicGroupBTrialDivision; + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB); + } + + 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, includeGroupB: false); + } + + if (GpuContextPool.ForceCpu) + { + return HeuristicIsPrimeCpu(n, sqrtLimit, nMod10); + } + + bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); + return !compositeDetected; + } + + private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB) + { + ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; + int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); + int groupAIndex = 0; + + if (!includeGroupB) + { + while (groupAIndex < groupADivisors.Length) + { + int processed = 0; + while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) + { + ulong divisor = (ulong)groupADivisors[groupAIndex]; + if (divisor > sqrtLimit) + { + groupAIndex = groupADivisors.Length; + break; + } + + if (n % divisor == 0UL) + { + return false; + } + + groupAIndex++; + processed++; + } + + if (groupAIndex >= groupADivisors.Length) + { + break; + } + + if ((ulong)groupADivisors[groupAIndex] > sqrtLimit) + { + groupAIndex = groupADivisors.Length; + break; + } + } + + return EvaluateWithOpenNumericFallback(n); + } + + ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); + if (endingOrder.IsEmpty) + { + return EvaluateWithOpenNumericFallback(n); + } + + Span indices = endingOrder.Length <= 8 + ? stackalloc int[endingOrder.Length] + : new int[endingOrder.Length]; + + bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; + bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); + + while (groupAHasMore || groupBHasMore) + { + if (groupAHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) + { + ulong divisor = (ulong)groupADivisors[groupAIndex]; + if (divisor > sqrtLimit) + { + groupAIndex = groupADivisors.Length; + groupAHasMore = false; + break; + } + + if (n % divisor == 0UL) + { + return false; + } + + groupAIndex++; + processed++; + } + + if (groupAIndex >= groupADivisors.Length) + { + groupAHasMore = false; + } + else + { + groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; + if (!groupAHasMore) + { + groupAIndex = groupADivisors.Length; + } + } + } + + if (groupBHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize) + { + if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) + { + groupBHasMore = false; + break; + } + + if (n % divisor == 0UL) + { + return false; + } + + processed++; + } + + if (groupBHasMore) + { + groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); + } + } + } + + return true; + } + + + + private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) + { + int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); + int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); + var divisorPool = ThreadStaticPools.UlongPool; + var hitPool = ThreadStaticPools.BytePool; + + ulong[]? divisorArray = null; + byte[]? hitFlags = null; + + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(); + + bool compositeDetected = false; + + try + { + divisorArray = divisorPool.Rent(batchCapacity); + hitFlags = hitPool.Rent(batchCapacity); + + var accelerator = gpu.Accelerator; + var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); + + lock (gpu.ExecutionLock) + { + var scratch = state.RentScratch(batchCapacity, accelerator); + try + { + int count = 0; + + bool ProcessBatch(int length) + { + if (length <= 0) + { + return false; + } + + scratch.Input.View.CopyFromCPU(ref divisorArray![0], length); + state.HeuristicTrialDivisionKernel(length, scratch.Input.View, n, scratch.Output.View); + accelerator.Synchronize(); + scratch.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; + } + + bool QueueDivisor(ulong divisor) + { + divisorArray![count] = divisor; + count++; + + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + return true; + } + + count = 0; + } + + 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]; + + bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; + bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, sqrtLimit); + + while (!compositeDetected && (groupAHasMore || groupBHasMore)) + { + if (groupAHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) + { + ulong divisor = (ulong)groupADivisors[groupAIndex]; + if (divisor > sqrtLimit) + { + groupAIndex = groupADivisors.Length; + groupAHasMore = false; + break; + } + + compositeDetected = QueueDivisor(divisor); + if (compositeDetected) + { + break; + } + + groupAIndex++; + processed++; + } + + if (compositeDetected) + { + break; + } + + if (groupAIndex >= groupADivisors.Length) + { + groupAHasMore = false; + } + else + { + groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; + if (!groupAHasMore) + { + groupAIndex = groupADivisors.Length; + } + } + } + + if (groupBHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize) + { + if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) + { + groupBHasMore = false; + break; + } + + compositeDetected = QueueDivisor(divisor); + if (compositeDetected) + { + break; + } + + processed++; + } + + if (compositeDetected) + { + break; + } + + if (groupBHasMore) + { + groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); + } + } + } + + if (!compositeDetected && count > 0) + { + compositeDetected = ProcessBatch(count); + } + } + finally + { + state.ReturnScratch(scratch); + } + } + } + finally + { + if (divisorArray is not null) + { + divisorPool.Return(divisorArray); + } + + if (hitFlags is not null) + { + 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 GetGroupBEndingOrder(byte nMod10) => nMod10 switch + { + 1 => GroupBEndingOrderMod1, + 3 => GroupBEndingOrderMod3, + 7 => GroupBEndingOrderMod7, + 9 => GroupBEndingOrderMod9, + _ => ReadOnlySpan.Empty, + }; + + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool TrySelectNextGroupBDivisor(ReadOnlySpan endingOrder, Span indices, ulong sqrtLimit, out ulong divisor) + { + ulong bestCandidate = ulong.MaxValue; + int bestEndingIndex = -1; + + for (int i = 0; i < endingOrder.Length; i++) + { + ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + int index = indices[i]; + if ((uint)index >= (uint)divisors.Length) + { + continue; + } + + ulong candidate = (ulong)divisors[index]; + if (candidate > sqrtLimit) + { + indices[i] = divisors.Length; + continue; + } + + if (bestEndingIndex == -1 || candidate < bestCandidate) + { + bestCandidate = candidate; + bestEndingIndex = i; + } + } + + if (bestEndingIndex == -1) + { + divisor = 0UL; + return false; + } + + divisor = bestCandidate; + indices[bestEndingIndex]++; + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool HasGroupBCandidates(ReadOnlySpan endingOrder, Span indices, ulong sqrtLimit) + { + for (int i = 0; i < endingOrder.Length; i++) + { + ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + int index = indices[i]; + if ((uint)index >= (uint)divisors.Length) + { + continue; + } + + ulong candidate = (ulong)divisors[index]; + if (candidate > sqrtLimit) + { + indices[i] = divisors.Length; + continue; + } + + return true; + } + + return false; + } + + private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch + { + 1 => HeuristicPrimeSieves.GroupBDivisorsEnding1, + 3 => HeuristicPrimeSieves.GroupBDivisorsEnding3, + 7 => HeuristicPrimeSieves.GroupBDivisorsEnding7, + 9 => HeuristicPrimeSieves.GroupBDivisorsEnding9, + _ => ReadOnlySpan.Empty, + }; + + private static ReadOnlySpan GetWheel210ResiduesForEnding(byte ending) => ending switch + { + 1 => Wheel210ResiduesEnding1, + 7 => Wheel210ResiduesEnding7, + 9 => Wheel210ResiduesEnding9, + _ => ReadOnlySpan.Empty, + }; + + private static int InitializeGroupBStates(byte nMod10, Span buffer) + { + ReadOnlySpan endings = GetGroupBEndingOrder(nMod10); + int count = 0; + + for (int priorityIndex = 0; priorityIndex < endings.Length; priorityIndex++) + { + ReadOnlySpan residues = GetWheel210ResiduesForEnding(endings[priorityIndex]); + for (int residueIndex = 0; residueIndex < residues.Length; residueIndex++) + { + buffer[count++] = new HeuristicGroupBSequenceState(residues[residueIndex], endings[priorityIndex], (byte)priorityIndex); + } + } + + return count; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); + } + + [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 int groupAConstantIndex; + private ulong groupAWheelCandidate; + private int groupAIncrementIndex; + private bool groupAWheelActive; + private Span groupBStates; + + public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + 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++) + { + 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; + } + + 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; + + ulong next = value + increment; + if (next <= value) + { + groupAWheelActive = false; + } + else + { + groupAWheelCandidate = next; + } + + 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; + 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; + + public void TryNormalize() + { + ulong candidate = Candidate; + if (candidate == ulong.MaxValue) + { + return; + } + + while (candidate <= 13UL) + { + ulong next = candidate + Wheel210; + if (next <= candidate) + { + Candidate = ulong.MaxValue; + return; + } + + candidate = next; + } + + Candidate = candidate; + } + + public void Advance() + { + ulong candidate = Candidate; + if (candidate == ulong.MaxValue) + { + return; + } + + ulong next = candidate + Wheel210; + Candidate = next > candidate ? next : ulong.MaxValue; + } + + + } + + 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/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 772fab8..4988284 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -395,6 +395,7 @@ public class MersenneDivisorCycles } // TODO: Benchmark PrimeTester.IsPrimeInternal against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. + // if (HeuristicPrimeTester.Exclusive.IsPrime(remaining, CancellationToken.None)) if (PrimeTester.IsPrimeInternal(remaining, CancellationToken.None)) { AddFactor(counts, remaining); @@ -927,6 +928,7 @@ public class MersenneDivisorCycles return true; } + // if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrime(divisor, CancellationToken.None)) if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) { ulong computedOrder = PrimeOrderCalculator.Calculate( diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs index ecf9022..be53488 100644 --- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs +++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs @@ -25,6 +25,7 @@ public static class MersennePrimeFactorTester return false; } + // if (!HeuristicPrimeTester.Exclusive.IsPrime(q, ct)) if (!PrimeTester.IsPrimeInternal(q, ct)) { return false; @@ -233,6 +234,7 @@ public static class MersennePrimeFactorTester continue; } + // if (HeuristicPrimeTester.Exclusive.IsPrime(m, ct)) if (PrimeTester.IsPrimeInternal(m, ct)) { dict[m] = 1; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index b7a8665..e24d625 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -1124,7 +1124,7 @@ internal static partial class PrimeOrderCalculator // bool isPrime = primeTester.HeuristicIsPrimeGpu(composite); // bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); - bool isPrime = PrimeTester.LegacyIsPrimeInternal(composite, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) @@ -1187,7 +1187,7 @@ internal static partial class PrimeOrderCalculator { // bool isPrime = primeTester.HeuristicIsPrimeGpu(composite); // bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); - bool isPrime = PrimeTester.LegacyIsPrimeInternal(composite, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); @@ -1220,7 +1220,7 @@ internal static partial class PrimeOrderCalculator { // cofactorIsPrime = primeTester.HeuristicIsPrimeGpu(cofactor); // cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None); - cofactorIsPrime = PrimeTester.LegacyIsPrimeInternal(cofactor, CancellationToken.None); + cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1753,7 +1753,7 @@ internal static partial class PrimeOrderCalculator // bool isPrime = primeTester.HeuristicIsPrimeGpu(value); // bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None); - bool isPrime = PrimeTester.LegacyIsPrimeInternal(value, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value); if (!knownComposite && isPrime) { @@ -1766,7 +1766,7 @@ internal static partial class PrimeOrderCalculator // isPrime = primeTester.HeuristicIsPrimeGpu(factor); // isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None); - isPrime = PrimeTester.LegacyIsPrimeInternal(factor, CancellationToken.None); + isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor); if (isPrime) { @@ -1798,7 +1798,7 @@ internal static partial class PrimeOrderCalculator { // isPrime = primeTester.HeuristicIsPrimeGpu(quotient); // isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None); - isPrime = PrimeTester.LegacyIsPrimeInternal(quotient, CancellationToken.None); + isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(quotient); } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs index d385fb8..bd81fbb 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -955,6 +955,7 @@ internal static partial class PrimeOrderCalculator { if (value <= ulong.MaxValue) { + // return HeuristicPrimeTester.Exclusive.IsPrime((ulong)value, CancellationToken.None); return PrimeTester.IsPrimeInternal((ulong)value, CancellationToken.None); } diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 5676956..9ec21b5 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -1,10 +1,9 @@ using System; using System.Buffers; using System.Collections.Concurrent; -using Open.Numeric.Primes; -using System.Collections.Generic; using System.Numerics; using System.Runtime.CompilerServices; +using System.Threading; using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; @@ -12,73 +11,10 @@ using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; -public sealed class PrimeTester(bool useInternal = false) +public sealed class PrimeTester { - static PrimeTester() + public PrimeTester(bool useInternal = false) { - HeuristicPrimeSieves.EnsureInitialized(); - } - - // 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]; - private static readonly byte[] GroupBEndingOrderMod9 = [9, 7, 1]; - - // 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; - internal const int MaxGroupBSequences = 36; - - private const bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. - - private readonly bool _useLegacyPrimeTester = useInternal; - - private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - - 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 bool IsPrime(ulong n, CancellationToken ct) - { - if (_useLegacyPrimeTester) - { - return LegacyIsPrimeInternal(n, ct); - } - - return IsPrimeInternal(n, ct); } [ThreadStatic] @@ -86,6 +22,11 @@ public sealed class PrimeTester(bool useInternal = false) public static PrimeTester Exclusive => _tester ??= new(); + public bool IsPrime(ulong n, CancellationToken ct) + { + return IsPrimeInternal(n, ct); + } + public static bool IsPrimeGpu(ulong n) { return Exclusive.IsPrimeGpu(n, CancellationToken.None); @@ -93,48 +34,20 @@ public sealed class PrimeTester(bool useInternal = false) public static bool IsPrimeGpu(ulong n, ulong limit, byte nMod10) { - return Exclusive.HeuristicIsPrimeGpu(n, limit, nMod10); + return Exclusive.IsPrimeGpu(n, CancellationToken.None); } - // Optional GPU-assisted primality: batched small-prime sieve on device. [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool IsPrimeGpu(ulong n, CancellationToken ct) { - if (_useLegacyPrimeTester || !EnableHeuristicPrimeTesting) - { - return LegacyIsPrimeGpu(n, ct); - } - - return HeuristicIsPrimeGpu(n); + return IsPrimeGpuFallback(n, ct); } - public static bool LegacyIsPrimeGpu(ulong n, CancellationToken ct) + public static bool IsPrimeGpuFallback(ulong n, CancellationToken ct) { - // EvenPerfectBitScanner never enqueues exponents <= 3 on the GPU path; keep the legacy guard commented out - // so ad-hoc callers that expect that behavior must re-enable it intentionally. - // if (n <= 3UL) - // { - // return n >= 2UL; - // } - - // Candidate generation (Add/AddPrimes transforms plus ModResidueTracker) already strips even values and multiples of five. - // The downstream small-factor sweep and Pollard routines factor divisors q = 2kp + 1 with k >= 1 without feeding new - // exponents back into this path, so the GPU sees the same odd sequence that already passed the CPU residues. - - // EvenPerfectBitScanner seeds p at 136,279,841 and advances monotonically, so production scans never hit the GPU with - // exponents below 31. Leave this CPU redirect commented out to keep the wrapper branchless while documenting the guard. - // The ternary fallback below still routes sub-31 values through LegacyIsPrimeInternal for correctness when ad-hoc callers bypass - // the scanner pipeline. - // if (n < 31UL) - // { - // return LegacyIsPrimeInternal(n, ct); - // } - bool forceCpu = GpuContextPool.ForceCpu; Span one = stackalloc ulong[1]; Span outFlags = stackalloc byte[1]; - // TODO: Inline the single-value GPU sieve fast path from GpuModularArithmeticBenchmarks so this wrapper - // can skip stackalloc buffers and reuse the pinned upload span the benchmark identified as fastest. one[0] = n; outFlags[0] = 0; @@ -147,866 +60,13 @@ public sealed class PrimeTester(bool useInternal = false) bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0; bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime; - return requiresCpuFallback ? LegacyIsPrimeInternal(n, ct) : true; + return requiresCpuFallback ? IsPrimeInternal(n, ct) : true; } - public bool HeuristicIsPrimeGpu(ulong n) + public static bool IsPrimeInternal(ulong n, CancellationToken ct) { - byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); - return HeuristicIsPrimeGpu(n, sqrtLimit, nMod10); - } - - private bool HeuristicIsPrimeGpu(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, includeGroupB: false); - } - - if (GpuContextPool.ForceCpu) - { - return HeuristicIsPrimeCpu(n, sqrtLimit, nMod10); - } - - bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); - return !compositeDetected; - } - - internal static bool HeuristicIsPrimeCpu(ulong n) - { - return HeuristicIsPrimeCpu(n, 0UL, (byte)n.Mod10()); - } - - internal static bool HeuristicIsPrimeCpu(ulong n, ulong sqrtLimit, byte nMod10) - { - if (sqrtLimit == 0UL) - { - sqrtLimit = ComputeHeuristicSqrt(n); - } - - if (sqrtLimit < 3UL) - { - return EvaluateWithOpenNumericFallback(n); - } - - bool includeGroupB = UseHeuristicGroupBTrialDivision; - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB); - } - - private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB) - { - ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; - int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); - int groupAIndex = 0; - - if (!includeGroupB) - { - while (groupAIndex < groupADivisors.Length) - { - int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) - { - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (divisor > sqrtLimit) - { - groupAIndex = groupADivisors.Length; - break; - } - - if (n % divisor == 0UL) - { - return false; - } - - groupAIndex++; - processed++; - } - - if (groupAIndex >= groupADivisors.Length) - { - break; - } - - if ((ulong)groupADivisors[groupAIndex] > sqrtLimit) - { - groupAIndex = groupADivisors.Length; - break; - } - } - - return EvaluateWithOpenNumericFallback(n); - } - - ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); - if (endingOrder.IsEmpty) - { - return EvaluateWithOpenNumericFallback(n); - } - - Span indices = endingOrder.Length <= 8 - ? stackalloc int[endingOrder.Length] - : new int[endingOrder.Length]; - - bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; - bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); - - while (groupAHasMore || groupBHasMore) - { - if (groupAHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) - { - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (divisor > sqrtLimit) - { - groupAIndex = groupADivisors.Length; - groupAHasMore = false; - break; - } - - if (n % divisor == 0UL) - { - return false; - } - - groupAIndex++; - processed++; - } - - if (groupAIndex >= groupADivisors.Length) - { - groupAHasMore = false; - } - else - { - groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; - if (!groupAHasMore) - { - groupAIndex = groupADivisors.Length; - } - } - } - - if (groupBHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize) - { - if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) - { - groupBHasMore = false; - break; - } - - if (n % divisor == 0UL) - { - return false; - } - - processed++; - } - - if (groupBHasMore) - { - groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); - } - } - } - - return true; - } - - - - private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) - { - int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); - int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); - var divisorPool = ThreadStaticPools.UlongPool; - var hitPool = ThreadStaticPools.BytePool; - - ulong[]? divisorArray = null; - byte[]? hitFlags = null; - - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = PrimeTesterGpuContextPool.Rent(); - - bool compositeDetected = false; - - try - { - divisorArray = divisorPool.Rent(batchCapacity); - hitFlags = hitPool.Rent(batchCapacity); - - var accelerator = gpu.Accelerator; - var state = GpuKernelState.GetOrCreate(accelerator); - - lock (gpu.ExecutionLock) - { - var scratch = state.RentScratch(batchCapacity, accelerator); - try - { - int count = 0; - - bool ProcessBatch(int length) - { - if (length <= 0) - { - return false; - } - - scratch.Input.View.CopyFromCPU(ref divisorArray![0], length); - state.HeuristicTrialDivisionKernel(length, scratch.Input.View, n, scratch.Output.View); - accelerator.Synchronize(); - scratch.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; - } - - bool QueueDivisor(ulong divisor) - { - divisorArray![count] = divisor; - count++; - - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - return true; - } - - count = 0; - } - - 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]; - - bool groupAHasMore = groupAIndex < groupADivisors.Length && (ulong)groupADivisors[groupAIndex] <= sqrtLimit; - bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, sqrtLimit); - - while (!compositeDetected && (groupAHasMore || groupBHasMore)) - { - if (groupAHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) - { - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (divisor > sqrtLimit) - { - groupAIndex = groupADivisors.Length; - groupAHasMore = false; - break; - } - - compositeDetected = QueueDivisor(divisor); - if (compositeDetected) - { - break; - } - - groupAIndex++; - processed++; - } - - if (compositeDetected) - { - break; - } - - if (groupAIndex >= groupADivisors.Length) - { - groupAHasMore = false; - } - else - { - groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; - if (!groupAHasMore) - { - groupAIndex = groupADivisors.Length; - } - } - } - - if (groupBHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize) - { - if (!TrySelectNextGroupBDivisor(endingOrder, indices, sqrtLimit, out ulong divisor)) - { - groupBHasMore = false; - break; - } - - compositeDetected = QueueDivisor(divisor); - if (compositeDetected) - { - break; - } - - processed++; - } - - if (compositeDetected) - { - break; - } - - if (groupBHasMore) - { - groupBHasMore = HasGroupBCandidates(endingOrder, indices, sqrtLimit); - } - } - } - - if (!compositeDetected && count > 0) - { - compositeDetected = ProcessBatch(count); - } - } - finally - { - state.ReturnScratch(scratch); - } - } - } - finally - { - if (divisorArray is not null) - { - divisorPool.Return(divisorArray); - } - - if (hitFlags is not null) - { - 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 GetGroupBEndingOrder(byte nMod10) => nMod10 switch - { - 1 => GroupBEndingOrderMod1, - 3 => GroupBEndingOrderMod3, - 7 => GroupBEndingOrderMod7, - 9 => GroupBEndingOrderMod9, - _ => ReadOnlySpan.Empty, - }; - - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool TrySelectNextGroupBDivisor(ReadOnlySpan endingOrder, Span indices, ulong sqrtLimit, out ulong divisor) - { - ulong bestCandidate = ulong.MaxValue; - int bestEndingIndex = -1; - - for (int i = 0; i < endingOrder.Length; i++) - { - ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); - int index = indices[i]; - if ((uint)index >= (uint)divisors.Length) - { - continue; - } - - ulong candidate = (ulong)divisors[index]; - if (candidate > sqrtLimit) - { - indices[i] = divisors.Length; - continue; - } - - if (bestEndingIndex == -1 || candidate < bestCandidate) - { - bestCandidate = candidate; - bestEndingIndex = i; - } - } - - if (bestEndingIndex == -1) - { - divisor = 0UL; - return false; - } - - divisor = bestCandidate; - indices[bestEndingIndex]++; - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool HasGroupBCandidates(ReadOnlySpan endingOrder, Span indices, ulong sqrtLimit) - { - for (int i = 0; i < endingOrder.Length; i++) - { - ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); - int index = indices[i]; - if ((uint)index >= (uint)divisors.Length) - { - continue; - } - - ulong candidate = (ulong)divisors[index]; - if (candidate > sqrtLimit) - { - indices[i] = divisors.Length; - continue; - } - - return true; - } - - return false; - } - - private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch - { - 1 => HeuristicPrimeSieves.GroupBDivisorsEnding1, - 3 => HeuristicPrimeSieves.GroupBDivisorsEnding3, - 7 => HeuristicPrimeSieves.GroupBDivisorsEnding7, - 9 => HeuristicPrimeSieves.GroupBDivisorsEnding9, - _ => ReadOnlySpan.Empty, - }; - - private static ReadOnlySpan GetWheel210ResiduesForEnding(byte ending) => ending switch - { - 1 => Wheel210ResiduesEnding1, - 7 => Wheel210ResiduesEnding7, - 9 => Wheel210ResiduesEnding9, - _ => ReadOnlySpan.Empty, - }; - - private static int InitializeGroupBStates(byte nMod10, Span buffer) - { - ReadOnlySpan endings = GetGroupBEndingOrder(nMod10); - int count = 0; - - for (int priorityIndex = 0; priorityIndex < endings.Length; priorityIndex++) - { - ReadOnlySpan residues = GetWheel210ResiduesForEnding(endings[priorityIndex]); - for (int residueIndex = 0; residueIndex < residues.Length; residueIndex++) - { - buffer[count++] = new HeuristicGroupBSequenceState(residues[residueIndex], endings[priorityIndex], (byte)priorityIndex); - } - } - - return count; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) - { - return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); - } - - [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 int groupAConstantIndex; - private ulong groupAWheelCandidate; - private int groupAIncrementIndex; - private bool groupAWheelActive; - private Span groupBStates; - - public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) - { - 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++) - { - 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; - } - - 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; - - ulong next = value + increment; - if (next <= value) - { - groupAWheelActive = false; - } - else - { - groupAWheelCandidate = next; - } - - 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; - 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; - - public void TryNormalize() - { - ulong candidate = Candidate; - if (candidate == ulong.MaxValue) - { - return; - } - - while (candidate <= 13UL) - { - ulong next = candidate + Wheel210; - if (next <= candidate) - { - Candidate = ulong.MaxValue; - return; - } - - candidate = next; - } - - Candidate = candidate; - } - - public void Advance() - { - ulong candidate = Candidate; - if (candidate == ulong.MaxValue) - { - return; - } - - ulong next = candidate + Wheel210; - Candidate = next > candidate ? next : ulong.MaxValue; - } - - - } - - 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); - } - } - - - internal static bool IsPrimeInternal(ulong n, CancellationToken ct) - { - if (!EnableHeuristicPrimeTesting) - { - return LegacyIsPrimeInternal(n, ct); - } - - return HeuristicIsPrimeCpu(n); - } - - public static bool LegacyIsPrimeInternal(ulong n, CancellationToken ct) - { - // EvenPerfectBitScanner streams monotonically increasing odd exponents that start at 136,279,841 and already exclude - // multiples of five. Factorization helpers (PrimeOrderCalculator, Pollard routines, and divisor-cycle warmups) reuse - // this path for arbitrary residues and cofactors, so keep the boolean guards even though the scanner never trips them. bool isTwo = n == 2UL; bool isOdd = (n & 1UL) != 0UL; - // TODO: Replace this modulo check with ULongExtensions.Mod5 so the CPU hot path reuses the benchmarked helper instead of `%`. ulong mod5 = n % 5UL; bool divisibleByFive = n > 5UL && mod5 == 0UL; @@ -1030,9 +90,6 @@ public sealed class PrimeTester(bool useInternal = false) break; } - // TODO: Route this small-prime filtering through the shared divisor-cycle cache once - // PrimeTester can consult it directly; the cached cycles avoid the repeated `%` work - // that slows these hot loops when sieving hundreds of millions of candidates. if (n % smallPrimeDivisors[i] == 0) { result = false; @@ -1045,27 +102,120 @@ public sealed class PrimeTester(bool useInternal = false) return result; } - public static bool EnableHeuristicPrimeTesting { get; set; } = true; + internal readonly struct HeuristicDivisorCandidate + { + private readonly HeuristicPrimeTester.HeuristicDivisorCandidate _inner; + + internal HeuristicDivisorCandidate(HeuristicPrimeTester.HeuristicDivisorCandidate inner) + { + _inner = inner; + } + + public ulong Value => _inner.Value; + + public HeuristicPrimeTester.HeuristicDivisorGroup Group => _inner.Group; + + public byte Ending => _inner.Ending; + + public byte PriorityIndex => _inner.PriorityIndex; + + public ushort WheelResidue => _inner.WheelResidue; + + internal HeuristicPrimeTester.HeuristicDivisorCandidate AsHeuristic() => _inner; + } + + internal readonly struct HeuristicDivisorPreparation + { + private readonly HeuristicPrimeTester.HeuristicDivisorPreparation _inner; + + internal HeuristicDivisorPreparation(HeuristicPrimeTester.HeuristicDivisorPreparation inner) + { + _inner = inner; + } + + public HeuristicDivisorCandidate Candidate => new(_inner.Candidate); + + public MontgomeryDivisorData DivisorData => _inner.DivisorData; + + public ulong CycleLengthHint => _inner.CycleLengthHint; + + public bool HasCycleLengthHint => _inner.HasCycleLengthHint; + + public bool RequiresCycleComputation => _inner.RequiresCycleComputation; + + internal HeuristicPrimeTester.HeuristicDivisorPreparation AsHeuristic() => _inner; + } + + internal ref struct MersenneHeuristicDivisorEnumerator + { + private HeuristicPrimeTester.MersenneHeuristicDivisorEnumerator _inner; + + public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + _inner = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(exponent, maxDivisor); + } + + public ulong ProcessedCount => _inner.ProcessedCount; + + public ulong LastDivisor => _inner.LastDivisor; + + public bool Exhausted => _inner.Exhausted; + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + if (_inner.TryGetNext(out var innerCandidate)) + { + candidate = new HeuristicDivisorCandidate(innerCandidate); + return true; + } + + candidate = default; + return false; + } + } + + [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) + { + var innerCandidate = candidate.AsHeuristic(); + var preparation = HeuristicPrimeTester.PrepareHeuristicDivisor(in innerCandidate); + return new HeuristicDivisorPreparation(preparation); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong ResolveHeuristicCycleLength( + ulong exponent, + in HeuristicDivisorPreparation preparation, + out bool cycleFromHint, + out bool cycleComputed, + out bool primeOrderFailed) + { + var innerPreparation = preparation.AsHeuristic(); + return HeuristicPrimeTester.ResolveHeuristicCycleLength( + exponent, + in innerPreparation, + out cycleFromHint, + out cycleComputed, + out primeOrderFailed); + } public static int GpuBatchSize { get; set; } = 262_144; - public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096; - - public static int HeuristicDivisorInterleaveBatchSize { get; set; } = 64; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void IsPrimeBatchGpu(ReadOnlySpan values, Span results) { - // Limit concurrency and declare variables outside loops for performance and reuse var limiter = GpuPrimeWorkLimiter.Acquire(); var gpu = PrimeTesterGpuContextPool.Rent(); try { var accelerator = gpu.Accelerator; - - // Get per-accelerator cached kernel and device primes buffer. var state = GpuKernelState.GetOrCreate(accelerator); int totalLength = values.Length; int batchSize = Math.Max(1, GpuBatchSize); @@ -1076,9 +226,6 @@ public sealed class PrimeTester(bool useInternal = false) var input = scratch.Input; var output = scratch.Output; ulong[] temp = ArrayPool.Shared.Rent(batchSize); - // TODO: Replace this ad-hoc ArrayPool buffer with the pinned span cache from - // PrimeSieveGpuBenchmarks so batch uploads reuse preallocated GPU-friendly - // memory and avoid the extra copy before every kernel launch. try { @@ -1089,7 +236,7 @@ public sealed class PrimeTester(bool useInternal = false) int count = remaining > batchSize ? batchSize : remaining; values.Slice(pos, count).CopyTo(temp); - input.View.CopyFromCPU(ref temp[0], count); + input.View.CopyFromCPU(ref temp[0], count); state.Kernel(count, input.View, state.DevicePrimes.View, output.View); accelerator.Synchronize(); @@ -1112,7 +259,7 @@ public sealed class PrimeTester(bool useInternal = false) } } - private static class PrimeTesterGpuContextPool + internal static class PrimeTesterGpuContextPool { internal sealed class PooledContext : IDisposable { @@ -1199,7 +346,7 @@ public sealed class PrimeTester(bool useInternal = false) } // Per-accelerator GPU state for prime sieve (kernel + uploaded primes). - private sealed class KernelState + internal sealed class KernelState { public Action, ArrayView, ArrayView> Kernel { get; } public Action, ulong, ArrayView> HeuristicTrialDivisionKernel { get; } @@ -1276,7 +423,7 @@ public sealed class PrimeTester(bool useInternal = false) } } - private static class GpuKernelState + internal static class GpuKernelState { // Map accelerator to cached state; use Lazy to serialize kernel creation private static readonly System.Collections.Concurrent.ConcurrentDictionary> States = new(); diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md index 96da69e..46868da 100644 --- a/docs/HeuristicIsPrimePlan.md +++ b/docs/HeuristicIsPrimePlan.md @@ -109,7 +109,7 @@ internal static bool IsPrimeInternal(ulong n, CancellationToken ct) } ``` -## 3. Implementation Plan for `PrimeTester.HeuristicIsPrime*` +## 3. Implementation Plan for `HeuristicPrimeTester` The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, residue priorities) and connects them to both CPU and GPU execution, including every by-divisor integration point. Items already completed are marked with `[done]`. @@ -117,12 +117,12 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, 2. [done] **Introduce placeholder heuristic entry points** – add `HeuristicIsPrimeCpu` and `HeuristicIsPrimeGpu` by copying the current internal and GPU methods so subsequent commits can iterate without disrupting legacy behavior. 3. [done] **Design shared heuristic scaffolding** * [done] Added `HeuristicPrimeSieves` to precompute 4M-entry Group A/B divisor tables so heuristic prime checks reuse cached sequences without regenerating wheel steps at runtime; the enumerator remains available for residue-driven consumers. - * [done] Surface helpers that accept the precomputed `R = ⌊√n⌋` and `nMod10` values from callers via the new overload of `HeuristicIsPrimeCpu`; PrimeTester still computes them locally until call sites forward the parameters. + * [done] Surface helpers that accept the precomputed `R = ⌊√n⌋` and `nMod10` values from callers via the new overload of `HeuristicIsPrimeCpu`; HeuristicPrimeTester still computes them locally until call sites forward the parameters. * [done] Simplified the heuristic execution path so it runs without cancellation tokens or callback plumbing while keeping the shared enumerator available to other components. 4. **Share residue and wheel state with by-divisor sessions** * [done] Extracted the rolling-modulus helpers (mod 10/8/5/3/7/11) and decimal masks from `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` into the shared `MersenneDivisorResidueStepper` so heuristic prime tests and by-divisor sessions consume the same filters. * [done] Audited how `MersenneCpuDivisorScanSession` advances residues via `ExponentRemainderStepperCpu`/`CycleRemainderStepper` and updated the GPU by-divisor path to reuse `MersenneDivisorResidueStepper`, keeping CPU/GPU residue stepping aligned. - * [done] Publish wheel iterators that can stream candidates directly into `DivisorCycleCache`, GPU hit buffers, or CPU Montgomery checks without recomputing offsets by exposing `PrimeTester.CreateHeuristicDivisorEnumerator` and the reusable `HeuristicGroupBSequenceState` buffer. + * [done] Publish wheel iterators that can stream candidates directly into `DivisorCycleCache`, GPU hit buffers, or CPU Montgomery checks without recomputing offsets by exposing `HeuristicPrimeTester.CreateHeuristicDivisorEnumerator` and the reusable `HeuristicGroupBSequenceState` buffer. 5. **Expose cycle-preparation helpers** * [done] Extend the heuristic scaffolding so it returns `MontgomeryDivisorData` descriptors and cycle-length hints for each candidate divisor via `HeuristicDivisorPreparation` and `PrepareHeuristicDivisor`. * Wire these outputs into the existing cycle calculators to avoid duplicate cache lookups inside the CPU/GPU divisor scans. @@ -136,16 +136,16 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, * [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. 8. [done] **Refactor by-divisor CPU/GPU testers to reuse heuristics** - * [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `PrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop while preserving the existing status counters. + * [done] Drove `MersenneNumberDivisorByDivisorCpuTester.CheckDivisors` from `HeuristicPrimeTester.CreateMersenneDivisorEnumerator`, eliminating the bespoke residue loop while preserving the existing status counters. * [done] Streamed GPU batches from the same enumerator and Montgomery preparation data, removing the filtered-divisor scratch arrays while continuing to feed the existing kernel and hit bookkeeping. 9. [done] **Route legacy prime checks through heuristics** - * [done] Update `IsPrimeInternal` and `IsPrimeGpu` to delegate to the heuristic implementations when `EnableHeuristicPrimeTesting` is enabled, while allowing `PrimeTester(useInternal: true)` to force the legacy CPU/GPU paths. - * [done] Maintain legacy paths for benchmarking via `LegacyIsPrimeInternal`, `LegacyIsPrimeGpu`, and the global `EnableHeuristicPrimeTesting` switch so benchmark harnesses can compare behaviours. + * [done] Introduce `HeuristicPrimeTester` so heuristic CPU/GPU checks live alongside the legacy `PrimeTester` implementations while callers opt in explicitly. + * [done] Maintain legacy paths for benchmarking via `PrimeTester.IsPrimeInternal`, `PrimeTester.IsPrimeGpuFallback`, and the global `HeuristicPrimeTester.EnableHeuristicPrimeTesting` switch so benchmark harnesses can compare behaviours. 10. **Validation and benchmarking** * Create targeted unit tests covering boundary cases for Group A/B transitions, wheel enumerations, residue deltas, and the stop-on-hit behaviour. * Expand GPU test harnesses to assert equivalence between CPU/GPU heuristics on sampled composite and prime inputs and to verify the residue-stream contract with by-divisor sessions. * Benchmark against historical data (hundreds of thousands of composites) to confirm that heuristic ordering improves early divisor detection and that residue sharing reduces duplicate computation. - * [done] Build BenchmarkDotNet comparisons for `PrimeTester.IsPrimeGpu`, `PrimeTester.IsPrimeInternal`, and `Open.Numeric.Primes.Prime.Numbers.IsPrime` across representative ranges (≤100, ≤4,000,000, and ≥138,000,000) before routing heuristic code into production. Implemented in `EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs` to exercise both CPU and GPU heuristic paths alongside the external library. + * [done] Build BenchmarkDotNet comparisons for `HeuristicPrimeTester.IsPrimeGpu`, `HeuristicPrimeTester.HeuristicIsPrimeCpu`, and `Open.Numeric.Primes.Prime.Numbers.IsPrime` across representative ranges (≤100, ≤4,000,000, and ≥138,000,000) before routing heuristic code into production. Implemented in `EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs` to exercise both CPU and GPU heuristic paths alongside the external library. 11. **Documentation and rollout** * Update README/docs to explain heuristic prerequisites (odd, not divisible by 5) and how configuration toggles affect execution across prime testing and by-divisor scans. * Provide migration guidance for operators toggling between legacy and heuristic modes and describe how telemetry/metrics change after integration. diff --git a/docs/HeuristicIsPrimeSpecification.md b/docs/HeuristicIsPrimeSpecification.md index 30de040..c0282a7 100644 --- a/docs/HeuristicIsPrimeSpecification.md +++ b/docs/HeuristicIsPrimeSpecification.md @@ -1,6 +1,6 @@ # Heuristic Prime Tester Specification -This document describes the heuristic prime testing procedure implemented by `PrimeTester.HeuristicIsPrimeCpu` and `PrimeTester.HeuristicIsPrimeGpu`. The guidance consolidates the residue analysis and empirical statistics for odd inputs that are not divisible by five, based on the first 500,000 odd composite numbers after removing multiples of five. It also records the integration contracts required by the by-divisor pipeline so the heuristic covers every path that participates in divisor searches, residue sieves, and factorisation checks. +This document describes the heuristic prime testing procedure implemented by `HeuristicPrimeTester.HeuristicIsPrimeCpu` and `HeuristicPrimeTester.IsPrimeGpu`. The guidance consolidates the residue analysis and empirical statistics for odd inputs that are not divisible by five, based on the first 500,000 odd composite numbers after removing multiples of five. It also records the integration contracts required by the by-divisor pipeline so the heuristic covers every path that participates in divisor searches, residue sieves, and factorisation checks. ## Input Assumptions @@ -61,9 +61,9 @@ No other endings are allowed or reordered. ## Default Prime Tester Entry Points -* `PrimeTester.IsPrime` and `PrimeTester.IsPrimeGpu` now delegate to the heuristic CPU/GPU implementations by default. -* Set `PrimeTester.EnableHeuristicPrimeTesting = false` or instantiate `new PrimeTester(useInternal: true)` to force the legacy trial-division and GPU-sieve logic without removing the heuristic helpers. -* Benchmark harnesses can call `PrimeTester.LegacyIsPrimeInternal` (and `PrimeTester(useInternal: true).IsPrimeGpu`) to compare the legacy behaviour with the heuristic flow when collecting performance baselines. +* `PrimeTester.IsPrime` and `PrimeTester.IsPrimeGpu` continue to run the legacy trial-division and GPU sieve paths. Use `HeuristicPrimeTester` when the heuristic CPU/GPU flow is required. +* Toggle `HeuristicPrimeTester.EnableHeuristicPrimeTesting` or fall back to `PrimeTester.IsPrimeInternal` / `PrimeTester.IsPrimeGpuFallback` to control heuristic rollouts in test harnesses. +* Benchmark harnesses can compare behaviours by calling the legacy helpers on `PrimeTester` alongside the heuristic equivalents on `HeuristicPrimeTester`. ## Integration with the By-Divisor Pipeline @@ -72,9 +72,9 @@ The by-divisor executors (`MersenneNumberDivisorByDivisorCpuTester` and `Mersenn * **Residue filters** – The CPU scan maintains rolling residues modulo 10/8/5/3/7/11 to match the filters in `MersenneCpuDivisorScanSession.CheckDivisors`. `HeuristicIsPrime*` must either share the residue update helpers or produce identical stepping deltas so the divisor scanners can reuse the pre-filtered candidate stream without recomputing moduli. The shared `MersenneDivisorResidueStepper` now centralises those modulo calculations so both paths reuse the same state transitions, and the GPU by-divisor tester instantiates the same stepper before filtering candidates. * **Stepper integration** – CPU sessions advance residues and unity checks through `ExponentRemainderStepperCpu` together with `CycleRemainderStepper`; GPU kernels mirror this contract via `ExponentRemainderStepperGpu`. The heuristic implementation must call into these steppers (or emit the same remainder deltas) so the shared stepping state remains aligned across CPU/GPU scans and factor validation. * **Wheel alignment** – `HeuristicPrimeSieves` precomputes the wheel 30/210 sequences once at startup so both the heuristic prime tester and the by-divisor executors consume identical Group A/B streams without regenerating residues on every call. -* **Divisor streaming API** – `PrimeTester.CreateHeuristicDivisorEnumerator` exposes the same precomputed sieves so by-divisor scans can pull candidates directly into Montgomery preparation and GPU kernels without bespoke residue loops. -* **GPU batching** – `HeuristicIsPrimeGpu` streams Group A/B candidates from `HeuristicPrimeSieves`, packs them into batches controlled by `PrimeTester.HeuristicGpuDivisorBatchSize`, and launches `PrimeTesterKernels.HeuristicTrialDivisionKernel` to test each divisor with standard `%` arithmetic. When GPUs are disabled, the method delegates to the CPU heuristic so fallbacks reuse the same candidate ordering without additional plumbing. -* **Cycle preparation** – When the heuristic discovers a divisor candidate, it should populate `MontgomeryDivisorData` and expose the computed cycle length (or the need to recompute it) so that `DivisorCycleCache` and `MersenneDivisorCycles` do not duplicate work for the same `(prime, divisor)` pair during factor verification. `PrimeTester.PrepareHeuristicDivisor` together with `PrimeTester.ResolveHeuristicCycleLength` now provides this payload to both CPU and GPU by-divisor scans, and the trial-division summaries capture the data whenever a divisor is found. +* **Divisor streaming API** – `HeuristicPrimeTester.CreateHeuristicDivisorEnumerator` exposes the same precomputed sieves so by-divisor scans can pull candidates directly into Montgomery preparation and GPU kernels without bespoke residue loops. +* **GPU batching** – `HeuristicIsPrimeGpu` streams Group A/B candidates from `HeuristicPrimeSieves`, packs them into batches controlled by `HeuristicPrimeTester.HeuristicGpuDivisorBatchSize`, and launches `PrimeTesterKernels.HeuristicTrialDivisionKernel` to test each divisor with standard `%` arithmetic. When GPUs are disabled, the method delegates to the CPU heuristic so fallbacks reuse the same candidate ordering without additional plumbing. +* **Cycle preparation** – When the heuristic discovers a divisor candidate, it should populate `MontgomeryDivisorData` and expose the computed cycle length (or the need to recompute it) so that `DivisorCycleCache` and `MersenneDivisorCycles` do not duplicate work for the same `(prime, divisor)` pair during factor verification. `HeuristicPrimeTester.PrepareHeuristicDivisor` together with `HeuristicPrimeTester.ResolveHeuristicCycleLength` now provides this payload to both CPU and GPU by-divisor scans, and the trial-division summaries capture the data whenever a divisor is found. * **Caching contract** – The by-divisor GPU implementation pools `DivisorScanSession` instances and reuses GPU buffers. The heuristic must accept externally allocated spans/buffers for divisor candidates, cycle lengths, and hits so the existing resource pools remain valid. @@ -103,5 +103,5 @@ The following summaries cover the first 500,000 odd composite numbers after remo ## Benchmark Reporting -Benchmark results comparing `PrimeTester.IsPrimeGpu`, `PrimeTester.IsPrimeInternal`, and `Open.Numeric.Primes.Prime.Numbers.IsPrime` will be documented here once the dedicated benchmarks covering ≤100, ≤4,000,000, and ≥138,000,000 candidates are executed. The specification will be updated with the measured timings and analysis so operators can evaluate heuristic rollouts against legacy and external implementations. Benchmark harnesses are available via `EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs`, which runs the heuristic CPU and GPU paths alongside `Open.Numeric.Primes` for the three target ranges. +Benchmark results comparing `HeuristicPrimeTester.IsPrimeGpu`, `HeuristicPrimeTester.HeuristicIsPrimeCpu`, and `Open.Numeric.Primes.Prime.Numbers.IsPrime` will be documented here once the dedicated benchmarks covering ≤100, ≤4,000,000, and ≥138,000,000 candidates are executed. The specification will be updated with the measured timings and analysis so operators can evaluate heuristic rollouts against legacy and external implementations. Benchmark harnesses are available via `EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs`, which runs the heuristic CPU and GPU paths alongside `Open.Numeric.Primes` for the three target ranges.