diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index 5b9bf5c..ac4bd39 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -70,8 +70,17 @@ internal static class PrimeTesterKernels public static void HeuristicTrialDivisionKernel(Index1D index, ArrayView divisors, ulong n, ArrayView resultFlag) { - ulong divisor = divisors[index]; - Atomic.Or(ref resultFlag[0], divisor <= 1UL ? 0 : (n % divisor == 0UL ? 1 : 0)); + ulong divisor = divisors[index]; + if (divisor <= 1UL) + { + return; + } + + ulong residue = n % divisor; + if (residue == 0UL) + { + resultFlag[0] = 1; + } } public static void SharesFactorKernel(Index1D index, ArrayView numbers, ArrayView results) diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 72add2f..24d52a7 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -1,11 +1,5 @@ -using System; 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; @@ -48,8 +42,6 @@ public sealed class HeuristicCombinedPrimeTester private const ulong Wheel210 = 210UL; - private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. - private static readonly uint[] CombinedDivisorsEnding1TwoAOneB; private static readonly uint[] CombinedDivisorsEnding3TwoAOneB; private static readonly uint[] CombinedDivisorsEnding7TwoAOneB; @@ -141,8 +133,6 @@ public sealed class HeuristicCombinedPrimeTester } - internal static int HeuristicGpuDivisorBatchSize = 4_096; - internal static int HeuristicDivisorInterleaveBatchSize = 64; [ThreadStatic] private static HeuristicCombinedPrimeTester? _tester; @@ -205,76 +195,84 @@ public sealed class HeuristicCombinedPrimeTester return EvaluateWithOpenNumericFallback(n); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static int GetBatchSize(ulong n) => n switch + { + < 1_000 => 8, + < 10_000 => 16, + < 100_000 => 32, + < 1_000_000 => 128, + _ => 256, + }; + private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10) { - int batchCapacity = HeuristicGpuDivisorBatchSize; - var divisorPool = ThreadStaticPools.UlongPool; - - ulong[] divisorArray = divisorPool.Rent(batchCapacity); - - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity); - var accelerator = gpu.Accelerator; - var kernel = gpu.HeuristicTrialDivisionKernel; - - bool compositeDetected = false; - int count = 0; - - var inputView = gpu.Input.View; - var flagView = gpu.HeuristicFlag.View; - - bool ProcessBatch(int length) - { - inputView.CopyFromCPU(ref divisorArray[0], length); - int compositeFlag = 0; - flagView.CopyFromCPU(ref compositeFlag, 1); - kernel(length, inputView, n, flagView); - accelerator.Synchronize(); - flagView.CopyToCPU(ref compositeFlag, 1); - - return compositeFlag != 0; - } - - ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); - ReadOnlySpan combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); - - int length = combinedDivisors.Length; - - for (int i = 0; i < length; i++) - { - ulong divisorSquare = combinedDivisorSquares[i]; - if (divisorSquare > maxDivisorSquare) - { - break; - } - - divisorArray[count] = combinedDivisors[i]; - count++; - - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; - } - } - - if (count > 0 && ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } + int batchCapacity = GetBatchSize(n); -Cleanup: - divisorPool.Return(divisorArray); - gpu.Dispose(); - limiter.Dispose(); - return compositeDetected; + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity); + var accelerator = gpu.Accelerator; + var kernel = gpu.HeuristicTrialDivisionKernel; + + bool compositeDetected = false; + int count = 0; + + var inputView = gpu.Input.View; + var flagView = gpu.HeuristicFlag.View; + ulong[] divisorArray = gpu.DivisorArray; + + bool ProcessBatch(int length) + { + inputView.CopyFromCPU(ref divisorArray[0], length); + int compositeFlag = 0; + flagView.CopyFromCPU(ref compositeFlag, 1); + kernel(length, inputView, n, flagView); + accelerator.Synchronize(); + flagView.CopyToCPU(ref compositeFlag, 1); + + return compositeFlag != 0; + } + + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + ReadOnlySpan combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); + + int length = combinedDivisors.Length; + + for (int i = 0; i < length; i++) + { + ulong divisorSquare = combinedDivisorSquares[i]; + if (divisorSquare > maxDivisorSquare) + { + break; + } + + divisorArray[count] = combinedDivisors[i]; + count++; + + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + + count = 0; + } + } + + if (count > 0 && ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + + + Cleanup: + gpu.Dispose(); + limiter.Dispose(); + return compositeDetected; } @@ -282,10 +280,10 @@ Cleanup: { return Prime.Numbers.IsPrime(n); } - internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n) - { - return n; - } + internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n) + { + return n; + } private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { @@ -451,7 +449,7 @@ Cleanup: } } - uint[] result = combined.ToArray(); + uint[] result = [.. combined]; ulong[] resultSquares = new ulong[result.Length]; for (int i = 0; i < result.Length; i++) { @@ -488,10 +486,10 @@ Cleanup: private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch { 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, - 3 => ReadOnlySpan.Empty, + 3 => [], 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, - _ => ReadOnlySpan.Empty, + _ => [], }; [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -578,50 +576,50 @@ Cleanup: } internal ref struct HeuristicDivisorEnumerator - { - private readonly ulong maxDivisorSquare; - private readonly ReadOnlySpan combinedDivisors; - private readonly ReadOnlySpan combinedDivisorSquares; - private int index; + { + private readonly ulong maxDivisorSquare; + private readonly ReadOnlySpan combinedDivisors; + private readonly ReadOnlySpan combinedDivisorSquares; + private int index; - public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) - { - this.maxDivisorSquare = maxDivisorSquare; - _ = groupBBuffer; - combinedDivisors = GetCombinedDivisors(nMod10); - combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); - index = 0; - } + public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) + { + this.maxDivisorSquare = maxDivisorSquare; + _ = groupBBuffer; + combinedDivisors = GetCombinedDivisors(nMod10); + combinedDivisorSquares = GetCombinedDivisorSquares(nMod10); + index = 0; + } - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - while (index < combinedDivisors.Length) - { - int currentIndex = index; - ulong divisorSquare = combinedDivisorSquares[currentIndex]; - if (divisorSquare > maxDivisorSquare) - { - index = combinedDivisors.Length; - break; - } + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (index < combinedDivisors.Length) + { + int currentIndex = index; + ulong divisorSquare = combinedDivisorSquares[currentIndex]; + if (divisorSquare > maxDivisorSquare) + { + index = combinedDivisors.Length; + break; + } - uint entry = combinedDivisors[currentIndex]; - index = currentIndex + 1; + uint entry = combinedDivisors[currentIndex]; + index = currentIndex + 1; - HeuristicDivisorGroup group = ResolveGroup(entry); - ulong value = entry; - byte ending = (byte)(value % 10UL); - ushort residue = (ushort)(value % Wheel210); - candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); - return true; - } + HeuristicDivisorGroup group = ResolveGroup(entry); + ulong value = entry; + byte ending = (byte)(value % 10UL); + ushort residue = (ushort)(value % Wheel210); + candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); + return true; + } - candidate = default; - return false; - } - } + candidate = default; + return false; + } + } - internal ref struct MersenneHeuristicDivisorEnumerator + internal ref struct MersenneHeuristicDivisorEnumerator { private readonly GpuUInt128 step; private readonly GpuUInt128 limit; @@ -690,11 +688,11 @@ Cleanup: residueStepper.Advance(); } - public ulong ProcessedCount => processedCount; + public readonly ulong ProcessedCount => processedCount; - public ulong LastDivisor => lastDivisor; + public readonly ulong LastDivisor => lastDivisor; - public bool Exhausted => !active; + public readonly bool Exhausted => !active; private static HeuristicDivisorCandidate CreateCandidate(ulong value) { diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index 2f491df..930de3d 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -1,13 +1,5 @@ -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; @@ -15,828 +7,809 @@ namespace PerfectNumbers.Core; public sealed class HeuristicPrimeTester { - static HeuristicPrimeTester() - { - HeuristicPrimeSieves.EnsureInitialized(); - } + static HeuristicPrimeTester() + { + HeuristicPrimeSieves.EnsureInitialized(); + } - // TODO: Don't use these tiny arrays. Hard-code the values / checks instead - 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: Don't use these tiny arrays. Hard-code the values / checks instead + 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 const ulong Wheel210 = 210UL; - private const int GroupAConstantCount = 4; - internal const int MaxGroupBSequences = 36; + // TODO: Would it speed up things if we implement something like ExponentRemainderStepper, + // MersenneDivisorResidueStepper, or CycleRemainderStepper work with these? + private const ulong Wheel210 = 210UL; + private const int GroupAConstantCount = 4; + internal const int MaxGroupBSequences = 36; - private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. + private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. - private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - internal enum HeuristicDivisorGroup : byte - { - None = 0, - GroupAConstant = 1, - GroupAWheel = 2, - GroupB = 3, - } + internal enum HeuristicDivisorGroup : byte + { + None = 0, + GroupAConstant = 1, + GroupAWheel = 2, + GroupB = 3, + } - internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue) - { - public readonly ulong Value = value; - public readonly HeuristicDivisorGroup Group = group; - public readonly byte Ending = ending; - public readonly byte PriorityIndex = priorityIndex; - public readonly ushort WheelResidue = wheelResidue; - } + internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue) + { + public readonly ulong Value = value; + public readonly HeuristicDivisorGroup Group = group; + public readonly byte Ending = ending; + public readonly byte PriorityIndex = priorityIndex; + public readonly ushort WheelResidue = wheelResidue; + } - internal readonly struct HeuristicDivisorPreparation( - in HeuristicDivisorCandidate candidate, - in MontgomeryDivisorData divisorData, - ulong cycleLengthHint, - bool hasCycleLengthHint) - { - public readonly HeuristicDivisorCandidate Candidate = candidate; - public readonly MontgomeryDivisorData DivisorData = divisorData; - public readonly ulong CycleLengthHint = cycleLengthHint; - public readonly bool HasCycleLengthHint = hasCycleLengthHint; - public readonly bool RequiresCycleComputation = !hasCycleLengthHint; - } + internal readonly struct HeuristicDivisorPreparation( + in HeuristicDivisorCandidate candidate, + in MontgomeryDivisorData divisorData, + ulong cycleLengthHint, + bool hasCycleLengthHint) + { + public readonly HeuristicDivisorCandidate Candidate = candidate; + public readonly MontgomeryDivisorData DivisorData = divisorData; + public readonly ulong CycleLengthHint = cycleLengthHint; + public readonly bool HasCycleLengthHint = hasCycleLengthHint; + public readonly bool RequiresCycleComputation = !hasCycleLengthHint; + } - public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096; - public static int HeuristicDivisorInterleaveBatchSize { get; set; } = 64; + [ThreadStatic] + private static HeuristicPrimeTester? _tester; - [ThreadStatic] - private static HeuristicPrimeTester? _tester; + public static HeuristicPrimeTester Exclusive => _tester ??= new(); - public static HeuristicPrimeTester Exclusive => _tester ??= new(); + public bool IsPrimeCpu(ulong n) + { + byte nMod10 = (byte)n.Mod10(); + ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); - public bool IsPrimeCpu(ulong n) - { - byte nMod10 = (byte)n.Mod10(); - ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); + if (maxDivisorSquare < 9UL) + { + return EvaluateWithOpenNumericFallback(n); + } - if (maxDivisorSquare < 9UL) - { - return EvaluateWithOpenNumericFallback(n); - } + bool includeGroupB = UseHeuristicGroupBTrialDivision; + return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB); + } - bool includeGroupB = UseHeuristicGroupBTrialDivision; - return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB); - } + public bool IsPrimeGpu(ulong n) + { + byte nMod10 = (byte)n.Mod10(); + ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); + return IsPrimeGpu(n, maxDivisorSquare, nMod10); + } - public bool IsPrimeGpu(ulong n) - { - byte nMod10 = (byte)n.Mod10(); - ulong maxDivisorSquare = ComputeHeuristicDivisorSquareLimit(n); - return IsPrimeGpu(n, maxDivisorSquare, nMod10); - } + public bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10) + { + return HeuristicIsPrimeGpuCore(n, maxDivisorSquare, nMod10); + } - public bool IsPrimeGpu(ulong n, ulong maxDivisorSquare, byte nMod10) - { - return HeuristicIsPrimeGpuCore(n, maxDivisorSquare, nMod10); - } + private bool HeuristicIsPrimeGpuCore(ulong n, ulong maxDivisorSquare, byte nMod10) + { + // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? + if (maxDivisorSquare < 9UL) + { + return EvaluateWithOpenNumericFallback(n); + } - private bool HeuristicIsPrimeGpuCore(ulong n, ulong maxDivisorSquare, byte nMod10) - { - // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? - if (maxDivisorSquare < 9UL) - { - return EvaluateWithOpenNumericFallback(n); - } + // if (!UseHeuristicGroupBTrialDivision) + // { + // return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB: false); + // } - // if (!UseHeuristicGroupBTrialDivision) - // { - // return HeuristicTrialDivisionCpu(n, maxDivisorSquare, nMod10, includeGroupB: false); - // } + bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10); + return !compositeDetected; + } - bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, maxDivisorSquare, nMod10); - return !compositeDetected; - } + private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10, bool includeGroupB) + { + ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; + ReadOnlySpan groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; + int interleaveBatchSize = HeuristicCombinedPrimeTester.GetBatchSize(n); + int groupAIndex = 0; + if (!includeGroupB) + { + int groupALength = groupADivisors.Length; + while (groupAIndex < groupALength) + { + int processed = 0; + while (processed < interleaveBatchSize && groupAIndex < groupALength) + { + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) + { + groupAIndex = groupALength; + break; + } - private static bool HeuristicTrialDivisionCpu(ulong n, ulong maxDivisorSquare, byte nMod10, bool includeGroupB) - { - ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; - ReadOnlySpan groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; - int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); - int groupAIndex = 0; + ulong divisor = (ulong)groupADivisors[groupAIndex]; + if (n % divisor == 0UL) + { + return false; + } - if (!includeGroupB) - { - int groupALength = groupADivisors.Length; - while (groupAIndex < groupALength) - { - int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupALength) - { - ulong divisorSquare = groupADivisorSquares[groupAIndex]; - if (divisorSquare > maxDivisorSquare) - { - groupAIndex = groupALength; - break; - } + groupAIndex++; + processed++; + } - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (n % divisor == 0UL) - { - return false; - } + if (groupAIndex >= groupALength) + { + break; + } - groupAIndex++; - processed++; - } + if (groupADivisorSquares[groupAIndex] > maxDivisorSquare) + { + groupAIndex = groupALength; + break; + } + } - if (groupAIndex >= groupALength) - { - break; - } + return EvaluateWithOpenNumericFallback(n); + } - if (groupADivisorSquares[groupAIndex] > maxDivisorSquare) - { - groupAIndex = groupALength; - break; - } - } + ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); + if (endingOrder.IsEmpty) + { + return EvaluateWithOpenNumericFallback(n); + } - return EvaluateWithOpenNumericFallback(n); - } + Span indices = endingOrder.Length <= 8 + ? stackalloc int[endingOrder.Length] + : new int[endingOrder.Length]; - ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); - if (endingOrder.IsEmpty) - { - return EvaluateWithOpenNumericFallback(n); - } + bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare; + bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); - Span indices = endingOrder.Length <= 8 - ? stackalloc int[endingOrder.Length] - : new int[endingOrder.Length]; + while (groupAHasMore || groupBHasMore) + { + if (groupAHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) + { + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) + { + groupAIndex = groupADivisors.Length; + groupAHasMore = false; + break; + } - bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare; - bool groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); + ulong divisor = (ulong)groupADivisors[groupAIndex]; + if (n % divisor == 0UL) + { + return false; + } - while (groupAHasMore || groupBHasMore) - { - if (groupAHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) - { - ulong divisorSquare = groupADivisorSquares[groupAIndex]; - if (divisorSquare > maxDivisorSquare) - { - groupAIndex = groupADivisors.Length; - groupAHasMore = false; - break; - } + groupAIndex++; + processed++; + } - ulong divisor = (ulong)groupADivisors[groupAIndex]; - if (n % divisor == 0UL) - { - return false; - } + if (groupAIndex >= groupADivisors.Length) + { + groupAHasMore = false; + } + else + { + groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare; + if (!groupAHasMore) + { + groupAIndex = groupADivisors.Length; + } + } + } - groupAIndex++; - processed++; - } - - if (groupAIndex >= groupADivisors.Length) - { - groupAHasMore = false; - } - else - { - groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare; - if (!groupAHasMore) - { - groupAIndex = groupADivisors.Length; - } - } - } - - if (groupBHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize) - { + if (groupBHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize) + { if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor)) - { - groupBHasMore = false; - break; - } + { + groupBHasMore = false; + break; + } - if (n % divisor == 0UL) - { - return false; - } + if (n % divisor == 0UL) + { + return false; + } - processed++; - } + processed++; + } - if (groupBHasMore) - { - groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); - } - } - } + if (groupBHasMore) + { + groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); + } + } + } - return true; - } + return true; + } + + private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10) + { + int batchCapacity = HeuristicCombinedPrimeTester.GetBatchSize(n); + int interleaveBatchSize = batchCapacity; + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity); + var accelerator = gpu.Accelerator; + var kernel = gpu.HeuristicTrialDivisionKernel; + bool compositeDetected = false; + int count = 0; + var inputView = gpu.Input.View; + var flagView = gpu.HeuristicFlag.View; + ulong[] divisorArray = gpu.DivisorArray; + bool ProcessBatch(int length) + { + if (length <= 0) + { + return false; + } - - private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong maxDivisorSquare, byte nMod10) - { - int batchCapacity = HeuristicGpuDivisorBatchSize; - int interleaveBatchSize = HeuristicDivisorInterleaveBatchSize; - var divisorPool = ThreadStaticPools.UlongPool; - - ulong[] divisorArray = divisorPool.Rent(batchCapacity); - - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(batchCapacity); - var accelerator = gpu.Accelerator; - var kernel = gpu.HeuristicTrialDivisionKernel; - - bool compositeDetected = false; - int count = 0; - - var inputView = gpu.Input.View; - var flagView = gpu.HeuristicFlag.View; - - bool ProcessBatch(int length) - { - if (length <= 0) - { - return false; - } - - inputView.CopyFromCPU(ref divisorArray[0], length); - int compositeFlag = 0; - flagView.CopyFromCPU(ref compositeFlag, 1); + inputView.CopyFromCPU(ref divisorArray[0], length); + int compositeFlag = 0; + flagView.CopyFromCPU(ref compositeFlag, 1); kernel(length, inputView, n, flagView); - accelerator.Synchronize(); - flagView.CopyToCPU(ref compositeFlag, 1); - return compositeFlag != 0; - } + accelerator.Synchronize(); + flagView.CopyToCPU(ref compositeFlag, 1); + return compositeFlag != 0; + } - ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; - ReadOnlySpan groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; - int groupAIndex = 0; - ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); - Span indices = endingOrder.Length <= 8 - ? stackalloc int[endingOrder.Length] - : new int[endingOrder.Length]; + ReadOnlySpan groupADivisors = HeuristicPrimeSieves.GroupADivisors; + ReadOnlySpan groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; + int groupAIndex = 0; + ReadOnlySpan endingOrder = GetGroupBEndingOrder(nMod10); + Span indices = endingOrder.Length <= 8 + ? stackalloc int[endingOrder.Length] + : new int[endingOrder.Length]; - bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare; - bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); + bool groupAHasMore = groupAIndex < groupADivisors.Length && groupADivisorSquares[groupAIndex] <= maxDivisorSquare; + bool groupBHasMore = !endingOrder.IsEmpty && HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); - while (!compositeDetected && (groupAHasMore || groupBHasMore)) - { - if (groupAHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) - { - ulong divisorSquare = groupADivisorSquares[groupAIndex]; - if (divisorSquare > maxDivisorSquare) - { - groupAIndex = groupADivisors.Length; - groupAHasMore = false; - break; - } + while (!compositeDetected && (groupAHasMore || groupBHasMore)) + { + if (groupAHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize && groupAIndex < groupADivisors.Length) + { + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) + { + groupAIndex = groupADivisors.Length; + groupAHasMore = false; + break; + } - divisorArray[count] = (ulong)groupADivisors[groupAIndex]; - count++; + divisorArray[count] = (ulong)groupADivisors[groupAIndex]; + count++; - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - break; - } + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + break; + } - count = 0; - } + count = 0; + } - groupAIndex++; - processed++; - } + groupAIndex++; + processed++; + } - if (compositeDetected) - { - break; - } + if (compositeDetected) + { + break; + } - if (groupAIndex >= groupADivisors.Length) - { - groupAHasMore = false; - } - else - { - groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare; - if (!groupAHasMore) - { - groupAIndex = groupADivisors.Length; - } - } - } + if (groupAIndex >= groupADivisors.Length) + { + groupAHasMore = false; + } + else + { + groupAHasMore = groupADivisorSquares[groupAIndex] <= maxDivisorSquare; + if (!groupAHasMore) + { + groupAIndex = groupADivisors.Length; + } + } + } - if (!compositeDetected && groupBHasMore) - { - int processed = 0; - while (processed < interleaveBatchSize) - { + if (!compositeDetected && groupBHasMore) + { + int processed = 0; + while (processed < interleaveBatchSize) + { if (!TrySelectNextGroupBDivisor(endingOrder, indices, maxDivisorSquare, out ulong divisor)) - { - groupBHasMore = false; - break; - } + { + groupBHasMore = false; + break; + } - divisorArray[count] = divisor; - count++; + divisorArray[count] = divisor; + count++; - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - break; - } + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + break; + } - count = 0; - } + count = 0; + } - processed++; - } + processed++; + } - if (!compositeDetected && groupBHasMore) - { - groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); - } - } - } + if (!compositeDetected && groupBHasMore) + { + groupBHasMore = HasGroupBCandidates(endingOrder, indices, maxDivisorSquare); + } + } + } - if (!compositeDetected && count > 0 && ProcessBatch(count)) - { - compositeDetected = true; - } + if (!compositeDetected && count > 0 && ProcessBatch(count)) + { + compositeDetected = true; + } - divisorPool.Return(divisorArray); - gpu.Dispose(); - limiter.Dispose(); - return compositeDetected; - } + gpu.Dispose(); + limiter.Dispose(); + return compositeDetected; + } + private static bool EvaluateWithOpenNumericFallback(ulong n) + { + return Prime.Numbers.IsPrime(n); + } + internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n) + { + return n; + } + private static ReadOnlySpan GetGroupBEndingOrder(byte nMod10) => nMod10 switch + { + 1 => GroupBEndingOrderMod1, + 3 => GroupBEndingOrderMod3, + 7 => GroupBEndingOrderMod7, + 9 => GroupBEndingOrderMod9, + _ => ReadOnlySpan.Empty, + }; - private static bool EvaluateWithOpenNumericFallback(ulong n) - { - return Prime.Numbers.IsPrime(n); - } - internal static ulong ComputeHeuristicDivisorSquareLimit(ulong n) - { - return n; - } - - private static ReadOnlySpan GetGroupBEndingOrder(byte nMod10) => nMod10 switch - { - 1 => GroupBEndingOrderMod1, - 3 => GroupBEndingOrderMod3, - 7 => GroupBEndingOrderMod7, - 9 => GroupBEndingOrderMod9, - _ => ReadOnlySpan.Empty, - }; - - - - [MethodImpl(MethodImplOptions.AggressiveInlining)] + [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TrySelectNextGroupBDivisor(ReadOnlySpan endingOrder, Span indices, ulong maxDivisorSquare, out ulong divisor) - { - ulong bestCandidate = ulong.MaxValue; - int bestEndingIndex = -1; + { + ulong bestCandidate = ulong.MaxValue; + int bestEndingIndex = -1; - for (int i = 0; i < endingOrder.Length; i++) - { - ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); - ReadOnlySpan squares = GetGroupBDivisorSquares(endingOrder[i]); - int index = indices[i]; - if ((uint)index >= (uint)divisors.Length) - { - continue; - } + for (int i = 0; i < endingOrder.Length; i++) + { + ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + ReadOnlySpan squares = GetGroupBDivisorSquares(endingOrder[i]); + int index = indices[i]; + if ((uint)index >= (uint)divisors.Length) + { + continue; + } - ulong candidateSquare = squares[index]; - if (candidateSquare > maxDivisorSquare) - { - indices[i] = divisors.Length; - continue; - } + ulong candidateSquare = squares[index]; + if (candidateSquare > maxDivisorSquare) + { + indices[i] = divisors.Length; + continue; + } - ulong candidate = (ulong)divisors[index]; - if (bestEndingIndex == -1 || candidate < bestCandidate) - { - bestCandidate = candidate; - bestEndingIndex = i; - } - } + ulong candidate = divisors[index]; + if (bestEndingIndex == -1 || candidate < bestCandidate) + { + bestCandidate = candidate; + bestEndingIndex = i; + } + } - if (bestEndingIndex == -1) - { - divisor = 0UL; - return false; - } + if (bestEndingIndex == -1) + { + divisor = 0UL; + return false; + } - ReadOnlySpan selectedDivisors = GetGroupBDivisors(endingOrder[bestEndingIndex]); - divisor = (ulong)selectedDivisors[indices[bestEndingIndex]]; - indices[bestEndingIndex]++; - return true; - } + ReadOnlySpan selectedDivisors = GetGroupBDivisors(endingOrder[bestEndingIndex]); + divisor = selectedDivisors[indices[bestEndingIndex]]; + indices[bestEndingIndex]++; + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static bool HasGroupBCandidates(ReadOnlySpan endingOrder, Span indices, ulong maxDivisorSquare) + { + for (int i = 0; i < endingOrder.Length; i++) + { + ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + ReadOnlySpan squares = GetGroupBDivisorSquares(endingOrder[i]); + int index = indices[i]; + if ((uint)index >= (uint)divisors.Length) + { + continue; + } + + if (squares[index] > maxDivisorSquare) + { + indices[i] = divisors.Length; + continue; + } + + return true; + } + + return false; + } + private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch + { + 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, + 3 => ReadOnlySpan.Empty, + 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, + 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, + _ => ReadOnlySpan.Empty, + }; - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool HasGroupBCandidates(ReadOnlySpan endingOrder, Span indices, ulong maxDivisorSquare) - { - for (int i = 0; i < endingOrder.Length; i++) - { - ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); - ReadOnlySpan squares = GetGroupBDivisorSquares(endingOrder[i]); - int index = indices[i]; - if ((uint)index >= (uint)divisors.Length) - { - continue; - } + private static ReadOnlySpan GetGroupBDivisorSquares(byte ending) => ending switch + { + 1 => DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree, + 3 => ReadOnlySpan.Empty, + 7 => DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree, + 9 => DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree, + _ => ReadOnlySpan.Empty, + }; - if (squares[index] > maxDivisorSquare) - { - indices[i] = divisors.Length; - continue; - } + private static int InitializeGroupBStates(byte nMod10, Span buffer) + { + ReadOnlySpan endings = GetGroupBEndingOrder(nMod10); + int count = endings.Length; - return true; - } + for (int i = 0; i < count; i++) + { + buffer[i] = new HeuristicGroupBSequenceState(endings[i], (byte)i); + } - return false; - } + return count; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) + { + return new HeuristicDivisorEnumerator(maxDivisorSquare, 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 maxDivisorSquare; + private readonly ReadOnlySpan groupADivisors; + private readonly ReadOnlySpan groupADivisorSquares; + private int groupAIndex; + private Span groupBStates; + + public HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) + { + this.maxDivisorSquare = maxDivisorSquare; + groupADivisors = HeuristicPrimeSieves.GroupADivisors; + groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; + groupAIndex = 0; + + int count = InitializeGroupBStates(nMod10, groupBBuffer); + groupBStates = count == 0 ? Span.Empty : groupBBuffer[..count]; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (groupAIndex < groupADivisors.Length) + { + ulong divisorSquare = groupADivisorSquares[groupAIndex]; + if (divisorSquare > maxDivisorSquare) + { + groupAIndex = groupADivisors.Length; + break; + } + + int currentIndex = groupAIndex; + ulong divisor = (ulong)groupADivisors[currentIndex]; + groupAIndex++; + + HeuristicDivisorGroup group = currentIndex < GroupAConstantCount + ? HeuristicDivisorGroup.GroupAConstant + : HeuristicDivisorGroup.GroupAWheel; + + candidate = new HeuristicDivisorCandidate( + divisor, + group, + (byte)(divisor % 10UL), + 0, + (ushort)(divisor % Wheel210)); + return true; + } + + ulong bestCandidate = ulong.MaxValue; + int bestIndex = -1; + + for (int i = 0; i < groupBStates.Length; i++) + { + ref HeuristicGroupBSequenceState state = ref groupBStates[i]; + if (!state.TryPeek(maxDivisorSquare, out ulong candidateValue)) + { + continue; + } + + if (bestIndex == -1 || candidateValue < bestCandidate || + (candidateValue == bestCandidate && state.PriorityIndex < groupBStates[bestIndex].PriorityIndex)) + { + bestCandidate = candidateValue; + bestIndex = i; + } + } + + if (bestIndex == -1) + { + candidate = default; + return false; + } + + ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex]; + bestState.Advance(); + candidate = new HeuristicDivisorCandidate( + bestCandidate, + HeuristicDivisorGroup.GroupB, + bestState.Ending, + bestState.PriorityIndex, + (ushort)(bestCandidate % Wheel210)); + return true; + } + } - private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch - { - 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, - 3 => ReadOnlySpan.Empty, - 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, - 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, - _ => ReadOnlySpan.Empty, - }; + internal struct HeuristicGroupBSequenceState + { + public byte Ending; + public byte PriorityIndex; + public int Index; - private static ReadOnlySpan GetGroupBDivisorSquares(byte ending) => ending switch - { - 1 => DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree, - 3 => ReadOnlySpan.Empty, - 7 => DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree, - 9 => DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree, - _ => ReadOnlySpan.Empty, - }; + public HeuristicGroupBSequenceState(byte ending, byte priorityIndex) + { + Ending = ending; + PriorityIndex = priorityIndex; + Index = 0; + } - private static int InitializeGroupBStates(byte nMod10, Span buffer) - { - ReadOnlySpan endings = GetGroupBEndingOrder(nMod10); - int count = endings.Length; + public bool TryPeek(ulong maxDivisorSquare, out ulong divisor) + { + ReadOnlySpan divisors = GetGroupBDivisors(Ending); + ReadOnlySpan squares = GetGroupBDivisorSquares(Ending); + int index = Index; + if ((uint)index >= (uint)divisors.Length) + { + divisor = 0UL; + return false; + } - for (int i = 0; i < count; i++) - { - buffer[i] = new HeuristicGroupBSequenceState(endings[i], (byte)i); - } + if (squares[index] > maxDivisorSquare) + { + Index = divisors.Length; + divisor = 0UL; + return false; + } - return count; - } + divisor = (ulong)divisors[index]; + return true; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) - { - return new HeuristicDivisorEnumerator(maxDivisorSquare, 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; - } + public void Advance() + { + int next = Index + 1; + ReadOnlySpan divisors = GetGroupBDivisors(Ending); + if ((uint)next >= (uint)divisors.Length) + { + Index = divisors.Length; + } + else + { + Index = next; + } + } + } -internal ref struct HeuristicDivisorEnumerator -{ - private readonly ulong maxDivisorSquare; - private readonly ReadOnlySpan groupADivisors; - private readonly ReadOnlySpan groupADivisorSquares; - private int groupAIndex; - private Span groupBStates; + 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 HeuristicDivisorEnumerator(ulong maxDivisorSquare, byte nMod10, Span groupBBuffer) - { - this.maxDivisorSquare = maxDivisorSquare; - groupADivisors = HeuristicPrimeSieves.GroupADivisors; - groupADivisorSquares = HeuristicPrimeSieves.GroupADivisorSquares; - groupAIndex = 0; + public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + var stepLocal = new GpuUInt128(exponent); + stepLocal.ShiftLeft(1); + step = stepLocal; - int count = InitializeGroupBStates(nMod10, groupBBuffer); - groupBStates = count == 0 ? Span.Empty : groupBBuffer[..count]; - } + var limitLocal = new GpuUInt128(maxDivisor); + limit = limitLocal; - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - while (groupAIndex < groupADivisors.Length) - { - ulong divisorSquare = groupADivisorSquares[groupAIndex]; - if (divisorSquare > maxDivisorSquare) - { - groupAIndex = groupADivisors.Length; - break; - } + var firstDivisor = stepLocal; + firstDivisor.Add(1UL); - int currentIndex = groupAIndex; - ulong divisor = (ulong)groupADivisors[currentIndex]; - groupAIndex++; + bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0; - HeuristicDivisorGroup group = currentIndex < GroupAConstantCount - ? HeuristicDivisorGroup.GroupAConstant - : HeuristicDivisorGroup.GroupAWheel; + current = hasCandidates ? firstDivisor : GpuUInt128.Zero; + residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default; + active = hasCandidates; + processedCount = 0UL; + lastDivisor = 0UL; + } - candidate = new HeuristicDivisorCandidate( - divisor, - group, - (byte)(divisor % 10UL), - 0, - (ushort)(divisor % Wheel210)); - return true; - } + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (active) + { + ulong value = current.Low; + processedCount++; + lastDivisor = value; - ulong bestCandidate = ulong.MaxValue; - int bestIndex = -1; + bool admissible = residueStepper.IsAdmissible(); - for (int i = 0; i < groupBStates.Length; i++) - { - ref HeuristicGroupBSequenceState state = ref groupBStates[i]; - if (!state.TryPeek(maxDivisorSquare, out ulong candidateValue)) - { - continue; - } + Advance(); - if (bestIndex == -1 || candidateValue < bestCandidate || - (candidateValue == bestCandidate && state.PriorityIndex < groupBStates[bestIndex].PriorityIndex)) - { - bestCandidate = candidateValue; - bestIndex = i; - } - } + if (admissible) + { + candidate = CreateCandidate(value); + return true; + } + } - if (bestIndex == -1) - { - candidate = default; - return false; - } + candidate = default; + return false; + } - ref HeuristicGroupBSequenceState bestState = ref groupBStates[bestIndex]; - bestState.Advance(); - candidate = new HeuristicDivisorCandidate( - bestCandidate, - HeuristicDivisorGroup.GroupB, - bestState.Ending, - bestState.PriorityIndex, - (ushort)(bestCandidate % Wheel210)); - return true; - } -} - - - -internal struct HeuristicGroupBSequenceState -{ - public byte Ending; - public byte PriorityIndex; - public int Index; - - public HeuristicGroupBSequenceState(byte ending, byte priorityIndex) - { - Ending = ending; - PriorityIndex = priorityIndex; - Index = 0; - } - - public bool TryPeek(ulong maxDivisorSquare, out ulong divisor) - { - ReadOnlySpan divisors = GetGroupBDivisors(Ending); - ReadOnlySpan squares = GetGroupBDivisorSquares(Ending); - int index = Index; - if ((uint)index >= (uint)divisors.Length) - { - divisor = 0UL; - return false; - } - - if (squares[index] > maxDivisorSquare) - { - Index = divisors.Length; - divisor = 0UL; - return false; - } - - divisor = (ulong)divisors[index]; - return true; - } - - public void Advance() - { - int next = Index + 1; - ReadOnlySpan divisors = GetGroupBDivisors(Ending); - if ((uint)next >= (uint)divisors.Length) - { - Index = divisors.Length; - } - else - { - Index = next; - } - } -} - - - internal ref struct MersenneHeuristicDivisorEnumerator - { - 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); - } - } + 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/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 52bb5bf..0b46256 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -273,6 +273,7 @@ public sealed class PrimeTester public MemoryBuffer1D Output; public MemoryBuffer1D HeuristicFlag; public int BufferCapacity; + public ulong[] DivisorArray; internal PrimeTesterGpuContextLease(int minBufferCapacity) { @@ -321,16 +322,18 @@ public sealed class PrimeTester DevicePrimesPow2LastNine = Accelerator.Allocate1D(primesPow2LastNine.Length); DevicePrimesPow2LastNine.View.CopyFromCPU(primesPow2LastNine); - BufferCapacity = minBufferCapacity < 1 ? 1 : minBufferCapacity; - Input = Accelerator.Allocate1D(BufferCapacity); - Output = Accelerator.Allocate1D(BufferCapacity); - HeuristicFlag = Accelerator.Allocate1D(BufferCapacity); + BufferCapacity = minBufferCapacity; + Input = Accelerator.Allocate1D(minBufferCapacity); + Output = Accelerator.Allocate1D(minBufferCapacity); + HeuristicFlag = Accelerator.Allocate1D(minBufferCapacity); + DivisorArray = ThreadStaticPools.UlongPool.Rent(minBufferCapacity); + _returned = false; } public void EnsureCapacity(int minCapacity) { - int required = minCapacity < 1 ? 1 : minCapacity; + int required = minCapacity; if (required <= BufferCapacity) { return; @@ -339,9 +342,11 @@ public sealed class PrimeTester Input.Dispose(); Output.Dispose(); HeuristicFlag.Dispose(); + ThreadStaticPools.UlongPool.Return(DivisorArray); Input = Accelerator.Allocate1D(required); Output = Accelerator.Allocate1D(required); HeuristicFlag = Accelerator.Allocate1D(required); + DivisorArray = ThreadStaticPools.UlongPool.Rent(required); BufferCapacity = required; }