From 43c7bafd5d1c7bcf3bd53daa6dea11cb5f863bd0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 15:19:49 +0100 Subject: [PATCH 1/9] Restore CPU PrimeTester usage and rename fallback helpers --- .../PrimeTesterBenchmarks.cs | 19 +- ...MersenneNumberDivisorByDivisorCpuTester.cs | 4 + ...MersenneNumberDivisorByDivisorGpuTester.cs | 4 + PerfectNumbers.Core/HeuristicPrimeTester.cs | 968 +++++++++++++++ PerfectNumbers.Core/MersenneDivisorCycles.cs | 2 + .../MersennePrimeFactorTester.cs | 2 + .../PrimeOrderCalculator.Cpu.cs | 12 +- .../PrimeOrderCalculator.Shared.cs | 1 + PerfectNumbers.Core/PrimeTester.cs | 1091 ++--------------- docs/HeuristicIsPrimePlan.md | 14 +- docs/HeuristicIsPrimeSpecification.md | 16 +- 11 files changed, 1130 insertions(+), 1003 deletions(-) create mode 100644 PerfectNumbers.Core/HeuristicPrimeTester.cs 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. From a83ae56caf3fc4dbbf6ca8c10b3069d11d3200cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 16:04:05 +0100 Subject: [PATCH 2/9] Use LastDigit enum across Mersenne scanners --- ...MersenneNumberIncrementalCpuTesterTests.cs | 5 ++- ...MersenneNumberIncrementalGpuTesterTests.cs | 5 ++- .../MersenneNumberOrderCpuTesterTests.cs | 5 ++- .../MersenneNumberOrderGpuTesterTests.cs | 7 ++-- .../MersenneNumberResidueCpuTesterTests.cs | 5 ++- .../MersenneNumberResidueGpuTesterTests.cs | 5 ++- .../Cpu/MersenneNumberIncrementalCpuTester.cs | 2 +- .../Cpu/MersenneNumberOrderCpuTester.cs | 4 +- .../Cpu/MersenneNumberResidueCpuTester.cs | 8 ++-- .../Gpu/MersenneNumberIncrementalGpuTester.cs | 4 +- .../Gpu/MersenneNumberOrderGpuTester.cs | 4 +- .../Gpu/MersenneNumberResidueGpuTester.cs | 4 +- PerfectNumbers.Core/LastDigit.cs | 10 +++++ PerfectNumbers.Core/MersenneNumberTester.cs | 26 ++++++------- PerfectNumbers.Core/PrimesGenerator.cs | 38 ++++++++++++++----- 15 files changed, 83 insertions(+), 49 deletions(-) create mode 100644 PerfectNumbers.Core/LastDigit.cs diff --git a/PerfectNumbers.Core.Tests/MersenneNumberIncrementalCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberIncrementalCpuTesterTests.cs index 9edcee3..dd7307d 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberIncrementalCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberIncrementalCpuTesterTests.cs @@ -1,4 +1,5 @@ using FluentAssertions; +using PerfectNumbers.Core; using PerfectNumbers.Core.Cpu; using Xunit; using UInt128 = System.UInt128; @@ -7,7 +8,7 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberIncrementalCpuTesterTests { - private static bool LastDigitIsSeven(ulong exponent) => (exponent & 3UL) == 3UL; + private static LastDigit GetLastDigit(ulong exponent) => (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; [Theory] [InlineData(GpuKernelType.Incremental)] @@ -30,7 +31,7 @@ public class MersenneNumberIncrementalCpuTesterTests private static void RunCase(MersenneNumberIncrementalCpuTester tester, ulong exponent, ulong maxK, bool expectedPrime) { bool isPrime = true; - tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime); + tester.Scan(exponent, (UInt128)exponent << 1, GetLastDigit(exponent), (UInt128)maxK, ref isPrime); isPrime.Should().Be(expectedPrime); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberIncrementalGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberIncrementalGpuTesterTests.cs index e17d04f..b990ece 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberIncrementalGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberIncrementalGpuTesterTests.cs @@ -1,4 +1,5 @@ using FluentAssertions; +using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; using Xunit; using UInt128 = System.UInt128; @@ -7,7 +8,7 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberIncrementalGpuTesterTests { - private static bool LastDigitIsSeven(ulong exponent) => (exponent & 3UL) == 3UL; + private static LastDigit GetLastDigit(ulong exponent) => (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; [Theory(Skip = "Disabled until the GPU incremental cycle heuristics are ported.")] [InlineData(GpuKernelType.Pow2Mod, false)] @@ -32,7 +33,7 @@ public class MersenneNumberIncrementalGpuTesterTests private static void RunCase(MersenneNumberIncrementalGpuTester tester, ulong exponent, ulong maxK, bool expectedPrime) { bool isPrime = true; - tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime); + tester.Scan(exponent, (UInt128)exponent << 1, GetLastDigit(exponent), (UInt128)maxK, ref isPrime); isPrime.Should().Be(expectedPrime); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberOrderCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberOrderCpuTesterTests.cs index 3348ceb..2d2756d 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberOrderCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberOrderCpuTesterTests.cs @@ -1,4 +1,5 @@ using FluentAssertions; +using PerfectNumbers.Core; using PerfectNumbers.Core.Cpu; using Xunit; using UInt128 = System.UInt128; @@ -7,7 +8,7 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberOrderCpuTesterTests { - private static bool LastDigitIsSeven(ulong exponent) => (exponent & 3UL) == 3UL; + private static LastDigit GetLastDigit(ulong exponent) => (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; [Theory] [InlineData(GpuKernelType.Incremental)] @@ -30,7 +31,7 @@ public class MersenneNumberOrderCpuTesterTests private static void RunCase(MersenneNumberOrderCpuTester tester, ulong exponent, ulong maxK, bool expectedPrime) { bool isPrime = true; - tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime); + tester.Scan(exponent, (UInt128)exponent << 1, GetLastDigit(exponent), (UInt128)maxK, ref isPrime); isPrime.Should().Be(expectedPrime); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberOrderGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberOrderGpuTesterTests.cs index 5f9a69b..2063482 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberOrderGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberOrderGpuTesterTests.cs @@ -1,5 +1,6 @@ using System.Reflection; using FluentAssertions; +using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; using Xunit; using UInt128 = System.UInt128; @@ -8,7 +9,7 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberOrderGpuTesterTests { - private static bool LastDigitIsSeven(ulong exponent) => (exponent & 3UL) == 3UL; + private static LastDigit GetLastDigit(ulong exponent) => (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; [Theory] [InlineData(GpuKernelType.Pow2Mod, false)] @@ -64,7 +65,7 @@ public class MersenneNumberOrderGpuTesterTests tester.Scan( exponent, (UInt128)exponent << 1, - LastDigitIsSeven(exponent), + GetLastDigit(exponent), (UInt128)10UL, ref isPrime); @@ -86,7 +87,7 @@ public class MersenneNumberOrderGpuTesterTests private static void RunCase(MersenneNumberOrderGpuTester tester, ulong exponent, ulong maxK, bool expectedPrime) { bool isPrime = true; - tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime); + tester.Scan(exponent, (UInt128)exponent << 1, GetLastDigit(exponent), (UInt128)maxK, ref isPrime); isPrime.Should().Be(expectedPrime); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs index 6c43e0d..aad010d 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberResidueCpuTesterTests.cs @@ -1,4 +1,5 @@ using FluentAssertions; +using PerfectNumbers.Core; using PerfectNumbers.Core.Cpu; using Xunit; using UInt128 = System.UInt128; @@ -7,7 +8,7 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberResidueCpuTesterTests { - private static bool LastDigitIsSeven(ulong exponent) => (exponent & 3UL) == 3UL; + private static LastDigit GetLastDigit(ulong exponent) => (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; [Fact] [Trait("Category", "Fast")] @@ -24,7 +25,7 @@ public class MersenneNumberResidueCpuTesterTests private static void RunCase(MersenneNumberResidueCpuTester tester, ulong exponent, ulong maxK, bool expectedPrime) { bool isPrime = true; - tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime); + tester.Scan(exponent, (UInt128)exponent << 1, GetLastDigit(exponent), (UInt128)maxK, ref isPrime); isPrime.Should().Be(expectedPrime); } } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs index 13ffde9..6063de1 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberResidueGpuTesterTests.cs @@ -1,4 +1,5 @@ using FluentAssertions; +using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; using Xunit; using UInt128 = System.UInt128; @@ -7,7 +8,7 @@ namespace PerfectNumbers.Core.Tests; public class MersenneNumberResidueGpuTesterTests { - private static bool LastDigitIsSeven(ulong exponent) => (exponent & 3UL) == 3UL; + private static LastDigit GetLastDigit(ulong exponent) => (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; [Theory] [InlineData(false)] @@ -27,7 +28,7 @@ public class MersenneNumberResidueGpuTesterTests private static void RunCase(MersenneNumberResidueGpuTester tester, ulong exponent, ulong maxK, bool expectedPrime) { bool isPrime = true; - tester.Scan(exponent, (UInt128)exponent << 1, LastDigitIsSeven(exponent), (UInt128)maxK, ref isPrime); + tester.Scan(exponent, (UInt128)exponent << 1, GetLastDigit(exponent), (UInt128)maxK, ref isPrime); isPrime.Should().Be(expectedPrime); } } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs index 7506cc1..c794e07 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberIncrementalCpuTester.cs @@ -5,7 +5,7 @@ public class MersenneNumberIncrementalCpuTester(GpuKernelType kernelType) private static readonly MersenneDivisorCycles _cycles = MersenneDivisorCycles.Shared; private readonly GpuKernelType _kernelType = kernelType; - public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) + public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime) { UInt128 k = UInt128.One; UInt128 one = k; diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs index cc1996f..14c3638 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberOrderCpuTester.cs @@ -4,7 +4,7 @@ public class MersenneNumberOrderCpuTester(GpuKernelType kernelType) { private readonly GpuKernelType _kernelType = kernelType; - public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) + public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime) { UInt128 k = 1UL; var auto = new MersenneResidueAutomaton(exponent); @@ -17,7 +17,7 @@ public class MersenneNumberOrderCpuTester(GpuKernelType kernelType) // TODO: When this lookup misses the snapshot, invoke the configured device to compute the // single required cycle on demand without persisting it so the CPU order path retains the // cycle-stepping speedups while honoring the no-extra-cache constraint for large divisors. - int allowMask2 = lastIsSeven ? CpuConstants.LastSevenMask10 : CpuConstants.LastOneMask10; + int allowMask2 = lastDigit == LastDigit.Seven ? CpuConstants.LastSevenMask10 : CpuConstants.LastOneMask10; bool shouldCheck = ((allowMask2 >> (int)auto.Mod10R) & 1) != 0; if (shouldCheck) { diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs index 8ee835d..fd19c49 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberResidueCpuTester.cs @@ -5,7 +5,7 @@ public class MersenneNumberResidueCpuTester private ModResidueTracker? _mersenneResidueTracker; // CPU residue variant using tracker + unrolled residue updates (no method calls in the loop). - public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) + public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime) { // Initialize/update Mersenne residue tracker and start a merge walk for ascending q divisors @@ -29,9 +29,9 @@ public class MersenneNumberResidueCpuTester step10 = step10.Mod10(); // Allowed last-digit sets for q depending on last digit of M_p - // lastIsSeven == true => allow {7,9} - // lastIsSeven == false => allow {1,3,7,9} - int allowMask = lastIsSeven ? CpuConstants.LastSevenMask10 : CpuConstants.LastOneMask10; + // LastDigit.Seven => allow {7,9} + // otherwise => allow {1,3,7,9} + int allowMask = lastDigit == LastDigit.Seven ? CpuConstants.LastSevenMask10 : CpuConstants.LastOneMask10; // Predeclare temps to avoid redeclaration overhead in the loop UInt128 remaining; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs index 871c201..ccc1bd9 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberIncrementalGpuTester.cs @@ -12,7 +12,7 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u private readonly GpuKernelType _kernelType = kernelType; private readonly bool _useGpuOrder = useGpuOrder; - public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) + public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime) { throw new NotImplementedException($"GPU incremental scanning requires the device cycle heuristics implementation (kernel: {_kernelType}, GPU order: {_useGpuOrder})."); @@ -23,7 +23,7 @@ public class MersenneNumberIncrementalGpuTester(GpuKernelType kernelType, bool u int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy UInt128 kStart = 1UL; ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL; - byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters + byte last = lastDigit == LastDigit.Seven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters var pow2Kernel = gpuLease.Pow2ModKernel; var incKernel = gpuLease.IncrementalKernel; diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs index d52970f..a4d0a60 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberOrderGpuTester.cs @@ -7,7 +7,7 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO private readonly GpuKernelType _kernelType = kernelType; private readonly bool _useGpuOrder = useGpuOrder; - public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) + public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime) { var gpuLease = GpuKernelPool.GetKernel(_useGpuOrder); var execution = gpuLease.EnterExecutionScope(); @@ -16,7 +16,7 @@ public class MersenneNumberOrderGpuTester(GpuKernelType kernelType, bool useGpuO int batchSize = GpuConstants.ScanBatchSize; // large batch improves GPU occupancy and avoids TDR UInt128 kStart = UInt128.One; ulong divMul = (ulong)((((UInt128)1 << 64) - UInt128.One) / exponent) + 1UL; - byte last = lastIsSeven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters + byte last = lastDigit == LastDigit.Seven ? (byte)1 : (byte)0; // ILGPU kernels do not support bool parameters var kernel = _kernelType switch { diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index 91a6b48..dac5e5a 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -24,7 +24,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) // requested by the caller, skip queuing additional block generation, and keep the snapshot cache untouched. // GPU residue variant: check 2^p % q == 1 for q = 2*p*k + 1. - public void Scan(ulong exponent, UInt128 twoP, bool lastIsSeven, UInt128 maxK, ref bool isPrime) + public void Scan(ulong exponent, UInt128 twoP, LastDigit lastDigit, UInt128 maxK, ref bool isPrime) { var limiter = GpuPrimeWorkLimiter.Acquire(); var gpuLease = RentAccelerator(); @@ -46,7 +46,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) UInt128 kStart = 1UL; UInt128 limit = maxK + UInt128.One; - byte last = lastIsSeven ? (byte)1 : (byte)0; + byte last = lastDigit == LastDigit.Seven ? (byte)1 : (byte)0; twoP.Mod10_8_5_3(out ulong step10, out ulong step8, out ulong step5, out ulong step3); step10 = step10.Mod10(); diff --git a/PerfectNumbers.Core/LastDigit.cs b/PerfectNumbers.Core/LastDigit.cs new file mode 100644 index 0000000..c53093f --- /dev/null +++ b/PerfectNumbers.Core/LastDigit.cs @@ -0,0 +1,10 @@ +namespace PerfectNumbers.Core; + +public enum LastDigit : byte +{ + One, + Three, + Five, + Seven, + Nine, +} diff --git a/PerfectNumbers.Core/MersenneNumberTester.cs b/PerfectNumbers.Core/MersenneNumberTester.cs index 0a5f763..cce3324 100644 --- a/PerfectNumbers.Core/MersenneNumberTester.cs +++ b/PerfectNumbers.Core/MersenneNumberTester.cs @@ -58,7 +58,7 @@ public sealed class MersenneNumberTester( } UInt128 twoP = (UInt128)exponent << 1; // 2 * p - bool lastIsSeven = (exponent & 3UL) == 3UL; + LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; // CPU path: compute orders directly and cache if (!_useGpuOrder) { @@ -71,7 +71,7 @@ public sealed class MersenneNumberTester( remainder = q.Mod10(); // TODO: Replace these repeated Mod10/Mod8/Mod3/Mod5 calls with a residue automaton walk so // order warmups reuse the optimized cycle stepping validated in the residue benchmarks. - bool shouldCheck = lastIsSeven + bool shouldCheck = lastDigit == LastDigit.Seven ? remainder == 7UL || remainder == 9UL || remainder == 3UL : remainder == 1UL || remainder == 3UL || remainder == 9UL; @@ -112,7 +112,7 @@ public sealed class MersenneNumberTester( remainder2 = q2.Mod10(); // TODO: Reuse the same residue-automaton fast path here so the GPU warmup staging avoids `%` and // branches the benchmarks showed slower than the tracked stepping helper. - bool shouldCheck2 = lastIsSeven + bool shouldCheck2 = lastDigit == LastDigit.Seven ? remainder2 == 7UL || remainder2 == 9UL || remainder2 == 3UL : remainder2 == 1UL || remainder2 == 3UL || remainder2 == 9UL; @@ -230,7 +230,7 @@ public sealed class MersenneNumberTester( UInt128 twoP = (UInt128)exponent << 1; // 2 * p // UInt128 maxDivisor = new(ulong.MaxValue, ulong.MaxValue); UInt128 maxK = UInt128Numbers.Two.Pow(exponent) / 2 + 1; - bool lastIsSeven = (exponent & 3UL) == 3UL; + LastDigit lastDigit = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; // In residue mode the residue scan is the final primality check. if (_useResidue) @@ -242,11 +242,11 @@ public sealed class MersenneNumberTester( if (_useGpuScan) { - _residueGpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime); + _residueGpuTester!.Scan(exponent, twoP, lastDigit, maxK, ref prePrime); } else { - _residueCpuTester!.Scan(exponent, twoP, lastIsSeven, maxK, ref prePrime); + _residueCpuTester!.Scan(exponent, twoP, lastDigit, maxK, ref prePrime); } return prePrime; @@ -258,14 +258,14 @@ public sealed class MersenneNumberTester( // Scan a limited range of k to quickly reject obvious cases. UInt128 twoPPre = (UInt128)exponent << 1; UInt128 maxKPre = (UInt128)PreLucasPreScanK; - bool lastIsSevenPre = (exponent & 3UL) == 3UL; + LastDigit lastDigitPre = (exponent & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; if (_useGpuScan) { - _incrementalGpuTester.Scan(exponent, twoPPre, lastIsSevenPre, maxKPre, ref prePrime); + _incrementalGpuTester.Scan(exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime); } else { - _incrementalCpuTester.Scan(exponent, twoPPre, lastIsSevenPre, maxKPre, ref prePrime); + _incrementalCpuTester.Scan(exponent, twoPPre, lastDigitPre, maxKPre, ref prePrime); } if (!prePrime) @@ -284,22 +284,22 @@ public sealed class MersenneNumberTester( { if (_useGpuOrder) { - _orderGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); + _orderGpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); } else { - _orderCpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); + _orderCpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); } } else { if (_useGpuScan) { - _incrementalGpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); + _incrementalGpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); } else { - _incrementalCpuTester.Scan(exponent, twoP, lastIsSeven, maxK, ref isPrime); + _incrementalCpuTester.Scan(exponent, twoP, lastDigit, maxK, ref isPrime); } } diff --git a/PerfectNumbers.Core/PrimesGenerator.cs b/PerfectNumbers.Core/PrimesGenerator.cs index 754d302..da1f050 100644 --- a/PerfectNumbers.Core/PrimesGenerator.cs +++ b/PerfectNumbers.Core/PrimesGenerator.cs @@ -1,3 +1,4 @@ +using System; using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; @@ -85,14 +86,16 @@ public static class PrimesGenerator if (candidate != 2U) { - if (lastOneCount < targetInt && IsAllowedForLastOne(candidate)) + LastDigit lastDigit = GetLastDigit(candidate); + + if (lastOneCount < targetInt && IsAllowedForLastOne(candidate, lastDigit)) { lastOne[lastOneCount] = candidate; lastOnePow2[lastOneCount] = checked((ulong)candidate * candidate); lastOneCount++; } - if (lastSevenCount < targetInt && IsAllowedForLastSeven(candidate)) + if (lastSevenCount < targetInt && IsAllowedForLastSeven(candidate, lastDigit)) { lastSeven[lastSevenCount] = candidate; lastSevenPow2[lastSevenCount] = checked((ulong)candidate * candidate); @@ -108,25 +111,40 @@ public static class PrimesGenerator } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool IsAllowedForLastOne(uint prime) + private static LastDigit GetLastDigit(uint prime) { - // TODO: Swap the `% 10` usage for ULongExtensions.Mod10 so the hot classification path - // reuses the benchmarked residue helper instead of repeated divisions. return (prime % 10U) switch { - 1U or 3U or 9U => true, - _ => prime == 7U || prime == 11U, + 1U => LastDigit.One, + 3U => LastDigit.Three, + 5U => LastDigit.Five, + 7U => LastDigit.Seven, + 9U => LastDigit.Nine, + _ => throw new ArgumentOutOfRangeException(nameof(prime)), }; } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool IsAllowedForLastSeven(uint prime) + private static bool IsAllowedForLastOne(uint prime, LastDigit lastDigit) + { + // TODO: Swap the `% 10` usage for ULongExtensions.Mod10 so the hot classification path + // reuses the benchmarked residue helper instead of repeated divisions. + return lastDigit switch + { + LastDigit.One or LastDigit.Three or LastDigit.Nine => true, + LastDigit.Seven => prime == 7U, + _ => prime == 11U, + }; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool IsAllowedForLastSeven(uint prime, LastDigit lastDigit) { // TODO: Route this `% 10` classification through ULongExtensions.Mod10 to match the faster // residue helper used elsewhere in the scanner. - return (prime % 10U) switch + return lastDigit switch { - 3U or 7U or 9U => true, + LastDigit.Three or LastDigit.Seven or LastDigit.Nine => true, _ => prime == 11U, }; } From c01abbd03acc035e6fb361491eba9f2e976c218f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 16:04:09 +0100 Subject: [PATCH 3/9] Restore legacy by-divisor scans for Mersenne testers --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 191 ++++++++++++--- ...MersenneNumberDivisorByDivisorGpuTester.cs | 224 +++++++++++++++--- 2 files changed, 351 insertions(+), 64 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index c061da1..bfe2808 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -1,5 +1,6 @@ using System.Collections.Generic; using System.Runtime.CompilerServices; +using PerfectNumbers.Core; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core.Cpu; @@ -12,6 +13,16 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private ulong _lastStatusDivisor; private int _batchSize = 1_024; + [ThreadStatic] + private static GpuUInt128WorkSet _divisorScanGpuWorkSet; + + private struct GpuUInt128WorkSet + { + public GpuUInt128 Step; + public GpuUInt128 Divisor; + public GpuUInt128 Limit; + } + public int BatchSize { get => _batchSize; @@ -103,41 +114,152 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv out bool processedAll, out ulong processedCount) { - processedAll = false; processedCount = 0UL; + processedAll = false; - // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); - var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); + // The EvenPerfectBitScanner feeds primes >= 138,000,000 here, so allowedMax >= 3 in production runs. + // Keeping the guard commented out documents the reasoning for benchmarks and tests. + // if (allowedMax < 3UL) + // { + // return false; + // } - // while (enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate)) - while (enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate)) + ref GpuUInt128WorkSet workSet = ref _divisorScanGpuWorkSet; + + ref GpuUInt128 step = ref workSet.Step; + step.High = 0UL; + step.Low = prime; + step.ShiftLeft(1); + + ref GpuUInt128 limit = ref workSet.Limit; + limit.High = 0UL; + limit.Low = allowedMax; + + ref GpuUInt128 divisor = ref workSet.Divisor; + divisor.High = step.High; + divisor.Low = step.Low; + divisor.Add(1UL); + if (divisor.CompareTo(limit) > 0) { - 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, - out _, - out _, - out _); - - if (divisorCycle == prime) - { - processedAll = true; - return true; - } - - if (divisorCycle == 0UL) - { - Console.WriteLine($"Divisor cycle was not calculated for {prime}"); - } + processedAll = true; + divisor.High = 0UL; + divisor.Low = 0UL; + step.High = 0UL; + step.Low = 0UL; + limit.High = 0UL; + limit.Low = 0UL; + return false; } - processedCount = enumerator.ProcessedCount; - processedAll = enumerator.Exhausted; + // Intentionally recomputes factorizations without a per-thread cache. + // The previous factor cache recorded virtually no hits and only slowed down the scan. + DivisorCycleCache cycleCache = DivisorCycleCache.Shared; + + ulong stepHigh = step.High; + ulong stepLow = step.Low; + byte step10 = (byte)((((stepHigh % 10UL) * 6UL) + (stepLow % 10UL)) % 10UL); + byte step8 = (byte)(stepLow % 8UL); + byte step5 = (byte)(((stepHigh % 5UL) + (stepLow % 5UL)) % 5UL); + byte step3 = (byte)(((stepHigh % 3UL) + (stepLow % 3UL)) % 3UL); + byte step7 = (byte)((((stepHigh % 7UL) * 2UL) + (stepLow % 7UL)) % 7UL); + byte step11 = (byte)((((stepHigh % 11UL) * 5UL) + (stepLow % 11UL)) % 11UL); + + ulong divisorHigh = divisor.High; + ulong divisorLow = divisor.Low; + byte remainder10 = (byte)((((divisorHigh % 10UL) * 6UL) + (divisorLow % 10UL)) % 10UL); + byte remainder8 = (byte)(divisorLow % 8UL); + byte remainder5 = (byte)(((divisorHigh % 5UL) + (divisorLow % 5UL)) % 5UL); + byte remainder3 = (byte)(((divisorHigh % 3UL) + (divisorLow % 3UL)) % 3UL); + byte remainder7 = (byte)((((divisorHigh % 7UL) * 2UL) + (divisorLow % 7UL)) % 7UL); + byte remainder11 = (byte)((((divisorHigh % 11UL) * 5UL) + (divisorLow % 11UL)) % 11UL); + + // Keep the divisibility filters aligned with the divisor-cycle generator so the + // CPU path never requests cycles that were skipped during cache creation. + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + bool lastIsSeven = lastDigit == LastDigit.Seven; + const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); + const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + + while (divisor.CompareTo(limit) <= 0) + { + ulong candidate = divisor.Low; + processedCount++; + + ushort decimalMask = lastIsSeven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + bool admissible = ((decimalMask >> remainder10) & 1) != 0; + + if (admissible + && (remainder8 == 1 || remainder8 == 7) + && remainder3 != 0 + && remainder5 != 0 + && remainder7 != 0 + && remainder11 != 0) + { + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(candidate); + ulong divisorCycle; + // Divisors generated from 2 * k * p + 1 exceed the small-cycle snapshot when p >= 138,000,000, so the short path below never runs. + // if (candidate <= PerfectNumberConstants.MaxQForDivisorCycles) + // { + // divisorCycle = cycleCache.GetCycleLength(candidate); + // } + // else + { + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + candidate, + prime, + divisorData, + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) + { + // Divisors produced by 2 * k * p + 1 always exceed PerfectNumberConstants.MaxQForDivisorCycles + // for the exponents scanned here, so skip the unused cache fallback and compute directly. + divisorCycle = MersenneDivisorCycles.CalculateCycleLength( + candidate, + divisorData, + skipPrimeOrderHeuristic: primeOrderFailed); + } + else + { + divisorCycle = computedCycle; + } + } + + if (divisorCycle == prime) + { + // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides + // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. + divisor.High = 0UL; + divisor.Low = 0UL; + step.High = 0UL; + step.Low = 0UL; + limit.High = 0UL; + limit.Low = 0UL; + processedAll = true; + return true; + } + + if (divisorCycle == 0UL) + { + Console.WriteLine($"Divisor cycle was not calculated for {prime}"); + } + } + + divisor.Add(step); + remainder10 = AddMod(remainder10, step10, (byte)10); + remainder8 = AddMod(remainder8, step8, (byte)8); + remainder5 = AddMod(remainder5, step5, (byte)5); + remainder3 = AddMod(remainder3, step3, (byte)3); + remainder7 = AddMod(remainder7, step7, (byte)7); + remainder11 = AddMod(remainder11, step11, (byte)11); + } + + processedAll = divisor.CompareTo(limit) > 0; + divisor.High = 0UL; + divisor.Low = 0UL; + step.High = 0UL; + step.Low = 0UL; + limit.High = 0UL; + limit.Low = 0UL; return false; } @@ -147,6 +269,17 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return residue == 1UL ? (byte)1 : (byte)0; } + private static byte AddMod(byte value, byte delta, byte modulus) + { + int sum = value + delta; + if (sum >= modulus) + { + sum -= modulus; + } + + return (byte)sum; + } + private void UpdateStatusUnsafe(ulong processedCount) { if (processedCount == 0UL) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index e1bc1de..e1b1fe5 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -112,6 +112,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv resources.HitIndexBuffer, resources.Divisors, resources.Exponents, + resources.FilteredDivisors, resources.DivisorData, resources.Offsets, resources.Counts, @@ -176,6 +177,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv MemoryBuffer1D hitIndexBuffer, ulong[] divisors, ulong[] exponents, + ulong[] filteredDivisors, GpuDivisorPartialData[] divisorData, int[] offsets, int[] counts, @@ -186,40 +188,131 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { int batchCapacity = (int)divisorDataBuffer.Length; bool composite = false; + bool processedAll = false; processedCount = 0UL; lastProcessed = 0UL; - // var enumerator = HeuristicPrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); - var enumerator = PrimeTester.CreateMersenneDivisorEnumerator(prime, allowedMax); + int chunkCapacity = Math.Max(1, batchCapacity); - Span divisorSpan = divisors.AsSpan(); - Span exponentSpan = exponents.AsSpan(); - Span divisorDataSpan = divisorData.AsSpan(); - Span offsetSpan = offsets.AsSpan(); - Span countSpan = counts.AsSpan(); - Span cycleSpan = cycles.AsSpan(); - - while (!composite) + UInt128 twoP128 = (UInt128)prime << 1; + UInt128 allowedMax128 = allowedMax; + UInt128 firstDivisor128 = twoP128 + UInt128.One; + bool invalidStride = twoP128 == UInt128.Zero; + bool outOfRange = firstDivisor128 > allowedMax128; + UInt128 numerator = allowedMax128 - UInt128.One; + UInt128 maxK128 = invalidStride || outOfRange ? UInt128.Zero : numerator / twoP128; + bool hasCandidates = !invalidStride && !outOfRange && maxK128 != UInt128.Zero; + if (!hasCandidates) { + coveredRange = true; + return false; + } + + ulong maxK = maxK128 > ulong.MaxValue ? ulong.MaxValue : (ulong)maxK128; + ulong currentK = 1UL; + + byte step10 = (byte)(twoP128 % 10UL); + byte step8 = (byte)(twoP128 % 8UL); + byte step5 = (byte)(twoP128 % 5UL); + byte step3 = (byte)(twoP128 % 3UL); + byte step7 = (byte)(twoP128 % 7UL); + byte step11 = (byte)(twoP128 % 11UL); + + UInt128 currentDivisor128 = firstDivisor128; + byte remainder10 = (byte)((ulong)(currentDivisor128 % 10UL)); + byte remainder8 = (byte)((ulong)(currentDivisor128 & 7UL)); + byte remainder5 = (byte)((ulong)(currentDivisor128 % 5UL)); + byte remainder3 = (byte)((ulong)(currentDivisor128 % 3UL)); + byte remainder7 = (byte)((ulong)(currentDivisor128 % 7UL)); + byte remainder11 = (byte)((ulong)(currentDivisor128 % 11UL)); + LastDigit lastDigit = (prime & 3UL) == 3UL ? LastDigit.Seven : LastDigit.One; + + Span filteredStorage = filteredDivisors.AsSpan(); + Span divisorStorage = divisors.AsSpan(); + Span exponentStorage = exponents.AsSpan(); + Span divisorDataStorage = divisorData.AsSpan(); + Span offsetStorage = offsets.AsSpan(); + Span countStorage = counts.AsSpan(); + Span cycleStorage = cycles.AsSpan(); + while (currentK <= maxK && !composite) + { + int chunkCount = Math.Min(chunkCapacity, batchCapacity); + ulong remainingK = maxK - currentK + 1UL; + if ((ulong)chunkCount > remainingK) + { + chunkCount = (int)remainingK; + } + + if (chunkCount <= 0) + { + processedAll = true; + break; + } + + int filteredCount = 0; + UInt128 nextDivisor128 = currentDivisor128; + + byte localRemainder10 = remainder10; + byte localRemainder8 = remainder8; + byte localRemainder5 = remainder5; + byte localRemainder3 = remainder3; + byte localRemainder7 = remainder7; + byte localRemainder11 = remainder11; + + for (int i = 0; i < chunkCount; i++) + { + ulong candidate = (ulong)nextDivisor128; + + if (CandidatePassesHeuristics(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastDigit)) + { + filteredStorage[filteredCount++] = candidate; + } + + nextDivisor128 += twoP128; + localRemainder10 = AddMod(localRemainder10, step10, 10); + localRemainder8 = AddMod(localRemainder8, step8, 8); + localRemainder5 = AddMod(localRemainder5, step5, 5); + localRemainder3 = AddMod(localRemainder3, step3, 3); + localRemainder7 = AddMod(localRemainder7, step7, 7); + localRemainder11 = AddMod(localRemainder11, step11, 11); + } + + remainder10 = localRemainder10; + remainder8 = localRemainder8; + remainder5 = localRemainder5; + remainder3 = localRemainder3; + remainder7 = localRemainder7; + remainder11 = localRemainder11; + + processedCount += (ulong)chunkCount; + + UInt128 lastDivisor128 = nextDivisor128 - twoP128; + lastProcessed = (ulong)lastDivisor128; + currentDivisor128 = nextDivisor128; + + currentK += (ulong)chunkCount; + processedAll = currentK > maxK; + + if (filteredCount == 0) + { + continue; + } + + Span filteredDivisorsSpan = filteredStorage[..filteredCount]; + + Span divisorSpan = divisorStorage; + Span exponentSpan = exponentStorage; + Span divisorDataSpan = divisorDataStorage; + Span offsetSpan = offsetStorage; + Span countSpan = countStorage; + Span cycleSpan = cycleStorage; int admissibleCount = 0; int exponentIndex = 0; - - // while (admissibleCount < batchCapacity && enumerator.TryGetNext(out HeuristicPrimeTester.HeuristicDivisorCandidate candidate)) - while (admissibleCount < batchCapacity && enumerator.TryGetNext(out PrimeTester.HeuristicDivisorCandidate candidate)) + for (int i = 0; i < filteredCount; i++) { - 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, - out _, - out _, - out _); - + ulong divisorValue = filteredDivisorsSpan[i]; + MontgomeryDivisorData montgomeryData = MontgomeryDivisorData.FromModulus(divisorValue); + ulong divisorCycle = ResolveDivisorCycle(divisorValue, prime, in montgomeryData); if (divisorCycle != prime) { continue; @@ -236,16 +329,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv exponentIndex++; } - processedCount = enumerator.ProcessedCount; - lastProcessed = enumerator.LastDivisor; - if (admissibleCount == 0) { - if (enumerator.Exhausted) - { - break; - } - continue; } @@ -277,10 +362,74 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv lastProcessed = hitFound ? divisorSpan[hitIndex] : divisorSpan[lastIndex]; } - coveredRange = composite || enumerator.Exhausted; + coveredRange = composite || processedAll || (currentDivisor128 > allowedMax128); return composite; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong ResolveDivisorCycle(ulong divisor, ulong prime, in MontgomeryDivisorData divisorData) + { + if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu(divisor, prime, divisorData, out ulong computedCycle, out bool primeOrderFailed) || computedCycle == 0UL) + { + return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); + } + + return computedCycle; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool CandidatePassesHeuristics( + byte remainder10, + byte remainder8, + byte remainder3, + byte remainder5, + byte remainder7, + byte remainder11, + LastDigit lastDigit) + { + const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); + const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + + ushort decimalMask = lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + if (((decimalMask >> remainder10) & 1) == 0) + { + return false; + } + + if (remainder8 != 1 && remainder8 != 7) + { + return false; + } + + if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0) + { + return false; + } + + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte AddMod(byte value, byte delta, byte modulus) + { + int result = value + delta; + if (result >= modulus) + { + result -= modulus; + if (result >= modulus) + { + result -= modulus; + } + } + + return (byte)result; + } + + + + + private static ulong ComputeDivisorLimitFromMaxPrimeGpu(ulong maxPrime) { // The --mersenne=bydivisor flow in EvenPerfectBitScanner only calls ConfigureFromMaxPrime with primes greater than 1, @@ -571,6 +720,8 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv internal ulong[] Exponents { get; private set; } = null!; + internal ulong[] FilteredDivisors { get; private set; } = null!; + internal GpuDivisorPartialData[] DivisorData { get; private set; } = null!; internal int[] Offsets { get; private set; } = null!; @@ -645,6 +796,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv Divisors = ulongPool.Rent(capacity); Exponents = ulongPool.Rent(capacity); + FilteredDivisors = ulongPool.Rent(capacity); DivisorData = partialPool.Rent(capacity); Offsets = intPool.Rent(capacity); Counts = intPool.Rent(capacity); @@ -683,12 +835,14 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ArrayPool partialPool = ThreadStaticPools.GpuDivisorPartialDataPool; ulongPool.Return(Divisors, clearArray: false); ulongPool.Return(Exponents, clearArray: false); + ulongPool.Return(FilteredDivisors, clearArray: false); partialPool.Return(DivisorData, clearArray: false); intPool.Return(Offsets, clearArray: false); intPool.Return(Counts, clearArray: false); ulongPool.Return(Cycles, clearArray: false); Divisors = Array.Empty(); Exponents = Array.Empty(); + FilteredDivisors = Array.Empty(); DivisorData = Array.Empty(); Offsets = Array.Empty(); Counts = Array.Empty(); From 9bc8fd50433bcec2fec59c29e19d35b713387f34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 16:18:46 +0100 Subject: [PATCH 4/9] Remove heuristic hooks from PrimeTester --- ...MersenneNumberDivisorByDivisorCpuTester.cs | 18 ---- PerfectNumbers.Core/PrimeTester.cs | 102 ------------------ 2 files changed, 120 deletions(-) diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index bfe2808..5254d7a 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -142,12 +142,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv if (divisor.CompareTo(limit) > 0) { processedAll = true; - divisor.High = 0UL; - divisor.Low = 0UL; - step.High = 0UL; - step.Low = 0UL; - limit.High = 0UL; - limit.Low = 0UL; return false; } @@ -228,12 +222,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv { // A cycle equal to the tested exponent (which is prime in this path) guarantees that the candidate divides // the corresponding Mersenne number because the order of 2 modulo the divisor is exactly p. - divisor.High = 0UL; - divisor.Low = 0UL; - step.High = 0UL; - step.Low = 0UL; - limit.High = 0UL; - limit.Low = 0UL; processedAll = true; return true; } @@ -254,12 +242,6 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } processedAll = divisor.CompareTo(limit) > 0; - divisor.High = 0UL; - divisor.Low = 0UL; - step.High = 0UL; - step.Low = 0UL; - limit.High = 0UL; - limit.Low = 0UL; return false; } diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 9ec21b5..d1c8d2a 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -103,108 +103,6 @@ public sealed class PrimeTester } - 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; [MethodImpl(MethodImplOptions.AggressiveInlining)] From a7fc34e37abdb76b4dc908c51b2c7879528d3598 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 16:45:18 +0100 Subject: [PATCH 5/9] Inline PrimeTester CPU path and reuse thread pools --- PerfectNumbers.Core/PrimeTester.cs | 91 ++++++++++++++++-------------- 1 file changed, 50 insertions(+), 41 deletions(-) diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index d1c8d2a..4e56414 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -13,7 +13,7 @@ namespace PerfectNumbers.Core; public sealed class PrimeTester { - public PrimeTester(bool useInternal = false) + public PrimeTester() { } @@ -24,7 +24,47 @@ public sealed class PrimeTester public bool IsPrime(ulong n, CancellationToken ct) { - return IsPrimeInternal(n, ct); + // EvenPerfectBitScanner feeds monotonically increasing odd exponents that start well above the + // small-prime cutoffs and already exclude multiples of five. These guards remain for targeted + // tests and diagnostics so ad-hoc callers can still probe small or even values when needed. + bool isTwo = n == 2UL; + bool isOdd = (n & 1UL) != 0UL; + ulong mod5 = n % 5UL; + bool divisibleByFive = n > 5UL && mod5 == 0UL; + + bool result = n >= 2UL && (isTwo || isOdd) && !divisibleByFive; + + // Production scans always fall through to trial division because n is already > 5 by the time this + // path executes. Keep the guard so synthetic callers that hammer tiny inputs can skip the loop. + bool requiresTrialDivision = result && n >= 7UL && !isTwo; + + if (requiresTrialDivision) + { + bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); + result &= !sharesMaxExponentFactor; + + if (result) + { + var smallPrimeDivisorsLength = PrimesGenerator.SmallPrimes.Length; + uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; + ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; + for (int i = 0; i < smallPrimeDivisorsLength; i++) + { + if (smallPrimeDivisorsMul[i] > n) + { + break; + } + + if (n % smallPrimeDivisors[i] == 0) + { + result = false; + break; + } + } + } + } + + return result; } public static bool IsPrimeGpu(ulong n) @@ -65,41 +105,8 @@ public sealed class PrimeTester public static bool IsPrimeInternal(ulong n, CancellationToken ct) { - bool isTwo = n == 2UL; - bool isOdd = (n & 1UL) != 0UL; - ulong mod5 = n % 5UL; - bool divisibleByFive = n > 5UL && mod5 == 0UL; - - bool result = n >= 2UL && (isTwo || isOdd) && !divisibleByFive; - bool requiresTrialDivision = result && n >= 7UL && !isTwo; - - if (requiresTrialDivision) - { - bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); - result &= !sharesMaxExponentFactor; - - if (result) - { - var smallPrimeDivisorsLength = PrimesGenerator.SmallPrimes.Length; - uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; - ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; - for (int i = 0; i < smallPrimeDivisorsLength; i++) - { - if (smallPrimeDivisorsMul[i] > n) - { - break; - } - - if (n % smallPrimeDivisors[i] == 0) - { - result = false; - break; - } - } - } - } - - return result; + // Preserve the legacy entry point for callers that bypass the thread-local caches. + return Exclusive.IsPrime(n, ct); } @@ -123,7 +130,8 @@ public sealed class PrimeTester var scratch = state.RentScratch(batchSize, accelerator); var input = scratch.Input; var output = scratch.Output; - ulong[] temp = ArrayPool.Shared.Rent(batchSize); + ArrayPool pool = ThreadStaticPools.UlongPool; + ulong[] temp = pool.Rent(batchSize); try { @@ -145,7 +153,7 @@ public sealed class PrimeTester } finally { - ArrayPool.Shared.Return(temp); + pool.Return(temp, clearArray: false); state.ReturnScratch(scratch); } } @@ -385,13 +393,14 @@ public sealed class PrimeTester var inputBuffer = accelerator.Allocate1D(length); var resultBuffer = accelerator.Allocate1D(length); - ulong[] temp = ArrayPool.Shared.Rent(length); + ArrayPool pool = ThreadStaticPools.UlongPool; + ulong[] temp = pool.Rent(length); values.CopyTo(temp); inputBuffer.View.CopyFromCPU(ref temp[0], length); kernel(length, inputBuffer.View, resultBuffer.View); accelerator.Synchronize(); resultBuffer.View.CopyToCPU(ref results[0], length); - ArrayPool.Shared.Return(temp); + pool.Return(temp, clearArray: false); resultBuffer.Dispose(); inputBuffer.Dispose(); gpu.Dispose(); From 55554d9565724d87942dae13a0a615827494b8b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 17:06:13 +0100 Subject: [PATCH 6/9] Inline legacy GPU sieve and update heuristic fallbacks --- .../PrimeTesterBenchmarks.cs | 2 +- PerfectNumbers.Core/HeuristicPrimeTester.cs | 4 ++-- PerfectNumbers.Core/PrimeTester.cs | 14 +++++--------- docs/HeuristicIsPrimePlan.md | 2 +- docs/HeuristicIsPrimeSpecification.md | 2 +- 5 files changed, 10 insertions(+), 14 deletions(-) diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index 981de9b..d3bae87 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -86,7 +86,7 @@ public class PrimeTesterBenchmarks for (int i = 0; i < values.Length; i++) { - if (PrimeTester.IsPrimeGpuFallback(values[i], cancellationToken)) + if (PrimeTester.Exclusive.IsPrimeGpu(values[i], cancellationToken)) { primeCount++; } diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index d5d70a6..ca6e9ea 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -102,7 +102,7 @@ public sealed class HeuristicPrimeTester { if (!EnableHeuristicPrimeTesting) { - return PrimeTester.IsPrimeGpuFallback(n, ct); + return PrimeTester.Exclusive.IsPrimeGpu(n, ct); } return IsPrimeGpu(n); @@ -119,7 +119,7 @@ public sealed class HeuristicPrimeTester { if (!EnableHeuristicPrimeTesting) { - return PrimeTester.IsPrimeGpuFallback(n, CancellationToken.None); + return PrimeTester.Exclusive.IsPrimeGpu(n, CancellationToken.None); } return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10); diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 4e56414..b7c49c8 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -72,18 +72,14 @@ public sealed class PrimeTester return Exclusive.IsPrimeGpu(n, CancellationToken.None); } - public static bool IsPrimeGpu(ulong n, ulong limit, byte nMod10) - { - return Exclusive.IsPrimeGpu(n, CancellationToken.None); - } + // Legacy overload retained for heuristic toggles; tests and benchmarks now call the thread-local instance directly. + // public static bool IsPrimeGpu(ulong n, ulong limit, byte nMod10) + // { + // return Exclusive.IsPrimeGpu(n, CancellationToken.None); + // } [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool IsPrimeGpu(ulong n, CancellationToken ct) - { - return IsPrimeGpuFallback(n, ct); - } - - public static bool IsPrimeGpuFallback(ulong n, CancellationToken ct) { bool forceCpu = GpuContextPool.ForceCpu; Span one = stackalloc ulong[1]; diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md index 46868da..fee5049 100644 --- a/docs/HeuristicIsPrimePlan.md +++ b/docs/HeuristicIsPrimePlan.md @@ -140,7 +140,7 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, * [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] 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. + * [done] Maintain legacy paths for benchmarking via `PrimeTester.IsPrimeInternal`, `PrimeTester.Exclusive.IsPrimeGpu`, 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. diff --git a/docs/HeuristicIsPrimeSpecification.md b/docs/HeuristicIsPrimeSpecification.md index c0282a7..df89eec 100644 --- a/docs/HeuristicIsPrimeSpecification.md +++ b/docs/HeuristicIsPrimeSpecification.md @@ -62,7 +62,7 @@ No other endings are allowed or reordered. ## Default Prime Tester Entry Points * `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. +* Toggle `HeuristicPrimeTester.EnableHeuristicPrimeTesting` or fall back to `PrimeTester.IsPrimeInternal` / `PrimeTester.Exclusive.IsPrimeGpu` 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 From 21d645fcf0d622f3c6e757148742c24f0e142bd3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 18:15:22 +0100 Subject: [PATCH 7/9] Comment out legacy PrimeTester guards --- .../PrimeTester/PrimeTesterIsPrimeTests.cs | 16 ++-- PerfectNumbers.Core/PrimeTester.cs | 86 ++++++++----------- 2 files changed, 45 insertions(+), 57 deletions(-) diff --git a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs index 51e76c4..d33dd32 100644 --- a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs @@ -8,16 +8,16 @@ public class PrimeTesterIsPrimeTests { [Theory] [Trait("Category", "Fast")] - // EvenPerfectBitScanner never routes these small values, but factoring helpers and residue checks still call into - // PrimeTester for sanity probes, so the regression suite keeps the tiny cases alongside large primes/composites. - [InlineData(1UL, false)] - [InlineData(2UL, true)] - [InlineData(3UL, true)] - [InlineData(4UL, false)] - [InlineData(5UL, true)] + // EvenPerfectBitScanner never routes the tiny values anymore; keep coverage focused on the admissible range handled by production scans so the guards in PrimeTester stay commented out without breaking validation. [InlineData(7UL, true)] + [InlineData(9UL, false)] + [InlineData(11UL, true)] + [InlineData(13UL, true)] + [InlineData(21UL, false)] [InlineData(97UL, true)] - [InlineData(121UL, false)] + [InlineData(133UL, false)] + [InlineData(137UL, true)] + [InlineData(341UL, false)] public void IsPrime_returns_expected_results(ulong n, bool expected) { var tester = new PrimeTester(); diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index b7c49c8..91bc672 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -27,16 +27,21 @@ public sealed class PrimeTester // EvenPerfectBitScanner feeds monotonically increasing odd exponents that start well above the // small-prime cutoffs and already exclude multiples of five. These guards remain for targeted // tests and diagnostics so ad-hoc callers can still probe small or even values when needed. - bool isTwo = n == 2UL; - bool isOdd = (n & 1UL) != 0UL; - ulong mod5 = n % 5UL; - bool divisibleByFive = n > 5UL && mod5 == 0UL; + if (n <= 1UL) + { + return false; + } - bool result = n >= 2UL && (isTwo || isOdd) && !divisibleByFive; + // bool isTwo = n == 2UL; + bool isOdd = (n & 1UL) != 0UL; + // ulong mod5 = n % 5UL; + // bool divisibleByFive = n > 5UL && mod5 == 0UL; + + bool result = /*n >= 2UL && (isTwo || */ (n == 2UL || isOdd) /*) && !divisibleByFive*/; // Production scans always fall through to trial division because n is already > 5 by the time this // path executes. Keep the guard so synthetic callers that hammer tiny inputs can skip the loop. - bool requiresTrialDivision = result && n >= 7UL && !isTwo; + bool requiresTrialDivision = result /*&& n >= 7UL && !isTwo*/ && n >= 7UL; if (requiresTrialDivision) { @@ -72,12 +77,6 @@ public sealed class PrimeTester return Exclusive.IsPrimeGpu(n, CancellationToken.None); } - // Legacy overload retained for heuristic toggles; tests and benchmarks now call the thread-local instance directly. - // public static bool IsPrimeGpu(ulong n, ulong limit, byte nMod10) - // { - // return Exclusive.IsPrimeGpu(n, CancellationToken.None); - // } - [MethodImpl(MethodImplOptions.AggressiveInlining)] public bool IsPrimeGpu(ulong n, CancellationToken ct) { @@ -113,52 +112,41 @@ public sealed class PrimeTester { var limiter = GpuPrimeWorkLimiter.Acquire(); var gpu = PrimeTesterGpuContextPool.Rent(); + var accelerator = gpu.Accelerator; + var state = GpuKernelState.GetOrCreate(accelerator); + int totalLength = values.Length; + int batchSize = Math.Max(1, GpuBatchSize); - try + lock (gpu.ExecutionLock) { - var accelerator = gpu.Accelerator; - var state = GpuKernelState.GetOrCreate(accelerator); - int totalLength = values.Length; - int batchSize = Math.Max(1, GpuBatchSize); + var scratch = state.RentScratch(batchSize, accelerator); + var input = scratch.Input; + var output = scratch.Output; + ArrayPool pool = ThreadStaticPools.UlongPool; + ulong[] temp = pool.Rent(batchSize); - lock (gpu.ExecutionLock) + int pos = 0; + while (pos < totalLength) { - var scratch = state.RentScratch(batchSize, accelerator); - var input = scratch.Input; - var output = scratch.Output; - ArrayPool pool = ThreadStaticPools.UlongPool; - ulong[] temp = pool.Rent(batchSize); + int remaining = totalLength - pos; + int count = remaining > batchSize ? batchSize : remaining; - try - { - int pos = 0; - while (pos < totalLength) - { - int remaining = totalLength - pos; - int count = remaining > batchSize ? batchSize : remaining; + values.Slice(pos, count).CopyTo(temp); + input.View.CopyFromCPU(ref temp[0], count); - values.Slice(pos, count).CopyTo(temp); - input.View.CopyFromCPU(ref temp[0], count); + state.Kernel(count, input.View, state.DevicePrimes.View, output.View); + accelerator.Synchronize(); + output.View.CopyToCPU(ref results[pos], count); - state.Kernel(count, input.View, state.DevicePrimes.View, output.View); - accelerator.Synchronize(); - output.View.CopyToCPU(ref results[pos], count); - - pos += count; - } - } - finally - { - pool.Return(temp, clearArray: false); - state.ReturnScratch(scratch); - } + pos += count; } + + pool.Return(temp, clearArray: false); + state.ReturnScratch(scratch); } - finally - { - gpu.Dispose(); - limiter.Dispose(); - } + + gpu.Dispose(); + limiter.Dispose(); } internal static class PrimeTesterGpuContextPool From 1e5eadac3461f1b0d0ceec89b8733e554c0bdb4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 18:15:27 +0100 Subject: [PATCH 8/9] Refine heuristic prime tester integration --- .../PrimeTesterBenchmarks.cs | 8 +- PerfectNumbers.Core/HeuristicPrimeTester.cs | 331 +++++++----------- PerfectNumbers.Core/MersenneDivisorCycles.cs | 10 +- .../MersennePrimeFactorTester.cs | 8 +- .../PrimeOrderCalculator.Cpu.cs | 26 +- .../PrimeOrderCalculator.Shared.cs | 4 +- PerfectNumbers.Core/PrimeTester.cs | 28 +- docs/HeuristicIsPrimePlan.md | 20 +- docs/HeuristicIsPrimeSpecification.md | 12 +- todo_list.txt | 2 +- 10 files changed, 186 insertions(+), 263 deletions(-) diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index d3bae87..e44c7b1 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -32,8 +32,8 @@ public class PrimeTesterBenchmarks } CancellationToken cancellationToken = CancellationToken.None; - _cpuTester.IsPrime(_candidates[0], cancellationToken); - _gpuTester.IsPrimeGpu(_candidates[0], cancellationToken); + _cpuTester.IsPrimeCpu(_candidates[0], cancellationToken); + _gpuTester.IsPrimeGpu(_candidates[0]); _ = Prime.Numbers.IsPrime(_candidates[0]); } @@ -47,7 +47,7 @@ public class PrimeTesterBenchmarks for (int i = 0; i < values.Length; i++) { - if (tester.IsPrime(values[i], cancellationToken)) + if (tester.IsPrimeCpu(values[i], cancellationToken)) { primeCount++; } @@ -105,7 +105,7 @@ public class PrimeTesterBenchmarks for (int i = 0; i < values.Length; i++) { - if (PrimeTester.IsPrimeInternal(values[i], cancellationToken)) + if (PrimeTester.IsPrimeCp(values[i], cancellationToken)) { primeCount++; } diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index ca6e9ea..2cbf881 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -71,7 +71,6 @@ public sealed class HeuristicPrimeTester public readonly bool RequiresCycleComputation = !hasCycleLengthHint; } - public static bool EnableHeuristicPrimeTesting { get; set; } = true; public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096; @@ -82,30 +81,25 @@ public sealed class HeuristicPrimeTester public static HeuristicPrimeTester Exclusive => _tester ??= new(); - public bool IsPrime(ulong n, CancellationToken ct) + public bool IsPrimeCpu(ulong n, CancellationToken ct) { - if (!EnableHeuristicPrimeTesting) + ct.ThrowIfCancellationRequested(); + + byte nMod10 = (byte)n.Mod10(); + ulong sqrtLimit = ComputeHeuristicSqrt(n); + + if (sqrtLimit < 3UL) { - return PrimeTester.IsPrimeInternal(n, ct); + return EvaluateWithOpenNumericFallback(n); } - return HeuristicIsPrimeCpu(n); + bool includeGroupB = UseHeuristicGroupBTrialDivision; + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB); } - public bool IsPrime(ulong n) + public bool IsPrimeCpu(ulong n) { - return IsPrime(n, CancellationToken.None); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool IsPrimeGpu(ulong n, CancellationToken ct) - { - if (!EnableHeuristicPrimeTesting) - { - return PrimeTester.Exclusive.IsPrimeGpu(n, ct); - } - - return IsPrimeGpu(n); + return IsPrimeCpu(n, CancellationToken.None); } public bool IsPrimeGpu(ulong n) @@ -117,45 +111,9 @@ public sealed class HeuristicPrimeTester public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) { - if (!EnableHeuristicPrimeTesting) - { - return PrimeTester.Exclusive.IsPrimeGpu(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? @@ -169,11 +127,6 @@ public sealed class HeuristicPrimeTester return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); } - if (GpuContextPool.ForceCpu) - { - return HeuristicIsPrimeCpu(n, sqrtLimit, nMod10); - } - bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); return !compositeDetected; } @@ -304,6 +257,7 @@ public sealed class HeuristicPrimeTester + private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) { int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); @@ -311,194 +265,165 @@ public sealed class HeuristicPrimeTester var divisorPool = ThreadStaticPools.UlongPool; var hitPool = ThreadStaticPools.BytePool; - ulong[]? divisorArray = null; - byte[]? hitFlags = null; + ulong[] divisorArray = divisorPool.Rent(batchCapacity); + byte[] hitFlags = hitPool.Rent(batchCapacity); var limiter = GpuPrimeWorkLimiter.Acquire(); var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(); + var accelerator = gpu.Accelerator; + var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); + var scratch = state.RentScratch(batchCapacity, accelerator); bool compositeDetected = false; + int count = 0; - try + lock (gpu.ExecutionLock) { - divisorArray = divisorPool.Rent(batchCapacity); - hitFlags = hitPool.Rent(batchCapacity); + var input = scratch.Input; + var output = scratch.Output; - var accelerator = gpu.Accelerator; - var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); - - lock (gpu.ExecutionLock) + bool ProcessBatch(int length) { - var scratch = state.RentScratch(batchCapacity, accelerator); - try + if (length <= 0) { - int count = 0; + return false; + } - 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); + + for (int i = 0; i < length; i++) + { + if (hitFlags[i] == 0) { - 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; + continue; } - bool QueueDivisor(ulong divisor) + ulong divisor = divisorArray[i]; + if (divisor > 1UL && n % divisor == 0UL) { - divisorArray![count] = divisor; + 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]; + + 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; + } + + divisorArray[count] = divisor; count++; if (count == batchCapacity) { if (ProcessBatch(count)) { - return true; + compositeDetected = true; + goto Cleanup; } count = 0; } - return false; + groupAIndex++; + processed++; } - 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 (groupAIndex >= groupADivisors.Length) { - 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); - } - } + groupAHasMore = false; } - - if (!compositeDetected && count > 0) + else { - compositeDetected = ProcessBatch(count); + groupAHasMore = (ulong)groupADivisors[groupAIndex] <= sqrtLimit; + if (!groupAHasMore) + { + groupAIndex = groupADivisors.Length; + } } } - finally + + if (!compositeDetected && groupBHasMore) { - state.ReturnScratch(scratch); + 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); + } } } - } - finally - { - if (divisorArray is not null) - { - divisorPool.Return(divisorArray); - } - if (hitFlags is not null) + if (!compositeDetected && count > 0) { - hitPool.Return(hitFlags); - } + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } - gpu.Dispose(); - limiter.Dispose(); + count = 0; + } } +Cleanup: + state.ReturnScratch(scratch); + divisorPool.Return(divisorArray); + hitPool.Return(hitFlags); + gpu.Dispose(); + limiter.Dispose(); return compositeDetected; } - - private static bool EvaluateWithOpenNumericFallback(ulong n) { return Prime.Numbers.IsPrime(n); diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 4988284..ea03d66 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -394,9 +394,9 @@ public class MersenneDivisorCycles return true; } - // 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)) + // TODO: Benchmark PrimeTester.IsPrimeCp against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. + // if (HeuristicPrimeTester.Exclusive.IsPrimeCpu(remaining, CancellationToken.None)) + if (PrimeTester.IsPrimeCp(remaining, CancellationToken.None)) { AddFactor(counts, remaining); return true; @@ -928,8 +928,8 @@ public class MersenneDivisorCycles return true; } - // if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrime(divisor, CancellationToken.None)) - if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) + // if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None)) + if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeCp(divisor, CancellationToken.None)) { ulong computedOrder = PrimeOrderCalculator.Calculate( divisor, diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs index be53488..d8a63cc 100644 --- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs +++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs @@ -25,8 +25,8 @@ public static class MersennePrimeFactorTester return false; } - // if (!HeuristicPrimeTester.Exclusive.IsPrime(q, ct)) - if (!PrimeTester.IsPrimeInternal(q, ct)) + // if (!HeuristicPrimeTester.Exclusive.IsPrimeCpu(q, ct)) + if (!PrimeTester.IsPrimeCp(q, ct)) { return false; } @@ -234,8 +234,8 @@ public static class MersennePrimeFactorTester continue; } - // if (HeuristicPrimeTester.Exclusive.IsPrime(m, ct)) - if (PrimeTester.IsPrimeInternal(m, ct)) + // if (HeuristicPrimeTester.Exclusive.IsPrimeCpu(m, ct)) + if (PrimeTester.IsPrimeCp(m, ct)) { dict[m] = 1; continue; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index e24d625..c85e6bf 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -516,7 +516,7 @@ internal static partial class PrimeOrderCalculator return false; } - // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None)) + // if (!PrimeTester.IsPrimeCp(orderFactors.Cofactor, CancellationToken.None)) if (!orderFactors.CofactorIsPrime) { return false; @@ -1123,8 +1123,8 @@ internal static partial class PrimeOrderCalculator // } // bool isPrime = primeTester.HeuristicIsPrimeGpu(composite); - // bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); + // bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) @@ -1186,8 +1186,8 @@ internal static partial class PrimeOrderCalculator if (!entry.HasKnownPrimality) { // bool isPrime = primeTester.HeuristicIsPrimeGpu(composite); - // bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); + // bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); @@ -1219,8 +1219,8 @@ internal static partial class PrimeOrderCalculator else { // cofactorIsPrime = primeTester.HeuristicIsPrimeGpu(cofactor); - // cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None); - cofactorIsPrime = PrimeTester.IsPrimeInternal(cofactor, CancellationToken.None); + // cofactorIsPrime = PrimeTester.IsPrimeCp(cofactor, CancellationToken.None); + cofactorIsPrime = PrimeTester.IsPrimeCp(cofactor, CancellationToken.None); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1752,8 +1752,8 @@ internal static partial class PrimeOrderCalculator PrimeTester primeTester = PrimeTester.Exclusive; // bool isPrime = primeTester.HeuristicIsPrimeGpu(value); - // bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeInternal(value, CancellationToken.None); + // bool isPrime = PrimeTester.IsPrimeCp(value, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeCp(value, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value); if (!knownComposite && isPrime) { @@ -1765,8 +1765,8 @@ internal static partial class PrimeOrderCalculator ulong quotient = value / factor; // isPrime = primeTester.HeuristicIsPrimeGpu(factor); - // isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None); - isPrime = PrimeTester.IsPrimeInternal(factor, CancellationToken.None); + // isPrime = PrimeTester.IsPrimeCp(factor, CancellationToken.None); + isPrime = PrimeTester.IsPrimeCp(factor, CancellationToken.None); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor); if (isPrime) { @@ -1797,8 +1797,8 @@ internal static partial class PrimeOrderCalculator else { // isPrime = primeTester.HeuristicIsPrimeGpu(quotient); - // isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None); - isPrime = PrimeTester.IsPrimeInternal(quotient, CancellationToken.None); + // isPrime = PrimeTester.IsPrimeCp(quotient, CancellationToken.None); + isPrime = PrimeTester.IsPrimeCp(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 bd81fbb..ee4b914 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -955,8 +955,8 @@ internal static partial class PrimeOrderCalculator { if (value <= ulong.MaxValue) { - // return HeuristicPrimeTester.Exclusive.IsPrime((ulong)value, CancellationToken.None); - return PrimeTester.IsPrimeInternal((ulong)value, CancellationToken.None); + // return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None); + return PrimeTester.IsPrimeCp((ulong)value, CancellationToken.None); } if ((value & UInt128.One) == UInt128.Zero) diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 91bc672..31933a0 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -24,29 +24,27 @@ public sealed class PrimeTester public bool IsPrime(ulong n, CancellationToken ct) { - // EvenPerfectBitScanner feeds monotonically increasing odd exponents that start well above the - // small-prime cutoffs and already exclude multiples of five. These guards remain for targeted - // tests and diagnostics so ad-hoc callers can still probe small or even values when needed. if (n <= 1UL) { return false; } - // bool isTwo = n == 2UL; + if (n == 2UL) + { + throw new InvalidOperationException("PrimeTester.IsPrime encountered the sentinel input 2."); + } + bool isOdd = (n & 1UL) != 0UL; - // ulong mod5 = n % 5UL; - // bool divisibleByFive = n > 5UL && mod5 == 0UL; + bool result = isOdd; - bool result = /*n >= 2UL && (isTwo || */ (n == 2UL || isOdd) /*) && !divisibleByFive*/; - - // Production scans always fall through to trial division because n is already > 5 by the time this - // path executes. Keep the guard so synthetic callers that hammer tiny inputs can skip the loop. - bool requiresTrialDivision = result /*&& n >= 7UL && !isTwo*/ && n >= 7UL; + bool requiresTrialDivision = result && n >= 7UL; if (requiresTrialDivision) { - bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); - result &= !sharesMaxExponentFactor; + // EvenPerfectBitScanner streams exponents starting at 136,279,841, so the Mod10/GCD guard never fires on the + // production path. Leave the logic commented out as instrumentation for diagnostic builds. + // bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); + // result &= !sharesMaxExponentFactor; if (result) { @@ -95,10 +93,10 @@ public sealed class PrimeTester bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0; bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime; - return requiresCpuFallback ? IsPrimeInternal(n, ct) : true; + return requiresCpuFallback ? IsPrimeCp(n, ct) : true; } - public static bool IsPrimeInternal(ulong n, CancellationToken ct) + public static bool IsPrimeCp(ulong n, CancellationToken ct) { // Preserve the legacy entry point for callers that bypass the thread-local caches. return Exclusive.IsPrime(n, ct); diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md index fee5049..ed33f39 100644 --- a/docs/HeuristicIsPrimePlan.md +++ b/docs/HeuristicIsPrimePlan.md @@ -64,12 +64,12 @@ public bool IsPrimeGpu(ulong n, CancellationToken ct) bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0; bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime; - return requiresCpuFallback ? IsPrimeInternal(n, ct) : true; + return requiresCpuFallback ? IsPrimeCp(n, ct) : true; } ``` ```csharp -internal static bool IsPrimeInternal(ulong n, CancellationToken ct) +internal static bool IsPrimeCp(ulong n, CancellationToken ct) { bool isTwo = n == 2UL; bool isOdd = (n & 1UL) != 0UL; @@ -113,11 +113,11 @@ internal static bool IsPrimeInternal(ulong n, CancellationToken ct) 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]`. -1. [done] **Snapshot current implementations** – preserve `IsPrimeGpu` and `IsPrimeInternal` as references for regression comparisons and migration checkpoints (see Section 2). -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. +1. [done] **Snapshot current implementations** – preserve `IsPrimeGpu` and `IsPrimeCp` as references for regression comparisons and migration checkpoints (see Section 2). +2. [done] **Introduce placeholder heuristic entry points** – add `IsPrimeCpu` and `IsPrimeGpu` 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`; HeuristicPrimeTester 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 `IsPrimeCpu`; 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. @@ -126,11 +126,11 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, 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. -6. [done] **Implement `HeuristicIsPrimeCpu`** - * Baseline Group A/B enumeration now executes inside `HeuristicIsPrimeCpu` via the cached `HeuristicPrimeSieves` tables, computing `⌊√n⌋` locally and short-circuiting on the first divisor. +6. [done] **Implement `IsPrimeCpu`** + * Baseline Group A/B enumeration now executes inside `IsPrimeCpu` via the cached `HeuristicPrimeSieves` tables, computing `⌊√n⌋` locally and short-circuiting on the first divisor. * [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 `HeuristicIsPrimeGpu`** +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] 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. @@ -140,12 +140,12 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, * [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] 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.Exclusive.IsPrimeGpu`, and the global `HeuristicPrimeTester.EnableHeuristicPrimeTesting` switch so benchmark harnesses can compare behaviours. + * [done] Maintain legacy paths for benchmarking via `PrimeTester.IsPrimeCp`, `PrimeTester.Exclusive.IsPrimeGpu`, and the explicit `HeuristicPrimeTester` entry points so benchmark harnesses can compare behaviours without runtime switches. 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 `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. + * [done] Build BenchmarkDotNet comparisons for `HeuristicPrimeTester.IsPrimeGpu`, `HeuristicPrimeTester.IsPrimeCpu`, 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 df89eec..2d02573 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 `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. +This document describes the heuristic prime testing procedure implemented by `HeuristicPrimeTester.IsPrimeCpu` 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 @@ -57,23 +57,23 @@ No other endings are allowed or reordered. 4. **Group B**: read the precomputed `n mod 10` (made available by the caller configuration such as `--gpu-prime-batch`). Iterate candidates via wheel 210, maintaining ascending order and applying the ending-priority sequence above. Any hit marks `n` as composite and terminates the search. 5. **Decision**: if no divisors were found up to `R`, report `n` as “prime by trial”. This outcome is final—no further confirmation runs are performed. -> **Temporary status:** Current builds disable Group B via `UseHeuristicGroupBTrialDivision`. After exhausting Group A, `HeuristicIsPrime*` delegates to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, preserving the Group B implementation for future use while ensuring production checks mirror the external library once Group A filters pass. +> **Temporary status:** Current builds disable Group B via `UseHeuristicGroupBTrialDivision`. After exhausting Group A, `HeuristicPrimeTester.IsPrime*` delegates to `Open.Numeric.Primes.Prime.Numbers.IsPrime`, preserving the Group B implementation for future use while ensuring production checks mirror the external library once Group A filters pass. ## Default Prime Tester Entry Points * `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.Exclusive.IsPrimeGpu` to control heuristic rollouts in test harnesses. +* `HeuristicPrimeTester` always executes the heuristic CPU/GPU implementations; opt into the legacy paths by calling `PrimeTester.IsPrimeCp` or `PrimeTester.Exclusive.IsPrimeGpu` directly. * Benchmark harnesses can compare behaviours by calling the legacy helpers on `PrimeTester` alongside the heuristic equivalents on `HeuristicPrimeTester`. ## Integration with the By-Divisor Pipeline The by-divisor executors (`MersenneNumberDivisorByDivisorCpuTester` and `MersenneNumberDivisorByDivisorGpuTester`) reuse prime heuristics to seed their divisor scans, residue sieves, and cycle computations. The heuristic implementations must therefore expose the following data and behaviours: -* **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. +* **Residue filters** – The CPU scan maintains rolling residues modulo 10/8/5/3/7/11 to match the filters in `MersenneCpuDivisorScanSession.CheckDivisors`. `HeuristicPrimeTester.IsPrime*` 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** – `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. +* **GPU batching** – `HeuristicPrimeTester.IsPrimeGpu` 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. The implementation now assumes the accelerator path is available and no longer routes through the CPU fallback. * **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 `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. +Benchmark results comparing `HeuristicPrimeTester.IsPrimeGpu`, `HeuristicPrimeTester.IsPrimeCpu`, 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. diff --git a/todo_list.txt b/todo_list.txt index 9eb6432..c65c9e4 100644 --- a/todo_list.txt +++ b/todo_list.txt @@ -18,7 +18,7 @@ Divisor Cycle Cache & Generation - PerfectNumbers.Core/UInt128Extensions.cs:372: If the cache lacks this cycle, immediately schedule the configured device (GPU by default) to compute it on the fly and skip inserting it into the cache so wide-order factoring can still leverage cycle stepping without breaching the memory cap or introducing extra synchronization. Primality & Order Screening -- PerfectNumbers.Core/MersenneDivisorCycles.cs:477: Benchmark PrimeTester.IsPrimeInternal against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. +- PerfectNumbers.Core/MersenneDivisorCycles.cs:477: Benchmark PrimeTester.IsPrimeCp against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. - PerfectNumbers.Core/ULongExtensions.cs:119: When the shared cycle snapshot cannot serve this divisor, trigger an on-demand GPU computation (respecting the configured device) without promoting the result into the cache so the order calculator still benefits from cycle stepping while keeping the single-block memory plan intact. - PerfectNumbers.Core/ULongExtensions.cs:132: Replace this `%` driven factor peeling with the divisor-cycle aware factoring helper so large orders reuse the cached remainders highlighted in the latest divisor-cycle benchmarks instead of recomputing slow modulo checks. - PerfectNumbers.Core/ULongExtensions.cs:279: Swap this modulo check for the shared small-prime cycle filter once the divisor-cycle cache is mandatory, matching the PrimeTester improvements noted in the CPU sieve benchmarks. From 37585244500061ff517ab2f3aa688038dedc4051 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Sun, 26 Oct 2025 18:59:07 +0100 Subject: [PATCH 9/9] Extract divisor validator and rename legacy CPU prime helper --- .../PrimeTesterBenchmarks.cs | 2 +- PerfectNumbers.Core/DivisorGenerator.cs | 38 +++++++++++++++++++ ...MersenneNumberDivisorByDivisorGpuTester.cs | 34 +---------------- PerfectNumbers.Core/MersenneDivisorCycles.cs | 6 +-- .../MersennePrimeFactorTester.cs | 4 +- .../PrimeOrderCalculator.Cpu.cs | 26 ++++++------- .../PrimeOrderCalculator.Shared.cs | 2 +- PerfectNumbers.Core/PrimeTester.cs | 4 +- docs/HeuristicIsPrimePlan.md | 8 ++-- docs/HeuristicIsPrimeSpecification.md | 2 +- todo_list.txt | 2 +- 11 files changed, 67 insertions(+), 61 deletions(-) create mode 100644 PerfectNumbers.Core/DivisorGenerator.cs diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index e44c7b1..04fb367 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -105,7 +105,7 @@ public class PrimeTesterBenchmarks for (int i = 0; i < values.Length; i++) { - if (PrimeTester.IsPrimeCp(values[i], cancellationToken)) + if (PrimeTester.IsPrimeCpu(values[i], cancellationToken)) { primeCount++; } diff --git a/PerfectNumbers.Core/DivisorGenerator.cs b/PerfectNumbers.Core/DivisorGenerator.cs new file mode 100644 index 0000000..0001267 --- /dev/null +++ b/PerfectNumbers.Core/DivisorGenerator.cs @@ -0,0 +1,38 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +internal static class DivisorGenerator +{ + private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); + private const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool IsValidDivisor( + byte remainder10, + byte remainder8, + byte remainder3, + byte remainder5, + byte remainder7, + byte remainder11, + LastDigit lastDigit) + { + ushort decimalMask = lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; + if (((decimalMask >> remainder10) & 1) == 0) + { + return false; + } + + if (remainder8 != 1 && remainder8 != 7) + { + return false; + } + + if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0) + { + return false; + } + + return true; + } +} diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index e1b1fe5..8d83ee2 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -263,7 +263,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { ulong candidate = (ulong)nextDivisor128; - if (CandidatePassesHeuristics(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastDigit)) + if (DivisorGenerator.IsValidDivisor(localRemainder10, localRemainder8, localRemainder3, localRemainder5, localRemainder7, localRemainder11, lastDigit)) { filteredStorage[filteredCount++] = candidate; } @@ -378,38 +378,6 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv return computedCycle; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool CandidatePassesHeuristics( - byte remainder10, - byte remainder8, - byte remainder3, - byte remainder5, - byte remainder7, - byte remainder11, - LastDigit lastDigit) - { - const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); - const ushort DecimalMaskOtherwise = (1 << 1) | (1 << 3) | (1 << 9); - - ushort decimalMask = lastDigit == LastDigit.Seven ? DecimalMaskWhenLastIsSeven : DecimalMaskOtherwise; - if (((decimalMask >> remainder10) & 1) == 0) - { - return false; - } - - if (remainder8 != 1 && remainder8 != 7) - { - return false; - } - - if (remainder3 == 0 || remainder5 == 0 || remainder7 == 0 || remainder11 == 0) - { - return false; - } - - return true; - } - [MethodImpl(MethodImplOptions.AggressiveInlining)] private static byte AddMod(byte value, byte delta, byte modulus) { diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index ea03d66..516c512 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -394,9 +394,9 @@ public class MersenneDivisorCycles return true; } - // TODO: Benchmark PrimeTester.IsPrimeCp against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. + // TODO: Benchmark PrimeTester.IsPrimeCpu against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. // if (HeuristicPrimeTester.Exclusive.IsPrimeCpu(remaining, CancellationToken.None)) - if (PrimeTester.IsPrimeCp(remaining, CancellationToken.None)) + if (PrimeTester.IsPrimeCpu(remaining, CancellationToken.None)) { AddFactor(counts, remaining); return true; @@ -929,7 +929,7 @@ public class MersenneDivisorCycles } // if (!skipPrimeOrderHeuristic && HeuristicPrimeTester.Exclusive.IsPrimeCpu(divisor, CancellationToken.None)) - if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeCp(divisor, CancellationToken.None)) + if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeCpu(divisor, CancellationToken.None)) { ulong computedOrder = PrimeOrderCalculator.Calculate( divisor, diff --git a/PerfectNumbers.Core/MersennePrimeFactorTester.cs b/PerfectNumbers.Core/MersennePrimeFactorTester.cs index d8a63cc..0bf2818 100644 --- a/PerfectNumbers.Core/MersennePrimeFactorTester.cs +++ b/PerfectNumbers.Core/MersennePrimeFactorTester.cs @@ -26,7 +26,7 @@ public static class MersennePrimeFactorTester } // if (!HeuristicPrimeTester.Exclusive.IsPrimeCpu(q, ct)) - if (!PrimeTester.IsPrimeCp(q, ct)) + if (!PrimeTester.IsPrimeCpu(q, ct)) { return false; } @@ -235,7 +235,7 @@ public static class MersennePrimeFactorTester } // if (HeuristicPrimeTester.Exclusive.IsPrimeCpu(m, ct)) - if (PrimeTester.IsPrimeCp(m, ct)) + if (PrimeTester.IsPrimeCpu(m, ct)) { dict[m] = 1; continue; diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index c85e6bf..00981e2 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -516,7 +516,7 @@ internal static partial class PrimeOrderCalculator return false; } - // if (!PrimeTester.IsPrimeCp(orderFactors.Cofactor, CancellationToken.None)) + // if (!PrimeTester.IsPrimeCpu(orderFactors.Cofactor, CancellationToken.None)) if (!orderFactors.CofactorIsPrime) { return false; @@ -1123,8 +1123,8 @@ internal static partial class PrimeOrderCalculator // } // bool isPrime = primeTester.HeuristicIsPrimeGpu(composite); - // bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); + // bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); if (isPrime) @@ -1186,8 +1186,8 @@ internal static partial class PrimeOrderCalculator if (!entry.HasKnownPrimality) { // bool isPrime = primeTester.HeuristicIsPrimeGpu(composite); - // bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeCp(composite, CancellationToken.None); + // bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); entry = entry.WithPrimality(isPrime); @@ -1219,8 +1219,8 @@ internal static partial class PrimeOrderCalculator else { // cofactorIsPrime = primeTester.HeuristicIsPrimeGpu(cofactor); - // cofactorIsPrime = PrimeTester.IsPrimeCp(cofactor, CancellationToken.None); - cofactorIsPrime = PrimeTester.IsPrimeCp(cofactor, CancellationToken.None); + // cofactorIsPrime = PrimeTester.IsPrimeCpu(cofactor, CancellationToken.None); + cofactorIsPrime = PrimeTester.IsPrimeCpu(cofactor, CancellationToken.None); // cofactorIsPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(cofactor); } @@ -1752,8 +1752,8 @@ internal static partial class PrimeOrderCalculator PrimeTester primeTester = PrimeTester.Exclusive; // bool isPrime = primeTester.HeuristicIsPrimeGpu(value); - // bool isPrime = PrimeTester.IsPrimeCp(value, CancellationToken.None); - bool isPrime = PrimeTester.IsPrimeCp(value, CancellationToken.None); + // bool isPrime = PrimeTester.IsPrimeCpu(value, CancellationToken.None); + bool isPrime = PrimeTester.IsPrimeCpu(value, CancellationToken.None); // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(value); if (!knownComposite && isPrime) { @@ -1765,8 +1765,8 @@ internal static partial class PrimeOrderCalculator ulong quotient = value / factor; // isPrime = primeTester.HeuristicIsPrimeGpu(factor); - // isPrime = PrimeTester.IsPrimeCp(factor, CancellationToken.None); - isPrime = PrimeTester.IsPrimeCp(factor, CancellationToken.None); + // isPrime = PrimeTester.IsPrimeCpu(factor, CancellationToken.None); + isPrime = PrimeTester.IsPrimeCpu(factor, CancellationToken.None); // isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factor); if (isPrime) { @@ -1797,8 +1797,8 @@ internal static partial class PrimeOrderCalculator else { // isPrime = primeTester.HeuristicIsPrimeGpu(quotient); - // isPrime = PrimeTester.IsPrimeCp(quotient, CancellationToken.None); - isPrime = PrimeTester.IsPrimeCp(quotient, CancellationToken.None); + // isPrime = PrimeTester.IsPrimeCpu(quotient, CancellationToken.None); + isPrime = PrimeTester.IsPrimeCpu(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 ee4b914..24ba684 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -956,7 +956,7 @@ internal static partial class PrimeOrderCalculator if (value <= ulong.MaxValue) { // return HeuristicPrimeTester.Exclusive.IsPrimeCpu((ulong)value, CancellationToken.None); - return PrimeTester.IsPrimeCp((ulong)value, CancellationToken.None); + return PrimeTester.IsPrimeCpu((ulong)value, CancellationToken.None); } if ((value & UInt128.One) == UInt128.Zero) diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 31933a0..dd15c4e 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -93,10 +93,10 @@ public sealed class PrimeTester bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0; bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime; - return requiresCpuFallback ? IsPrimeCp(n, ct) : true; + return requiresCpuFallback ? IsPrimeCpu(n, ct) : true; } - public static bool IsPrimeCp(ulong n, CancellationToken ct) + public static bool IsPrimeCpu(ulong n, CancellationToken ct) { // Preserve the legacy entry point for callers that bypass the thread-local caches. return Exclusive.IsPrime(n, ct); diff --git a/docs/HeuristicIsPrimePlan.md b/docs/HeuristicIsPrimePlan.md index ed33f39..fd6671e 100644 --- a/docs/HeuristicIsPrimePlan.md +++ b/docs/HeuristicIsPrimePlan.md @@ -64,12 +64,12 @@ public bool IsPrimeGpu(ulong n, CancellationToken ct) bool gpuReportedPrime = !forceCpu && !belowGpuRange && outFlags[0] != 0; bool requiresCpuFallback = forceCpu || belowGpuRange || !gpuReportedPrime; - return requiresCpuFallback ? IsPrimeCp(n, ct) : true; + return requiresCpuFallback ? IsPrimeCpu(n, ct) : true; } ``` ```csharp -internal static bool IsPrimeCp(ulong n, CancellationToken ct) +internal static bool IsPrimeCpu(ulong n, CancellationToken ct) { bool isTwo = n == 2UL; bool isOdd = (n & 1UL) != 0UL; @@ -113,7 +113,7 @@ internal static bool IsPrimeCp(ulong n, CancellationToken ct) 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]`. -1. [done] **Snapshot current implementations** – preserve `IsPrimeGpu` and `IsPrimeCp` as references for regression comparisons and migration checkpoints (see Section 2). +1. [done] **Snapshot current implementations** – preserve `IsPrimeGpu` and `IsPrimeCpu` as references for regression comparisons and migration checkpoints (see Section 2). 2. [done] **Introduce placeholder heuristic entry points** – add `IsPrimeCpu` and `IsPrimeGpu` 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. @@ -140,7 +140,7 @@ The plan incorporates the updated divisor-class heuristics (Groups A/B, wheels, * [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] 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.IsPrimeCp`, `PrimeTester.Exclusive.IsPrimeGpu`, and the explicit `HeuristicPrimeTester` entry points so benchmark harnesses can compare behaviours without runtime switches. + * [done] Maintain legacy paths for benchmarking via `PrimeTester.IsPrimeCpu`, `PrimeTester.Exclusive.IsPrimeGpu`, and the explicit `HeuristicPrimeTester` entry points so benchmark harnesses can compare behaviours without runtime switches. 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. diff --git a/docs/HeuristicIsPrimeSpecification.md b/docs/HeuristicIsPrimeSpecification.md index 2d02573..6d3d374 100644 --- a/docs/HeuristicIsPrimeSpecification.md +++ b/docs/HeuristicIsPrimeSpecification.md @@ -62,7 +62,7 @@ No other endings are allowed or reordered. ## Default Prime Tester Entry Points * `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. -* `HeuristicPrimeTester` always executes the heuristic CPU/GPU implementations; opt into the legacy paths by calling `PrimeTester.IsPrimeCp` or `PrimeTester.Exclusive.IsPrimeGpu` directly. +* `HeuristicPrimeTester` always executes the heuristic CPU/GPU implementations; opt into the legacy paths by calling `PrimeTester.IsPrimeCpu` or `PrimeTester.Exclusive.IsPrimeGpu` directly. * Benchmark harnesses can compare behaviours by calling the legacy helpers on `PrimeTester` alongside the heuristic equivalents on `HeuristicPrimeTester`. ## Integration with the By-Divisor Pipeline diff --git a/todo_list.txt b/todo_list.txt index c65c9e4..72016e5 100644 --- a/todo_list.txt +++ b/todo_list.txt @@ -18,7 +18,7 @@ Divisor Cycle Cache & Generation - PerfectNumbers.Core/UInt128Extensions.cs:372: If the cache lacks this cycle, immediately schedule the configured device (GPU by default) to compute it on the fly and skip inserting it into the cache so wide-order factoring can still leverage cycle stepping without breaching the memory cap or introducing extra synchronization. Primality & Order Screening -- PerfectNumbers.Core/MersenneDivisorCycles.cs:477: Benchmark PrimeTester.IsPrimeCp against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. +- PerfectNumbers.Core/MersenneDivisorCycles.cs:477: Benchmark PrimeTester.IsPrimeCpu against Open.Numeric.Primes on the EvenPerfectBitScanner workloads to confirm the faster option for wide composites. - PerfectNumbers.Core/ULongExtensions.cs:119: When the shared cycle snapshot cannot serve this divisor, trigger an on-demand GPU computation (respecting the configured device) without promoting the result into the cache so the order calculator still benefits from cycle stepping while keeping the single-block memory plan intact. - PerfectNumbers.Core/ULongExtensions.cs:132: Replace this `%` driven factor peeling with the divisor-cycle aware factoring helper so large orders reuse the cached remainders highlighted in the latest divisor-cycle benchmarks instead of recomputing slow modulo checks. - PerfectNumbers.Core/ULongExtensions.cs:279: Swap this modulo check for the shared small-prime cycle filter once the divisor-cycle cache is mandatory, matching the PrimeTester improvements noted in the CPU sieve benchmarks.