From c693a057083b860d3e8db5033583e3209089a5dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 28 Oct 2025 21:56:18 +0100 Subject: [PATCH 01/16] Add reusable group B divisor tables without trailing-three primes --- .../PrimesGeneratorTests.cs | 56 +++++++++++ PerfectNumbers.Core/DivisorGenerator.cs | 95 ++++++++++++++----- PerfectNumbers.Core/HeuristicPrimeTester.cs | 16 ++-- 3 files changed, 136 insertions(+), 31 deletions(-) diff --git a/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs b/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs index 95df8a2..8ccebe4 100644 --- a/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs +++ b/PerfectNumbers.Core.Tests/PrimesGeneratorTests.cs @@ -36,15 +36,30 @@ public class PrimesGeneratorTests var lastThreePow2 = DivisorGenerator.SmallPrimesPow2LastThree; var lastNine = DivisorGenerator.SmallPrimesLastNine; var lastNinePow2 = DivisorGenerator.SmallPrimesPow2LastNine; + var lastOneWithoutThree = DivisorGenerator.SmallPrimesLastOneWithoutLastThree; + var lastOneWithoutThreePow2 = DivisorGenerator.SmallPrimesPow2LastOneWithoutLastThree; + var lastSevenWithoutThree = DivisorGenerator.SmallPrimesLastSevenWithoutLastThree; + var lastSevenWithoutThreePow2 = DivisorGenerator.SmallPrimesPow2LastSevenWithoutLastThree; + var lastNineWithoutThree = DivisorGenerator.SmallPrimesLastNineWithoutLastThree; + var lastNineWithoutThreePow2 = DivisorGenerator.SmallPrimesPow2LastNineWithoutLastThree; lastOne.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); lastSeven.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); lastThree.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); lastNine.Length.Should().Be((int)PerfectNumberConstants.PrimesLimit); + lastOneWithoutThree.Length.Should().BeGreaterThan(0); + lastSevenWithoutThree.Length.Should().BeGreaterThan(0); + lastNineWithoutThree.Length.Should().BeGreaterThan(0); + lastOneWithoutThree.Length.Should().BeLessThan(lastOne.Length); + lastSevenWithoutThree.Length.Should().BeLessThan(lastSeven.Length); + lastNineWithoutThree.Length.Should().BeLessThan(lastNine.Length); lastOnePow2.Length.Should().Be(lastOne.Length); lastSevenPow2.Length.Should().Be(lastSeven.Length); lastThreePow2.Length.Should().Be(lastThree.Length); lastNinePow2.Length.Should().Be(lastNine.Length); + lastOneWithoutThreePow2.Length.Should().Be(lastOneWithoutThree.Length); + lastSevenWithoutThreePow2.Length.Should().Be(lastSevenWithoutThree.Length); + lastNineWithoutThreePow2.Length.Should().Be(lastNineWithoutThree.Length); lastOne[..5].Should().Equal([3U, 7U, 11U, 13U, 19U]); lastSeven[..5].Should().Equal([3U, 7U, 11U, 13U, 17U]); @@ -78,6 +93,47 @@ public class PrimesGeneratorTests lastThreePow2[i].Should().Be(prime * (ulong)prime); } + + lastOneWithoutThree.Should().NotBeEmpty(); + lastSevenWithoutThree.Should().NotBeEmpty(); + lastNineWithoutThree.Should().NotBeEmpty(); + + for (int i = 0; i < lastOneWithoutThree.Length; i++) + { + uint prime = lastOneWithoutThree[i]; + (prime % 10U != 3U).Should().BeTrue(); + if (i > 0) + { + prime.Should().BeGreaterThan(lastOneWithoutThree[i - 1]); + } + + lastOneWithoutThreePow2[i].Should().Be(prime * (ulong)prime); + } + + for (int i = 0; i < lastSevenWithoutThree.Length; i++) + { + uint prime = lastSevenWithoutThree[i]; + (prime % 10U != 3U).Should().BeTrue(); + if (i > 0) + { + prime.Should().BeGreaterThan(lastSevenWithoutThree[i - 1]); + } + + lastSevenWithoutThreePow2[i].Should().Be(prime * (ulong)prime); + } + + for (int i = 0; i < lastNineWithoutThree.Length; i++) + { + uint prime = lastNineWithoutThree[i]; + (prime % 10U != 3U).Should().BeTrue(); + if (i > 0) + { + prime.Should().BeGreaterThan(lastNineWithoutThree[i - 1]); + } + + lastNineWithoutThreePow2[i].Should().Be(prime * (ulong)prime); + } + for (int i = 0; i < lastNine.Length; i++) { uint prime = lastNine[i]; diff --git a/PerfectNumbers.Core/DivisorGenerator.cs b/PerfectNumbers.Core/DivisorGenerator.cs index 20b8c81..4cce836 100644 --- a/PerfectNumbers.Core/DivisorGenerator.cs +++ b/PerfectNumbers.Core/DivisorGenerator.cs @@ -16,6 +16,12 @@ internal static class DivisorGenerator public static readonly ulong[] SmallPrimesPow2LastThree; public static readonly uint[] SmallPrimesLastNine; public static readonly ulong[] SmallPrimesPow2LastNine; + public static readonly uint[] SmallPrimesLastOneWithoutLastThree; + public static readonly ulong[] SmallPrimesPow2LastOneWithoutLastThree; + public static readonly uint[] SmallPrimesLastSevenWithoutLastThree; + public static readonly ulong[] SmallPrimesPow2LastSevenWithoutLastThree; + public static readonly uint[] SmallPrimesLastNineWithoutLastThree; + public static readonly ulong[] SmallPrimesPow2LastNineWithoutLastThree; private const ushort DecimalMaskWhenLastIsSeven = (1 << 3) | (1 << 7) | (1 << 9); private const ushort DecimalMaskWhenLastIsThree = (1 << 1) | (1 << 3) | (1 << 7) | (1 << 9); @@ -34,7 +40,13 @@ internal static class DivisorGenerator out SmallPrimesLastThree, out SmallPrimesPow2LastThree, out SmallPrimesLastNine, - out SmallPrimesPow2LastNine); + out SmallPrimesPow2LastNine, + out SmallPrimesLastOneWithoutLastThree, + out SmallPrimesPow2LastOneWithoutLastThree, + out SmallPrimesLastSevenWithoutLastThree, + out SmallPrimesPow2LastSevenWithoutLastThree, + out SmallPrimesLastNineWithoutLastThree, + out SmallPrimesPow2LastNineWithoutLastThree); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -59,7 +71,13 @@ internal static class DivisorGenerator out uint[] lastThree, out ulong[] lastThreePow2, out uint[] lastNine, - out ulong[] lastNinePow2) + out ulong[] lastNinePow2, + out uint[] lastOneWithoutThree, + out ulong[] lastOneWithoutThreePow2, + out uint[] lastSevenWithoutThree, + out ulong[] lastSevenWithoutThreePow2, + out uint[] lastNineWithoutThree, + out ulong[] lastNineWithoutThreePow2) { uint target = PerfectNumberConstants.PrimesLimit; int targetInt = checked((int)target); @@ -75,6 +93,13 @@ internal static class DivisorGenerator lastNinePow2 = new ulong[targetInt]; List primes = new(targetInt * 4 / 3 + 16); + List lastOneWithoutThreeList = new(targetInt); + List lastOneWithoutThreePow2List = new(targetInt); + List lastSevenWithoutThreeList = new(targetInt); + List lastSevenWithoutThreePow2List = new(targetInt); + List lastNineWithoutThreeList = new(targetInt); + List lastNineWithoutThreePow2List = new(targetInt); + uint candidate = 2U; int allCount = 0; int lastOneCount = 0; @@ -82,24 +107,22 @@ internal static class DivisorGenerator int lastThreeCount = 0; int lastNineCount = 0; - while (allCount < targetInt || lastOneCount < targetInt || lastSevenCount < targetInt || lastThreeCount < targetInt || lastNineCount < targetInt) + while (allCount < targetInt || + lastOneCount < targetInt || + lastSevenCount < targetInt || + lastThreeCount < targetInt || + lastNineCount < targetInt) { bool isPrime = true; int primesCount = primes.Count; for (int i = 0; i < primesCount; i++) { uint p = primes[i]; - // TODO: Reuse a single squared local for `p` within this loop so we avoid recalculating - // `p * p` repeatedly without storing every squared prime globally (the cache must remain - // limited to the small-divisor tables to honor the memory constraints). if (p * p > candidate) { break; } - // TODO: Replace this trial-division `%` with the sieve-based generator that avoids - // per-candidate modulo work so building the small-prime tables stops dominating - // startup time for large scans. if (candidate % p == 0U) { isPrime = false; @@ -110,54 +133,80 @@ internal static class DivisorGenerator if (isPrime) { primes.Add(candidate); + ulong candidateSquared = checked((ulong)candidate * candidate); if (allCount < targetInt) { - // TODO: Reuse a single cached `candidateSquared` per iteration so we avoid issuing - // three separate 64-bit multiplications when populating the pow2 arrays here and - // in the last-one/last-seven branches below. all[allCount] = candidate; - allPow2[allCount] = checked((ulong)candidate * candidate); + allPow2[allCount] = candidateSquared; allCount++; } if (candidate != 2U) { LastDigit lastDigit = GetLastDigit(candidate); + bool allowedForLastOne = IsAllowedForLastOne(candidate, lastDigit); + bool allowedForLastSeven = IsAllowedForLastSeven(candidate, lastDigit); + bool allowedForLastThree = IsAllowedForLastThree(candidate, lastDigit); + bool allowedForLastNine = IsAllowedForLastNine(candidate, lastDigit); + bool lastDigitIsThree = lastDigit == LastDigit.Three; - if (lastOneCount < targetInt && IsAllowedForLastOne(candidate, lastDigit)) + if (allowedForLastOne && lastOneCount < targetInt) { lastOne[lastOneCount] = candidate; - lastOnePow2[lastOneCount] = checked((ulong)candidate * candidate); + lastOnePow2[lastOneCount] = candidateSquared; lastOneCount++; } - if (lastSevenCount < targetInt && IsAllowedForLastSeven(candidate, lastDigit)) + if (allowedForLastSeven && lastSevenCount < targetInt) { lastSeven[lastSevenCount] = candidate; - lastSevenPow2[lastSevenCount] = checked((ulong)candidate * candidate); + lastSevenPow2[lastSevenCount] = candidateSquared; lastSevenCount++; } - if (lastThreeCount < targetInt && IsAllowedForLastThree(candidate, lastDigit)) + if (allowedForLastThree && lastThreeCount < targetInt) { lastThree[lastThreeCount] = candidate; - lastThreePow2[lastThreeCount] = checked((ulong)candidate * candidate); + lastThreePow2[lastThreeCount] = candidateSquared; lastThreeCount++; } - if (lastNineCount < targetInt && IsAllowedForLastNine(candidate, lastDigit)) + if (allowedForLastNine && lastNineCount < targetInt) { lastNine[lastNineCount] = candidate; - lastNinePow2[lastNineCount] = checked((ulong)candidate * candidate); + lastNinePow2[lastNineCount] = candidateSquared; lastNineCount++; } + + if (!lastDigitIsThree && allowedForLastOne) + { + lastOneWithoutThreeList.Add(candidate); + lastOneWithoutThreePow2List.Add(candidateSquared); + } + + if (!lastDigitIsThree && allowedForLastSeven) + { + lastSevenWithoutThreeList.Add(candidate); + lastSevenWithoutThreePow2List.Add(candidateSquared); + } + + if (!lastDigitIsThree && allowedForLastNine) + { + lastNineWithoutThreeList.Add(candidate); + lastNineWithoutThreePow2List.Add(candidateSquared); + } } } - // TODO: Switch this increment to the Mod6 stepping table identified in the Mod6ComparisonBenchmarks - // so we skip the composite residues without relying on per-iteration `%` checks. candidate = candidate == 2U ? 3U : candidate + 2U; } + + lastOneWithoutThree = lastOneWithoutThreeList.ToArray(); + lastOneWithoutThreePow2 = lastOneWithoutThreePow2List.ToArray(); + lastSevenWithoutThree = lastSevenWithoutThreeList.ToArray(); + lastSevenWithoutThreePow2 = lastSevenWithoutThreePow2List.ToArray(); + lastNineWithoutThree = lastNineWithoutThreeList.ToArray(); + lastNineWithoutThreePow2 = lastNineWithoutThreePow2List.ToArray(); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index 2cbf881..c3983a1 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -470,7 +470,7 @@ Cleanup: for (int i = 0; i < endingOrder.Length; i++) { - ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); int index = indices[i]; if ((uint)index >= (uint)divisors.Length) { @@ -507,7 +507,7 @@ Cleanup: { for (int i = 0; i < endingOrder.Length; i++) { - ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); + ReadOnlySpan divisors = GetGroupBDivisors(endingOrder[i]); int index = indices[i]; if ((uint)index >= (uint)divisors.Length) { @@ -527,13 +527,13 @@ Cleanup: return false; } - private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch + private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch { - 1 => HeuristicPrimeSieves.GroupBDivisorsEnding1, - 3 => HeuristicPrimeSieves.GroupBDivisorsEnding3, - 7 => HeuristicPrimeSieves.GroupBDivisorsEnding7, - 9 => HeuristicPrimeSieves.GroupBDivisorsEnding9, - _ => ReadOnlySpan.Empty, + 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, + 3 => ReadOnlySpan.Empty, + 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, + 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, + _ => ReadOnlySpan.Empty, }; private static ReadOnlySpan GetWheel210ResiduesForEnding(byte ending) => ending switch From 41bff698e8c65e24672cf63840f8a45a00acf845 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 28 Oct 2025 22:40:51 +0100 Subject: [PATCH 02/16] Add combined heuristic prime tester with shared divisor arrays --- .../HeuristicCombinedPrimeTester.cs | 761 ++++++++++++++++++ 1 file changed, 761 insertions(+) create mode 100644 PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs new file mode 100644 index 0000000..666c2ed --- /dev/null +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -0,0 +1,761 @@ +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 HeuristicCombinedPrimeTester +{ + static HeuristicCombinedPrimeTester() + { + HeuristicPrimeSieves.EnsureInitialized(); + } + + private const ulong Wheel210 = 210UL; + + private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. + + private static readonly Lazy CombinedDivisorsEnding1 = new(() => BuildCombinedDivisors(1), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding3 = new(() => BuildCombinedDivisors(3), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding7 = new(() => BuildCombinedDivisors(7), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding9 = new(() => BuildCombinedDivisors(9), LazyThreadSafetyMode.ExecutionAndPublication); + + 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 int HeuristicGpuDivisorBatchSize { get; set; } = 4_096; + + public static int HeuristicDivisorInterleaveBatchSize { get; set; } = 64; + + [ThreadStatic] + private static HeuristicCombinedPrimeTester? _tester; + + public static HeuristicCombinedPrimeTester Exclusive => _tester ??= new(); + + public bool IsPrimeCpu(ulong n, CancellationToken ct) + { + ct.ThrowIfCancellationRequested(); + + byte nMod10 = (byte)n.Mod10(); + ulong sqrtLimit = ComputeHeuristicSqrt(n); + + if (sqrtLimit < 3UL) + { + return EvaluateWithOpenNumericFallback(n); + } + + bool includeGroupB = UseHeuristicGroupBTrialDivision; + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB); + } + + public bool IsPrimeCpu(ulong n) + { + return IsPrimeCpu(n, CancellationToken.None); + } + + public bool IsPrimeGpu(ulong n) + { + byte nMod10 = (byte)n.Mod10(); + ulong sqrtLimit = ComputeHeuristicSqrt(n); + return IsPrimeGpu(n, sqrtLimit, nMod10); + } + + public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) + { + return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10); + } + + private bool HeuristicIsPrimeGpuCore(ulong n, ulong sqrtLimit, byte nMod10) + { + // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? + if (sqrtLimit < 3UL) + { + return EvaluateWithOpenNumericFallback(n); + } + + if (!UseHeuristicGroupBTrialDivision) + { + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); + } + + bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); + return !compositeDetected; + } + + private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB) + { + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + if (combinedDivisors.IsEmpty) + { + return EvaluateWithOpenNumericFallback(n); + } + + int targetGroupACount = CountGroupAEntriesUpTo(sqrtLimit); + int targetGroupBCount = includeGroupB ? CountGroupBEntriesUpTo(sqrtLimit, nMod10) : 0; + + if (targetGroupACount == 0 && (!includeGroupB || targetGroupBCount == 0)) + { + return EvaluateWithOpenNumericFallback(n); + } + + int processedA = 0; + int processedB = 0; + + for (int i = 0; i < combinedDivisors.Length; i++) + { + uint entry = combinedDivisors[i]; + HeuristicDivisorGroup group = ResolveGroup(entry); + + if (group == HeuristicDivisorGroup.GroupB) + { + if (!includeGroupB) + { + continue; + } + + if (processedB >= targetGroupBCount) + { + if (processedA >= targetGroupACount) + { + break; + } + + continue; + } + + processedB++; + } + else + { + if (processedA >= targetGroupACount) + { + if (!includeGroupB || processedB >= targetGroupBCount) + { + break; + } + + continue; + } + + processedA++; + } + + ulong divisor = entry; + if (divisor > sqrtLimit) + { + continue; + } + + if (n % divisor == 0UL) + { + return false; + } + + if (processedA >= targetGroupACount && (!includeGroupB || processedB >= targetGroupBCount)) + { + break; + } + } + + return EvaluateWithOpenNumericFallback(n); + } + + private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) + { + int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); + var divisorPool = ThreadStaticPools.UlongPool; + var hitPool = ThreadStaticPools.BytePool; + + ulong[] divisorArray = divisorPool.Rent(batchCapacity); + byte[] hitFlags = hitPool.Rent(batchCapacity); + + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(); + var accelerator = gpu.Accelerator; + var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); + var scratch = state.RentScratch(batchCapacity, accelerator); + + bool compositeDetected = false; + int count = 0; + + lock (gpu.ExecutionLock) + { + var input = scratch.Input; + var output = scratch.Output; + + bool ProcessBatch(int length) + { + if (length <= 0) + { + return false; + } + + input.View.CopyFromCPU(ref divisorArray[0], length); + state.HeuristicTrialDivisionKernel(length, input.View, n, output.View); + accelerator.Synchronize(); + output.View.CopyToCPU(ref hitFlags[0], length); + + for (int i = 0; i < length; i++) + { + if (hitFlags[i] == 0) + { + continue; + } + + ulong divisor = divisorArray[i]; + if (divisor > 1UL && n % divisor == 0UL) + { + return true; + } + } + + return false; + } + + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + if (combinedDivisors.IsEmpty) + { + compositeDetected = !EvaluateWithOpenNumericFallback(n); + goto Cleanup; + } + + int targetGroupACount = CountGroupAEntriesUpTo(sqrtLimit); + int targetGroupBCount = CountGroupBEntriesUpTo(sqrtLimit, nMod10); + bool includeGroupB = targetGroupBCount > 0; + + if (targetGroupACount == 0 && !includeGroupB) + { + compositeDetected = !EvaluateWithOpenNumericFallback(n); + goto Cleanup; + } + + int processedA = 0; + int processedB = 0; + + for (int i = 0; i < combinedDivisors.Length; i++) + { + uint entry = combinedDivisors[i]; + HeuristicDivisorGroup group = ResolveGroup(entry); + + if (group == HeuristicDivisorGroup.GroupB) + { + if (!includeGroupB) + { + continue; + } + + if (processedB >= targetGroupBCount) + { + if (processedA >= targetGroupACount) + { + break; + } + + continue; + } + + processedB++; + } + else + { + if (processedA >= targetGroupACount) + { + if (!includeGroupB || processedB >= targetGroupBCount) + { + break; + } + + continue; + } + + processedA++; + } + + ulong divisor = entry; + if (divisor > sqrtLimit) + { + continue; + } + + divisorArray[count] = divisor; + count++; + + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + + count = 0; + } + + if (processedA >= targetGroupACount && (!includeGroupB || processedB >= targetGroupBCount)) + { + break; + } + } + + if (!compositeDetected && count > 0) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + + count = 0; + } + } + +Cleanup: + state.ReturnScratch(scratch); + divisorPool.Return(divisorArray); + hitPool.Return(hitFlags); + gpu.Dispose(); + limiter.Dispose(); + return compositeDetected; + } + + private static bool EvaluateWithOpenNumericFallback(ulong n) + { + return Prime.Numbers.IsPrime(n); + } + + internal static ulong ComputeHeuristicSqrt(ulong n) + { + ulong sqrt = (ulong)Math.Sqrt(n); + UInt128 square = (UInt128)sqrt * sqrt; + + while (square > n) + { + sqrt--; + square = (UInt128)sqrt * sqrt; + } + + ulong next = sqrt + 1UL; + UInt128 nextSquare = (UInt128)next * next; + while (nextSquare <= n) + { + sqrt = next; + next++; + nextSquare = (UInt128)next * next; + } + + return sqrt; + } + + private static ReadOnlySpan GetCombinedDivisors(byte nMod10) => nMod10 switch + { + 1 => CombinedDivisorsEnding1.Value, + 3 => CombinedDivisorsEnding3.Value, + 7 => CombinedDivisorsEnding7.Value, + 9 => CombinedDivisorsEnding9.Value, + _ => ReadOnlySpan.Empty, + }; + + private static uint[] BuildCombinedDivisors(byte nMod10) + { + ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; + ReadOnlySpan groupB = GetGroupBDivisors(nMod10); + + var combined = new List(groupA.Length + groupB.Length); + + int groupAIndex = 0; + while (groupAIndex < groupA.Length && groupA[groupAIndex] <= 19) + { + combined.Add((uint)groupA[groupAIndex]); + groupAIndex++; + } + + int groupBIndex = GetGroupBStartIndex(groupB); + + while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) + { + int addedA = 0; + while (addedA < 2 && groupAIndex < groupA.Length) + { + combined.Add((uint)groupA[groupAIndex]); + groupAIndex++; + addedA++; + } + + if (groupBIndex < groupB.Length) + { + combined.Add(groupB[groupBIndex]); + groupBIndex++; + } + } + + return combined.ToArray(); + } + + private static int CountGroupAEntriesUpTo(ulong sqrtLimit) + { + ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; + int count = 0; + + while (count < groupA.Length && (ulong)groupA[count] <= sqrtLimit) + { + count++; + } + + return count; + } + + private static int CountGroupBEntriesUpTo(ulong sqrtLimit, byte nMod10) + { + ReadOnlySpan divisors = GetGroupBDivisors(nMod10); + if (divisors.IsEmpty) + { + return 0; + } + + int index = GetGroupBStartIndex(divisors); + int count = 0; + + for (int i = index; i < divisors.Length; i++) + { + if ((ulong)divisors[i] > sqrtLimit) + { + break; + } + + count++; + } + + return count; + } + + private static int GetGroupBStartIndex(ReadOnlySpan divisors) + { + int index = 0; + while (index < divisors.Length && IsGroupAPrefixValue(divisors[index])) + { + index++; + } + + return index; + } + + private static bool IsGroupAPrefixValue(uint value) => value is 3U or 7U or 11U or 13U; + + private static HeuristicDivisorGroup ResolveGroup(uint value) => value switch + { + 3U or 7U or 11U or 13U => HeuristicDivisorGroup.GroupAConstant, + _ when value % 10U == 3U => HeuristicDivisorGroup.GroupAWheel, + _ => HeuristicDivisorGroup.GroupB, + }; + + private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch + { + 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, + 3 => ReadOnlySpan.Empty, + 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, + 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, + _ => ReadOnlySpan.Empty, + }; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); + } + + internal struct HeuristicGroupBSequenceState + { + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) + { + ulong divisor = candidate.Value; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); + return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong ResolveHeuristicCycleLength( + ulong exponent, + in HeuristicDivisorPreparation preparation, + out bool cycleFromHint, + out bool cycleComputed, + out bool primeOrderFailed) + { + if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL) + { + cycleFromHint = true; + cycleComputed = true; + primeOrderFailed = false; + return preparation.CycleLengthHint; + } + + ulong divisor = preparation.Candidate.Value; + MontgomeryDivisorData divisorData = preparation.DivisorData; + + bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + divisor, + exponent, + divisorData, + out ulong computedCycle, + out bool primeOrderFailedLocal); + + if (trySuccess && computedCycle != 0UL) + { + cycleFromHint = false; + cycleComputed = true; + primeOrderFailed = primeOrderFailedLocal; + return computedCycle; + } + + primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL; + ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); + cycleFromHint = false; + cycleComputed = resolvedCycle != 0UL; + return resolvedCycle; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength) + { + ulong[] snapshot = HeuristicSmallCycleSnapshot; + if (divisor < (ulong)snapshot.Length) + { + cycleLength = snapshot[(int)divisor]; + if (cycleLength != 0UL) + { + return true; + } + } + + cycleLength = 0UL; + return false; + } + + internal ref struct HeuristicDivisorEnumerator + { + private readonly ulong sqrtLimit; + private readonly ReadOnlySpan combinedDivisors; + private readonly bool includeGroupB; + private readonly int targetGroupACount; + private readonly int targetGroupBCount; + private int index; + private int processedA; + private int processedB; + + public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + this.sqrtLimit = sqrtLimit; + _ = groupBBuffer; + combinedDivisors = GetCombinedDivisors(nMod10); + targetGroupACount = CountGroupAEntriesUpTo(sqrtLimit); + targetGroupBCount = CountGroupBEntriesUpTo(sqrtLimit, nMod10); + includeGroupB = targetGroupBCount > 0; + index = 0; + processedA = 0; + processedB = 0; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (index < combinedDivisors.Length) + { + uint entry = combinedDivisors[index++]; + HeuristicDivisorGroup group = ResolveGroup(entry); + + if (group == HeuristicDivisorGroup.GroupB) + { + if (!includeGroupB) + { + continue; + } + + if (processedB >= targetGroupBCount) + { + if (processedA >= targetGroupACount) + { + break; + } + + continue; + } + + processedB++; + } + else + { + if (processedA >= targetGroupACount) + { + if (!includeGroupB || processedB >= targetGroupBCount) + { + break; + } + + continue; + } + + processedA++; + } + + ulong value = entry; + if (value > sqrtLimit) + { + continue; + } + + byte ending = (byte)(value % 10UL); + ushort residue = (ushort)(value % Wheel210); + candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); + return true; + } + + candidate = default; + return false; + } + } + + internal ref struct MersenneHeuristicDivisorEnumerator + { + private readonly GpuUInt128 step; + private readonly GpuUInt128 limit; + private GpuUInt128 current; + private MersenneDivisorResidueStepper residueStepper; + private bool active; + private ulong processedCount; + private ulong lastDivisor; + + public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + var stepLocal = new GpuUInt128(exponent); + stepLocal.ShiftLeft(1); + step = stepLocal; + + var limitLocal = new GpuUInt128(maxDivisor); + limit = limitLocal; + + var firstDivisor = stepLocal; + firstDivisor.Add(1UL); + + bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0; + + current = hasCandidates ? firstDivisor : GpuUInt128.Zero; + residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default; + active = hasCandidates; + processedCount = 0UL; + lastDivisor = 0UL; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (active) + { + ulong value = current.Low; + processedCount++; + lastDivisor = value; + + bool admissible = residueStepper.IsAdmissible(); + + Advance(); + + if (admissible) + { + candidate = CreateCandidate(value); + return true; + } + } + + candidate = default; + return false; + } + + private void Advance() + { + GpuUInt128 next = current; + next.Add(step); + + if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0) + { + active = false; + return; + } + + current = next; + residueStepper.Advance(); + } + + public ulong ProcessedCount => processedCount; + + public ulong LastDivisor => lastDivisor; + + public bool Exhausted => !active; + + private static HeuristicDivisorCandidate CreateCandidate(ulong value) + { + byte ending = (byte)(value % 10UL); + HeuristicDivisorGroup group; + + if (value == 3UL || value == 7UL || value == 11UL || value == 13UL) + { + group = HeuristicDivisorGroup.GroupAConstant; + } + else if (ending == 3) + { + group = HeuristicDivisorGroup.GroupAWheel; + } + else + { + group = HeuristicDivisorGroup.GroupB; + } + + ushort residue = (ushort)(value % Wheel210); + return new HeuristicDivisorCandidate(value, group, ending, 0, residue); + } + } +} From 4c445f2a42ba985ae29fd2d458f7ca16cd18ce9b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 28 Oct 2025 22:55:17 +0100 Subject: [PATCH 03/16] Switch combined heuristic divisors to 1A-1B schedule --- .../HeuristicCombinedPrimeTester.cs | 58 ++++++++++++++----- 1 file changed, 45 insertions(+), 13 deletions(-) diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 666c2ed..1942066 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -24,13 +24,24 @@ public sealed class HeuristicCombinedPrimeTester private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. - private static readonly Lazy CombinedDivisorsEnding1 = new(() => BuildCombinedDivisors(1), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding3 = new(() => BuildCombinedDivisors(3), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding7 = new(() => BuildCombinedDivisors(7), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding9 = new(() => BuildCombinedDivisors(9), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding1TwoAOneB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding3TwoAOneB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding7TwoAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding9TwoAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + + private static readonly Lazy CombinedDivisorsEnding1OneAOneB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding3OneAOneB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding7OneAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding9OneAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + internal enum CombinedDivisorPattern : byte + { + TwoAOneB = 0, + OneAOneB = 1, + } + internal enum HeuristicDivisorGroup : byte { None = 0, @@ -385,16 +396,36 @@ Cleanup: return sqrt; } - private static ReadOnlySpan GetCombinedDivisors(byte nMod10) => nMod10 switch + private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { - 1 => CombinedDivisorsEnding1.Value, - 3 => CombinedDivisorsEnding3.Value, - 7 => CombinedDivisorsEnding7.Value, - 9 => CombinedDivisorsEnding9.Value, - _ => ReadOnlySpan.Empty, - }; + return GetCombinedDivisors(nMod10, CombinedDivisorPattern.OneAOneB); + } - private static uint[] BuildCombinedDivisors(byte nMod10) + internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) + { + return pattern switch + { + CombinedDivisorPattern.TwoAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1TwoAOneB.Value, + 3 => CombinedDivisorsEnding3TwoAOneB.Value, + 7 => CombinedDivisorsEnding7TwoAOneB.Value, + 9 => CombinedDivisorsEnding9TwoAOneB.Value, + _ => ReadOnlySpan.Empty, + }, + CombinedDivisorPattern.OneAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1OneAOneB.Value, + 3 => CombinedDivisorsEnding3OneAOneB.Value, + 7 => CombinedDivisorsEnding7OneAOneB.Value, + 9 => CombinedDivisorsEnding9OneAOneB.Value, + _ => ReadOnlySpan.Empty, + }, + _ => ReadOnlySpan.Empty, + }; + } + + private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) { ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; ReadOnlySpan groupB = GetGroupBDivisors(nMod10); @@ -409,11 +440,12 @@ Cleanup: } int groupBIndex = GetGroupBStartIndex(groupB); + int groupAInterleaveCount = pattern == CombinedDivisorPattern.OneAOneB ? 1 : 2; while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) { int addedA = 0; - while (addedA < 2 && groupAIndex < groupA.Length) + while (addedA < groupAInterleaveCount && groupAIndex < groupA.Length) { combined.Add((uint)groupA[groupAIndex]); groupAIndex++; From 465b40bf353de7e1850f60e9820db94aad99fd53 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 28 Oct 2025 23:09:41 +0100 Subject: [PATCH 04/16] Add three-A-one-B schedule to combined tester --- .../HeuristicCombinedPrimeTesterTests.cs | 93 +++++++++++++++++++ .../HeuristicCombinedPrimeTester.cs | 24 ++++- 2 files changed, 115 insertions(+), 2 deletions(-) create mode 100644 PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs diff --git a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs new file mode 100644 index 0000000..d6bdd28 --- /dev/null +++ b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs @@ -0,0 +1,93 @@ +using System; +using System.Collections.Generic; +using PerfectNumbers.Core; +using FluentAssertions; +using Xunit; + +namespace PerfectNumbers.Core.Tests; + +public class HeuristicCombinedPrimeTesterTests +{ + [Theory] + [Trait("Category", "Fast")] + [InlineData((byte)1)] + [InlineData((byte)7)] + [InlineData((byte)9)] + public void Enumerator_uses_three_a_one_b_schedule_by_default(byte lastDigit) + { + Span groupBState = + stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4]; + + var enumerator = HeuristicCombinedPrimeTester.CreateHeuristicDivisorEnumerator(500, lastDigit, groupBState); + var results = new List(32); + + while (results.Count < 20 && enumerator.TryGetNext(out var candidate)) + { + results.Add(candidate); + } + + results.Count.Should().BeGreaterThanOrEqualTo(16); + results[0].Value.Should().Be(3UL); + results[0].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant); + results[1].Value.Should().Be(7UL); + results[1].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant); + results[2].Value.Should().Be(11UL); + results[2].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant); + results[3].Value.Should().Be(13UL); + results[3].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant); + + int index = 4; + const int segmentsToCheck = 3; + + for (int segment = 0; segment < segmentsToCheck; segment++) + { + for (int i = 0; i < 3; i++) + { + HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupACandidate = results[index++]; + groupACandidate.Value.Should().BeGreaterThan(13UL); + (groupACandidate.Value % 10UL).Should().Be(3UL); + groupACandidate.Group.Should().NotBe(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); + } + + HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupBCandidate = results[index++]; + groupBCandidate.Value.Should().BeGreaterThan(13UL); + (groupBCandidate.Value % 10UL).Should().NotBe(3UL); + groupBCandidate.Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); + } + } + + [Theory] + [Trait("Category", "Fast")] + [InlineData((byte)1)] + [InlineData((byte)7)] + [InlineData((byte)9)] + public void Three_a_one_b_pattern_interleaves_three_group_a_values_before_group_b(byte lastDigit) + { + ReadOnlySpan combined = HeuristicCombinedPrimeTester.GetCombinedDivisors( + lastDigit, + HeuristicCombinedPrimeTester.CombinedDivisorPattern.ThreeAOneB); + + combined.Length.Should().BeGreaterThan(24); + combined[0].Should().Be(3U); + combined[1].Should().Be(7U); + combined[2].Should().Be(11U); + combined[3].Should().Be(13U); + + int index = 4; + const int segmentsToCheck = 4; + + for (int segment = 0; segment < segmentsToCheck; segment++) + { + for (int i = 0; i < 3; i++) + { + uint groupAValue = combined[index++]; + groupAValue.Should().BeGreaterThan(13U); + (groupAValue % 10U).Should().Be(3U); + } + + uint groupBValue = combined[index++]; + groupBValue.Should().BeGreaterThan(13U); + (groupBValue % 10U).Should().NotBe(3U); + } + } +} diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 1942066..4a26976 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -34,12 +34,18 @@ public sealed class HeuristicCombinedPrimeTester private static readonly Lazy CombinedDivisorsEnding7OneAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); private static readonly Lazy CombinedDivisorsEnding9OneAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding1ThreeAOneB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding3ThreeAOneB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding7ThreeAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding9ThreeAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); internal enum CombinedDivisorPattern : byte { TwoAOneB = 0, OneAOneB = 1, + ThreeAOneB = 2, } internal enum HeuristicDivisorGroup : byte @@ -398,7 +404,7 @@ Cleanup: private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { - return GetCombinedDivisors(nMod10, CombinedDivisorPattern.OneAOneB); + return GetCombinedDivisors(nMod10, CombinedDivisorPattern.ThreeAOneB); } internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) @@ -421,6 +427,14 @@ Cleanup: 9 => CombinedDivisorsEnding9OneAOneB.Value, _ => ReadOnlySpan.Empty, }, + CombinedDivisorPattern.ThreeAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeAOneB.Value, + 3 => CombinedDivisorsEnding3ThreeAOneB.Value, + 7 => CombinedDivisorsEnding7ThreeAOneB.Value, + 9 => CombinedDivisorsEnding9ThreeAOneB.Value, + _ => ReadOnlySpan.Empty, + }, _ => ReadOnlySpan.Empty, }; } @@ -440,7 +454,13 @@ Cleanup: } int groupBIndex = GetGroupBStartIndex(groupB); - int groupAInterleaveCount = pattern == CombinedDivisorPattern.OneAOneB ? 1 : 2; + int groupAInterleaveCount = pattern switch + { + CombinedDivisorPattern.OneAOneB => 1, + CombinedDivisorPattern.TwoAOneB => 2, + CombinedDivisorPattern.ThreeAOneB => 3, + _ => 2, + }; while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) { From e8d163f595e09ceaf02483e2ffb404944ee04ead Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 28 Oct 2025 23:22:16 +0100 Subject: [PATCH 05/16] Add three-A-two-B combined divisor schedule --- .../HeuristicCombinedPrimeTesterTests.cs | 51 +++++++++++++++++-- .../HeuristicCombinedPrimeTester.cs | 50 +++++++++++++++--- 2 files changed, 88 insertions(+), 13 deletions(-) diff --git a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs index d6bdd28..a9e94a3 100644 --- a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs +++ b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs @@ -13,7 +13,7 @@ public class HeuristicCombinedPrimeTesterTests [InlineData((byte)1)] [InlineData((byte)7)] [InlineData((byte)9)] - public void Enumerator_uses_three_a_one_b_schedule_by_default(byte lastDigit) + public void Enumerator_uses_three_a_two_b_schedule_by_default(byte lastDigit) { Span groupBState = stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4]; @@ -49,10 +49,51 @@ public class HeuristicCombinedPrimeTesterTests groupACandidate.Group.Should().NotBe(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); } - HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupBCandidate = results[index++]; - groupBCandidate.Value.Should().BeGreaterThan(13UL); - (groupBCandidate.Value % 10UL).Should().NotBe(3UL); - groupBCandidate.Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); + for (int i = 0; i < 2; i++) + { + HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupBCandidate = results[index++]; + groupBCandidate.Value.Should().BeGreaterThan(13UL); + (groupBCandidate.Value % 10UL).Should().NotBe(3UL); + groupBCandidate.Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); + } + } + } + + [Theory] + [Trait("Category", "Fast")] + [InlineData((byte)1)] + [InlineData((byte)7)] + [InlineData((byte)9)] + public void Three_a_two_b_pattern_interleaves_three_group_a_values_before_two_group_b_values(byte lastDigit) + { + ReadOnlySpan combined = HeuristicCombinedPrimeTester.GetCombinedDivisors( + lastDigit, + HeuristicCombinedPrimeTester.CombinedDivisorPattern.ThreeATwoB); + + combined.Length.Should().BeGreaterThan(28); + combined[0].Should().Be(3U); + combined[1].Should().Be(7U); + combined[2].Should().Be(11U); + combined[3].Should().Be(13U); + + int index = 4; + const int segmentsToCheck = 4; + + for (int segment = 0; segment < segmentsToCheck; segment++) + { + for (int i = 0; i < 3; i++) + { + uint groupAValue = combined[index++]; + groupAValue.Should().BeGreaterThan(13U); + (groupAValue % 10U).Should().Be(3U); + } + + for (int i = 0; i < 2; i++) + { + uint groupBValue = combined[index++]; + groupBValue.Should().BeGreaterThan(13U); + (groupBValue % 10U).Should().NotBe(3U); + } } } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 4a26976..2ab9447 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -39,6 +39,11 @@ public sealed class HeuristicCombinedPrimeTester private static readonly Lazy CombinedDivisorsEnding7ThreeAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); private static readonly Lazy CombinedDivisorsEnding9ThreeAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding1ThreeATwoB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding3ThreeATwoB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding7ThreeATwoB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly Lazy CombinedDivisorsEnding9ThreeATwoB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); internal enum CombinedDivisorPattern : byte @@ -46,6 +51,7 @@ public sealed class HeuristicCombinedPrimeTester TwoAOneB = 0, OneAOneB = 1, ThreeAOneB = 2, + ThreeATwoB = 3, } internal enum HeuristicDivisorGroup : byte @@ -404,7 +410,7 @@ Cleanup: private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { - return GetCombinedDivisors(nMod10, CombinedDivisorPattern.ThreeAOneB); + return GetCombinedDivisors(nMod10, CombinedDivisorPattern.ThreeATwoB); } internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) @@ -435,6 +441,14 @@ Cleanup: 9 => CombinedDivisorsEnding9ThreeAOneB.Value, _ => ReadOnlySpan.Empty, }, + CombinedDivisorPattern.ThreeATwoB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeATwoB.Value, + 3 => CombinedDivisorsEnding3ThreeATwoB.Value, + 7 => CombinedDivisorsEnding7ThreeATwoB.Value, + 9 => CombinedDivisorsEnding9ThreeATwoB.Value, + _ => ReadOnlySpan.Empty, + }, _ => ReadOnlySpan.Empty, }; } @@ -454,13 +468,31 @@ Cleanup: } int groupBIndex = GetGroupBStartIndex(groupB); - int groupAInterleaveCount = pattern switch + int groupAInterleaveCount; + int groupBInterleaveCount; + switch (pattern) { - CombinedDivisorPattern.OneAOneB => 1, - CombinedDivisorPattern.TwoAOneB => 2, - CombinedDivisorPattern.ThreeAOneB => 3, - _ => 2, - }; + case CombinedDivisorPattern.OneAOneB: + groupAInterleaveCount = 1; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.TwoAOneB: + groupAInterleaveCount = 2; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.ThreeAOneB: + groupAInterleaveCount = 3; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.ThreeATwoB: + groupAInterleaveCount = 3; + groupBInterleaveCount = 2; + break; + default: + groupAInterleaveCount = 2; + groupBInterleaveCount = 1; + break; + } while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) { @@ -472,10 +504,12 @@ Cleanup: addedA++; } - if (groupBIndex < groupB.Length) + int addedB = 0; + while (addedB < groupBInterleaveCount && groupBIndex < groupB.Length) { combined.Add(groupB[groupBIndex]); groupBIndex++; + addedB++; } } From f993561cd5ec1fe7ed37736f9a27de3e59e0985f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Tue, 28 Oct 2025 23:43:15 +0100 Subject: [PATCH 06/16] Eagerly build combined divisors and restore two-A-one-B default --- .../HeuristicCombinedPrimeTesterTests.cs | 41 ++++++++- .../HeuristicCombinedPrimeTester.cs | 88 +++++++++++-------- 2 files changed, 91 insertions(+), 38 deletions(-) diff --git a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs index a9e94a3..a03ccc5 100644 --- a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs +++ b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs @@ -13,7 +13,7 @@ public class HeuristicCombinedPrimeTesterTests [InlineData((byte)1)] [InlineData((byte)7)] [InlineData((byte)9)] - public void Enumerator_uses_three_a_two_b_schedule_by_default(byte lastDigit) + public void Enumerator_uses_two_a_one_b_schedule_by_default(byte lastDigit) { Span groupBState = stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4]; @@ -41,7 +41,7 @@ public class HeuristicCombinedPrimeTesterTests for (int segment = 0; segment < segmentsToCheck; segment++) { - for (int i = 0; i < 3; i++) + for (int i = 0; i < 2; i++) { HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupACandidate = results[index++]; groupACandidate.Value.Should().BeGreaterThan(13UL); @@ -49,7 +49,7 @@ public class HeuristicCombinedPrimeTesterTests groupACandidate.Group.Should().NotBe(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); } - for (int i = 0; i < 2; i++) + for (int i = 0; i < 1; i++) { HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupBCandidate = results[index++]; groupBCandidate.Value.Should().BeGreaterThan(13UL); @@ -59,6 +59,41 @@ public class HeuristicCombinedPrimeTesterTests } } + [Theory] + [Trait("Category", "Fast")] + [InlineData((byte)1)] + [InlineData((byte)7)] + [InlineData((byte)9)] + public void Two_a_one_b_pattern_interleaves_two_group_a_values_before_group_b(byte lastDigit) + { + ReadOnlySpan combined = HeuristicCombinedPrimeTester.GetCombinedDivisors( + lastDigit, + HeuristicCombinedPrimeTester.CombinedDivisorPattern.TwoAOneB); + + combined.Length.Should().BeGreaterThan(20); + combined[0].Should().Be(3U); + combined[1].Should().Be(7U); + combined[2].Should().Be(11U); + combined[3].Should().Be(13U); + + int index = 4; + const int segmentsToCheck = 5; + + for (int segment = 0; segment < segmentsToCheck; segment++) + { + for (int i = 0; i < 2; i++) + { + uint groupAValue = combined[index++]; + groupAValue.Should().BeGreaterThan(13U); + (groupAValue % 10U).Should().Be(3U); + } + + uint groupBValue = combined[index++]; + groupBValue.Should().BeGreaterThan(13U); + (groupBValue % 10U).Should().NotBe(3U); + } + } + [Theory] [Trait("Category", "Fast")] [InlineData((byte)1)] diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 2ab9447..240180f 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -1,6 +1,4 @@ using System; -using System.Buffers; -using System.Collections.Concurrent; using Open.Numeric.Primes; using System.Collections.Generic; using System.Numerics; @@ -18,31 +16,51 @@ public sealed class HeuristicCombinedPrimeTester static HeuristicCombinedPrimeTester() { HeuristicPrimeSieves.EnsureInitialized(); + + CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB); + + CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB); + + CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB); + + CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB); } private const ulong Wheel210 = 210UL; private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. - private static readonly Lazy CombinedDivisorsEnding1TwoAOneB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding3TwoAOneB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding7TwoAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding9TwoAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly uint[] CombinedDivisorsEnding1TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding3TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding7TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding9TwoAOneB; - private static readonly Lazy CombinedDivisorsEnding1OneAOneB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding3OneAOneB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding7OneAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding9OneAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly uint[] CombinedDivisorsEnding1OneAOneB; + private static readonly uint[] CombinedDivisorsEnding3OneAOneB; + private static readonly uint[] CombinedDivisorsEnding7OneAOneB; + private static readonly uint[] CombinedDivisorsEnding9OneAOneB; - private static readonly Lazy CombinedDivisorsEnding1ThreeAOneB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding3ThreeAOneB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding7ThreeAOneB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding9ThreeAOneB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly uint[] CombinedDivisorsEnding1ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding3ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding7ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding9ThreeAOneB; - private static readonly Lazy CombinedDivisorsEnding1ThreeATwoB = new(() => BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding3ThreeATwoB = new(() => BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding7ThreeATwoB = new(() => BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); - private static readonly Lazy CombinedDivisorsEnding9ThreeATwoB = new(() => BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB), LazyThreadSafetyMode.ExecutionAndPublication); + private static readonly uint[] CombinedDivisorsEnding1ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding3ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding7ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding9ThreeATwoB; private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); @@ -410,7 +428,7 @@ Cleanup: private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { - return GetCombinedDivisors(nMod10, CombinedDivisorPattern.ThreeATwoB); + return GetCombinedDivisors(nMod10, CombinedDivisorPattern.TwoAOneB); } internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) @@ -419,34 +437,34 @@ Cleanup: { CombinedDivisorPattern.TwoAOneB => nMod10 switch { - 1 => CombinedDivisorsEnding1TwoAOneB.Value, - 3 => CombinedDivisorsEnding3TwoAOneB.Value, - 7 => CombinedDivisorsEnding7TwoAOneB.Value, - 9 => CombinedDivisorsEnding9TwoAOneB.Value, + 1 => CombinedDivisorsEnding1TwoAOneB, + 3 => CombinedDivisorsEnding3TwoAOneB, + 7 => CombinedDivisorsEnding7TwoAOneB, + 9 => CombinedDivisorsEnding9TwoAOneB, _ => ReadOnlySpan.Empty, }, CombinedDivisorPattern.OneAOneB => nMod10 switch { - 1 => CombinedDivisorsEnding1OneAOneB.Value, - 3 => CombinedDivisorsEnding3OneAOneB.Value, - 7 => CombinedDivisorsEnding7OneAOneB.Value, - 9 => CombinedDivisorsEnding9OneAOneB.Value, + 1 => CombinedDivisorsEnding1OneAOneB, + 3 => CombinedDivisorsEnding3OneAOneB, + 7 => CombinedDivisorsEnding7OneAOneB, + 9 => CombinedDivisorsEnding9OneAOneB, _ => ReadOnlySpan.Empty, }, CombinedDivisorPattern.ThreeAOneB => nMod10 switch { - 1 => CombinedDivisorsEnding1ThreeAOneB.Value, - 3 => CombinedDivisorsEnding3ThreeAOneB.Value, - 7 => CombinedDivisorsEnding7ThreeAOneB.Value, - 9 => CombinedDivisorsEnding9ThreeAOneB.Value, + 1 => CombinedDivisorsEnding1ThreeAOneB, + 3 => CombinedDivisorsEnding3ThreeAOneB, + 7 => CombinedDivisorsEnding7ThreeAOneB, + 9 => CombinedDivisorsEnding9ThreeAOneB, _ => ReadOnlySpan.Empty, }, CombinedDivisorPattern.ThreeATwoB => nMod10 switch { - 1 => CombinedDivisorsEnding1ThreeATwoB.Value, - 3 => CombinedDivisorsEnding3ThreeATwoB.Value, - 7 => CombinedDivisorsEnding7ThreeATwoB.Value, - 9 => CombinedDivisorsEnding9ThreeATwoB.Value, + 1 => CombinedDivisorsEnding1ThreeATwoB, + 3 => CombinedDivisorsEnding3ThreeATwoB, + 7 => CombinedDivisorsEnding7ThreeATwoB, + 9 => CombinedDivisorsEnding9ThreeATwoB, _ => ReadOnlySpan.Empty, }, _ => ReadOnlySpan.Empty, From b2cb17a96f714e1bbb153c05e8700fa3fe2fe1fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 29 Oct 2025 00:44:39 +0100 Subject: [PATCH 07/16] Return to one-A-one-B combined schedule --- .../HeuristicCombinedPrimeTesterTests.cs | 61 +++++++++++++------ .../HeuristicCombinedPrimeTester.cs | 2 +- 2 files changed, 45 insertions(+), 18 deletions(-) diff --git a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs index a03ccc5..307fdc7 100644 --- a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs +++ b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs @@ -13,7 +13,7 @@ public class HeuristicCombinedPrimeTesterTests [InlineData((byte)1)] [InlineData((byte)7)] [InlineData((byte)9)] - public void Enumerator_uses_two_a_one_b_schedule_by_default(byte lastDigit) + public void Enumerator_uses_one_a_one_b_schedule_by_default(byte lastDigit) { Span groupBState = stackalloc HeuristicCombinedPrimeTester.HeuristicGroupBSequenceState[4]; @@ -26,7 +26,8 @@ public class HeuristicCombinedPrimeTesterTests results.Add(candidate); } - results.Count.Should().BeGreaterThanOrEqualTo(16); + const int segmentsToCheck = 6; + results.Count.Should().BeGreaterThanOrEqualTo(4 + segmentsToCheck * 2); results[0].Value.Should().Be(3UL); results[0].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant); results[1].Value.Should().Be(7UL); @@ -37,25 +38,51 @@ public class HeuristicCombinedPrimeTesterTests results[3].Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupAConstant); int index = 4; - const int segmentsToCheck = 3; for (int segment = 0; segment < segmentsToCheck; segment++) { - for (int i = 0; i < 2; i++) - { - HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupACandidate = results[index++]; - groupACandidate.Value.Should().BeGreaterThan(13UL); - (groupACandidate.Value % 10UL).Should().Be(3UL); - groupACandidate.Group.Should().NotBe(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); - } + HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupACandidate = results[index++]; + groupACandidate.Value.Should().BeGreaterThan(13UL); + (groupACandidate.Value % 10UL).Should().Be(3UL); + groupACandidate.Group.Should().NotBe(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); - for (int i = 0; i < 1; i++) - { - HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupBCandidate = results[index++]; - groupBCandidate.Value.Should().BeGreaterThan(13UL); - (groupBCandidate.Value % 10UL).Should().NotBe(3UL); - groupBCandidate.Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); - } + HeuristicCombinedPrimeTester.HeuristicDivisorCandidate groupBCandidate = results[index++]; + groupBCandidate.Value.Should().BeGreaterThan(13UL); + (groupBCandidate.Value % 10UL).Should().NotBe(3UL); + groupBCandidate.Group.Should().Be(HeuristicCombinedPrimeTester.HeuristicDivisorGroup.GroupB); + } + + } + + [Theory] + [Trait("Category", "Fast")] + [InlineData((byte)1)] + [InlineData((byte)7)] + [InlineData((byte)9)] + public void One_a_one_b_pattern_interleaves_group_a_and_group_b_values(byte lastDigit) + { + ReadOnlySpan combined = HeuristicCombinedPrimeTester.GetCombinedDivisors( + lastDigit, + HeuristicCombinedPrimeTester.CombinedDivisorPattern.OneAOneB); + + combined.Length.Should().BeGreaterThan(18); + combined[0].Should().Be(3U); + combined[1].Should().Be(7U); + combined[2].Should().Be(11U); + combined[3].Should().Be(13U); + + int index = 4; + const int segmentsToCheck = 6; + + for (int segment = 0; segment < segmentsToCheck; segment++) + { + uint groupAValue = combined[index++]; + groupAValue.Should().BeGreaterThan(13U); + (groupAValue % 10U).Should().Be(3U); + + uint groupBValue = combined[index++]; + groupBValue.Should().BeGreaterThan(13U); + (groupBValue % 10U).Should().NotBe(3U); } } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 240180f..88d36c5 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -428,7 +428,7 @@ Cleanup: private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { - return GetCombinedDivisors(nMod10, CombinedDivisorPattern.TwoAOneB); + return GetCombinedDivisors(nMod10, CombinedDivisorPattern.OneAOneB); } internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) From b54f4b6bd609782a29c33cee88bc937b594d3f0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 29 Oct 2025 00:44:46 +0100 Subject: [PATCH 08/16] Simplify combined heuristic divisor selection --- .../HeuristicCombinedPrimeTester.cs | 249 +++--------------- 1 file changed, 39 insertions(+), 210 deletions(-) diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 88d36c5..1ae2880 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -36,6 +36,10 @@ public sealed class HeuristicCombinedPrimeTester CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB); CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB); CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB; + CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB; + CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB; + CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB; } private const ulong Wheel210 = 210UL; @@ -62,6 +66,11 @@ public sealed class HeuristicCombinedPrimeTester private static readonly uint[] CombinedDivisorsEnding7ThreeATwoB; private static readonly uint[] CombinedDivisorsEnding9ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding1; + private static readonly uint[] CombinedDivisorsEnding3; + private static readonly uint[] CombinedDivisorsEnding7; + private static readonly uint[] CombinedDivisorsEnding9; + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); internal enum CombinedDivisorPattern : byte @@ -124,8 +133,7 @@ public sealed class HeuristicCombinedPrimeTester return EvaluateWithOpenNumericFallback(n); } - bool includeGroupB = UseHeuristicGroupBTrialDivision; - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB); + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); } public bool IsPrimeCpu(ulong n) @@ -155,14 +163,14 @@ public sealed class HeuristicCombinedPrimeTester if (!UseHeuristicGroupBTrialDivision) { - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); } bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); return !compositeDetected; } - private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10, bool includeGroupB) + private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10) { ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); if (combinedDivisors.IsEmpty) @@ -170,71 +178,19 @@ public sealed class HeuristicCombinedPrimeTester return EvaluateWithOpenNumericFallback(n); } - int targetGroupACount = CountGroupAEntriesUpTo(sqrtLimit); - int targetGroupBCount = includeGroupB ? CountGroupBEntriesUpTo(sqrtLimit, nMod10) : 0; - - if (targetGroupACount == 0 && (!includeGroupB || targetGroupBCount == 0)) - { - return EvaluateWithOpenNumericFallback(n); - } - - int processedA = 0; - int processedB = 0; - - for (int i = 0; i < combinedDivisors.Length; i++) + int length = combinedDivisors.Length; + for (int i = 0; i < length; i++) { uint entry = combinedDivisors[i]; - HeuristicDivisorGroup group = ResolveGroup(entry); - - if (group == HeuristicDivisorGroup.GroupB) - { - if (!includeGroupB) - { - continue; - } - - if (processedB >= targetGroupBCount) - { - if (processedA >= targetGroupACount) - { - break; - } - - continue; - } - - processedB++; - } - else - { - if (processedA >= targetGroupACount) - { - if (!includeGroupB || processedB >= targetGroupBCount) - { - break; - } - - continue; - } - - processedA++; - } - - ulong divisor = entry; - if (divisor > sqrtLimit) - { - continue; - } - - if (n % divisor == 0UL) - { - return false; - } - - if (processedA >= targetGroupACount && (!includeGroupB || processedB >= targetGroupBCount)) + if (entry > sqrtLimit) { break; } + + if (n % entry == 0UL) + { + return false; + } } return EvaluateWithOpenNumericFallback(n); @@ -295,69 +251,19 @@ public sealed class HeuristicCombinedPrimeTester ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); if (combinedDivisors.IsEmpty) { - compositeDetected = !EvaluateWithOpenNumericFallback(n); goto Cleanup; } - int targetGroupACount = CountGroupAEntriesUpTo(sqrtLimit); - int targetGroupBCount = CountGroupBEntriesUpTo(sqrtLimit, nMod10); - bool includeGroupB = targetGroupBCount > 0; - - if (targetGroupACount == 0 && !includeGroupB) - { - compositeDetected = !EvaluateWithOpenNumericFallback(n); - goto Cleanup; - } - - int processedA = 0; - int processedB = 0; - - for (int i = 0; i < combinedDivisors.Length; i++) + int length = combinedDivisors.Length; + for (int i = 0; i < length; i++) { uint entry = combinedDivisors[i]; - HeuristicDivisorGroup group = ResolveGroup(entry); - - if (group == HeuristicDivisorGroup.GroupB) + if (entry > sqrtLimit) { - if (!includeGroupB) - { - continue; - } - - if (processedB >= targetGroupBCount) - { - if (processedA >= targetGroupACount) - { - break; - } - - continue; - } - - processedB++; - } - else - { - if (processedA >= targetGroupACount) - { - if (!includeGroupB || processedB >= targetGroupBCount) - { - break; - } - - continue; - } - - processedA++; + break; } - ulong divisor = entry; - if (divisor > sqrtLimit) - { - continue; - } - - divisorArray[count] = divisor; + divisorArray[count] = entry; count++; if (count == batchCapacity) @@ -370,11 +276,6 @@ public sealed class HeuristicCombinedPrimeTester count = 0; } - - if (processedA >= targetGroupACount && (!includeGroupB || processedB >= targetGroupBCount)) - { - break; - } } if (!compositeDetected && count > 0) @@ -428,7 +329,14 @@ Cleanup: private static ReadOnlySpan GetCombinedDivisors(byte nMod10) { - return GetCombinedDivisors(nMod10, CombinedDivisorPattern.OneAOneB); + return nMod10 switch + { + 1 => CombinedDivisorsEnding1, + 3 => CombinedDivisorsEnding3, + 7 => CombinedDivisorsEnding7, + 9 => CombinedDivisorsEnding9, + _ => ReadOnlySpan.Empty, + }; } internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) @@ -534,42 +442,7 @@ Cleanup: return combined.ToArray(); } - private static int CountGroupAEntriesUpTo(ulong sqrtLimit) - { - ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; - int count = 0; - while (count < groupA.Length && (ulong)groupA[count] <= sqrtLimit) - { - count++; - } - - return count; - } - - private static int CountGroupBEntriesUpTo(ulong sqrtLimit, byte nMod10) - { - ReadOnlySpan divisors = GetGroupBDivisors(nMod10); - if (divisors.IsEmpty) - { - return 0; - } - - int index = GetGroupBStartIndex(divisors); - int count = 0; - - for (int i = index; i < divisors.Length; i++) - { - if ((ulong)divisors[i] > sqrtLimit) - { - break; - } - - count++; - } - - return count; - } private static int GetGroupBStartIndex(ReadOnlySpan divisors) { @@ -687,24 +560,14 @@ Cleanup: { private readonly ulong sqrtLimit; private readonly ReadOnlySpan combinedDivisors; - private readonly bool includeGroupB; - private readonly int targetGroupACount; - private readonly int targetGroupBCount; private int index; - private int processedA; - private int processedB; public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) { this.sqrtLimit = sqrtLimit; _ = groupBBuffer; combinedDivisors = GetCombinedDivisors(nMod10); - targetGroupACount = CountGroupAEntriesUpTo(sqrtLimit); - targetGroupBCount = CountGroupBEntriesUpTo(sqrtLimit, nMod10); - includeGroupB = targetGroupBCount > 0; index = 0; - processedA = 0; - processedB = 0; } public bool TryGetNext(out HeuristicDivisorCandidate candidate) @@ -712,48 +575,14 @@ Cleanup: while (index < combinedDivisors.Length) { uint entry = combinedDivisors[index++]; + if (entry > sqrtLimit) + { + index = combinedDivisors.Length; + break; + } + HeuristicDivisorGroup group = ResolveGroup(entry); - - if (group == HeuristicDivisorGroup.GroupB) - { - if (!includeGroupB) - { - continue; - } - - if (processedB >= targetGroupBCount) - { - if (processedA >= targetGroupACount) - { - break; - } - - continue; - } - - processedB++; - } - else - { - if (processedA >= targetGroupACount) - { - if (!includeGroupB || processedB >= targetGroupBCount) - { - break; - } - - continue; - } - - processedA++; - } - ulong value = entry; - if (value > sqrtLimit) - { - continue; - } - byte ending = (byte)(value % 10UL); ushort residue = (ushort)(value % Wheel210); candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); From eeeff42be6d319559abefe7a00dc8694002673a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 29 Oct 2025 01:02:31 +0100 Subject: [PATCH 09/16] Replace heuristic sqrt with integer algorithm --- .../HeuristicCombinedPrimeTesterTests.cs | 47 +++++++++++++++++++ .../HeuristicCombinedPrimeTester.cs | 37 ++++++++++----- PerfectNumbers.Core/HeuristicPrimeTester.cs | 37 ++++++++++----- 3 files changed, 95 insertions(+), 26 deletions(-) diff --git a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs index 307fdc7..2d76983 100644 --- a/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs +++ b/PerfectNumbers.Core.Tests/HeuristicCombinedPrimeTesterTests.cs @@ -1,4 +1,5 @@ using System; +using System.Numerics; using System.Collections.Generic; using PerfectNumbers.Core; using FluentAssertions; @@ -193,4 +194,50 @@ public class HeuristicCombinedPrimeTesterTests (groupBValue % 10U).Should().NotBe(3U); } } + + [Theory] + [Trait("Category", "Fast")] + [InlineData(0UL)] + [InlineData(1UL)] + [InlineData(2UL)] + [InlineData(3UL)] + [InlineData(4UL)] + [InlineData(5UL)] + [InlineData(8UL)] + [InlineData(9UL)] + [InlineData(15UL)] + [InlineData(16UL)] + [InlineData(17UL)] + [InlineData(24UL)] + [InlineData(25UL)] + [InlineData(26UL)] + [InlineData(255UL)] + [InlineData(256UL)] + [InlineData(257UL)] + [InlineData(1_000UL)] + [InlineData(1_048_575UL)] + [InlineData(1_048_576UL)] + [InlineData(1_048_577UL)] + [InlineData(4_000_000UL)] + [InlineData(4_000_001UL)] + [InlineData(4_000_004UL)] + [InlineData(65_535_999UL)] + [InlineData(65_536_000UL)] + [InlineData(65_536_001UL)] + [InlineData(4_294_967_295UL)] + [InlineData(4_294_967_296UL)] + [InlineData(4_294_967_297UL)] + [InlineData(18_446_744_073_709_551_615UL)] + [InlineData(18_446_744_073_709_551_614UL)] + [InlineData(18_446_744_073_709_551_613UL)] + public void ComputeHeuristicSqrt_returns_floor_square_root(ulong value) + { + ulong sqrt = HeuristicCombinedPrimeTester.ComputeHeuristicSqrt(value); + + bool squareBelowOrEqual = (UInt128)sqrt * sqrt <= (UInt128)value; + bool nextSquareAbove = (UInt128)(sqrt + 1UL) * (UInt128)(sqrt + 1UL) > (UInt128)value; + + squareBelowOrEqual.Should().BeTrue(); + nextSquareAbove.Should().BeTrue(); + } } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 1ae2880..8d29efc 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -306,25 +306,36 @@ Cleanup: internal static ulong ComputeHeuristicSqrt(ulong n) { - ulong sqrt = (ulong)Math.Sqrt(n); - UInt128 square = (UInt128)sqrt * sqrt; - - while (square > n) + if (n < 2UL) { - sqrt--; - square = (UInt128)sqrt * sqrt; + return n; } - ulong next = sqrt + 1UL; - UInt128 nextSquare = (UInt128)next * next; - while (nextSquare <= n) + ulong bit = 1UL << 62; + while (bit > n) { - sqrt = next; - next++; - nextSquare = (UInt128)next * next; + bit >>= 2; } - return sqrt; + ulong result = 0UL; + ulong remainder = n; + while (bit != 0) + { + ulong candidate = result + bit; + if (remainder >= candidate) + { + remainder -= candidate; + result = (result >> 1) + bit; + } + else + { + result >>= 1; + } + + bit >>= 2; + } + + return result; } private static ReadOnlySpan GetCombinedDivisors(byte nMod10) diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index c3983a1..d8fb2cd 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -431,25 +431,36 @@ Cleanup: internal static ulong ComputeHeuristicSqrt(ulong n) { - ulong sqrt = (ulong)Math.Sqrt(n); - UInt128 square = (UInt128)sqrt * sqrt; - - while (square > n) + if (n < 2UL) { - sqrt--; - square = (UInt128)sqrt * sqrt; + return n; } - ulong next = sqrt + 1UL; - UInt128 nextSquare = (UInt128)next * next; - while (nextSquare <= n) + ulong bit = 1UL << 62; + while (bit > n) { - sqrt = next; - next++; - nextSquare = (UInt128)next * next; + bit >>= 2; } - return sqrt; + ulong result = 0UL; + ulong remainder = n; + while (bit != 0) + { + ulong candidate = result + bit; + if (remainder >= candidate) + { + remainder -= candidate; + result = (result >> 1) + bit; + } + else + { + result >>= 1; + } + + bit >>= 2; + } + + return result; } private static ReadOnlySpan GetGroupBEndingOrder(byte nMod10) => nMod10 switch From 5a9e103bd50908bbbd4f5b69235b766d32ad0c77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 29 Oct 2025 01:29:12 +0100 Subject: [PATCH 10/16] Replace heuristic sqrt with log2-based estimate --- .../HeuristicCombinedPrimeTester.cs | 29 +++++++++---------- PerfectNumbers.Core/HeuristicPrimeTester.cs | 29 +++++++++---------- 2 files changed, 26 insertions(+), 32 deletions(-) diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 8d29efc..05c22a8 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -311,31 +311,28 @@ Cleanup: return n; } - ulong bit = 1UL << 62; - while (bit > n) + int log2 = BitOperations.Log2(n); + ulong estimate = 1UL << ((log2 + 1) >> 1); + ulong quotient = n / estimate; + + while (estimate > quotient) { - bit >>= 2; + estimate = (estimate + quotient) >> 1; + quotient = n / estimate; } - ulong result = 0UL; - ulong remainder = n; - while (bit != 0) + while (true) { - ulong candidate = result + bit; - if (remainder >= candidate) + ulong nextEstimate = estimate + 1UL; + if (nextEstimate > n / nextEstimate) { - remainder -= candidate; - result = (result >> 1) + bit; - } - else - { - result >>= 1; + break; } - bit >>= 2; + estimate = nextEstimate; } - return result; + return estimate; } private static ReadOnlySpan GetCombinedDivisors(byte nMod10) diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index d8fb2cd..d752147 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -436,31 +436,28 @@ Cleanup: return n; } - ulong bit = 1UL << 62; - while (bit > n) + int log2 = BitOperations.Log2(n); + ulong estimate = 1UL << ((log2 + 1) >> 1); + ulong quotient = n / estimate; + + while (estimate > quotient) { - bit >>= 2; + estimate = (estimate + quotient) >> 1; + quotient = n / estimate; } - ulong result = 0UL; - ulong remainder = n; - while (bit != 0) + while (true) { - ulong candidate = result + bit; - if (remainder >= candidate) + ulong nextEstimate = estimate + 1UL; + if (nextEstimate > n / nextEstimate) { - remainder -= candidate; - result = (result >> 1) + bit; - } - else - { - result >>= 1; + break; } - bit >>= 2; + estimate = nextEstimate; } - return result; + return estimate; } private static ReadOnlySpan GetGroupBEndingOrder(byte nMod10) => nMod10 switch From f49963c4f69475dcbf8bc2b5427c7cb737ff84d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 29 Oct 2025 01:43:12 +0100 Subject: [PATCH 11/16] Restore Math.Sqrt heuristic bound --- .../HeuristicCombinedPrimeTester.cs | 34 +++++++------------ PerfectNumbers.Core/HeuristicPrimeTester.cs | 34 +++++++------------ 2 files changed, 26 insertions(+), 42 deletions(-) diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 05c22a8..1ae2880 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -306,33 +306,25 @@ Cleanup: internal static ulong ComputeHeuristicSqrt(ulong n) { - if (n < 2UL) + ulong sqrt = (ulong)Math.Sqrt(n); + UInt128 square = (UInt128)sqrt * sqrt; + + while (square > n) { - return n; + sqrt--; + square = (UInt128)sqrt * sqrt; } - int log2 = BitOperations.Log2(n); - ulong estimate = 1UL << ((log2 + 1) >> 1); - ulong quotient = n / estimate; - - while (estimate > quotient) + ulong next = sqrt + 1UL; + UInt128 nextSquare = (UInt128)next * next; + while (nextSquare <= n) { - estimate = (estimate + quotient) >> 1; - quotient = n / estimate; + sqrt = next; + next++; + nextSquare = (UInt128)next * next; } - while (true) - { - ulong nextEstimate = estimate + 1UL; - if (nextEstimate > n / nextEstimate) - { - break; - } - - estimate = nextEstimate; - } - - return estimate; + return sqrt; } private static ReadOnlySpan GetCombinedDivisors(byte nMod10) diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index d752147..c3983a1 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -431,33 +431,25 @@ Cleanup: internal static ulong ComputeHeuristicSqrt(ulong n) { - if (n < 2UL) + ulong sqrt = (ulong)Math.Sqrt(n); + UInt128 square = (UInt128)sqrt * sqrt; + + while (square > n) { - return n; + sqrt--; + square = (UInt128)sqrt * sqrt; } - int log2 = BitOperations.Log2(n); - ulong estimate = 1UL << ((log2 + 1) >> 1); - ulong quotient = n / estimate; - - while (estimate > quotient) + ulong next = sqrt + 1UL; + UInt128 nextSquare = (UInt128)next * next; + while (nextSquare <= n) { - estimate = (estimate + quotient) >> 1; - quotient = n / estimate; + sqrt = next; + next++; + nextSquare = (UInt128)next * next; } - while (true) - { - ulong nextEstimate = estimate + 1UL; - if (nextEstimate > n / nextEstimate) - { - break; - } - - estimate = nextEstimate; - } - - return estimate; + return sqrt; } private static ReadOnlySpan GetGroupBEndingOrder(byte nMod10) => nMod10 switch From 5d14efb7397167795f9efe3424155d3a88c3fb3a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Wed, 29 Oct 2025 04:31:45 +0100 Subject: [PATCH 12/16] Unroll divisor-by-divisor GPU kernel loops --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 184 ++++++++++++++---- 1 file changed, 149 insertions(+), 35 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index 42f64b5..e18874b 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -1,3 +1,4 @@ +using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; @@ -50,50 +51,136 @@ internal static class DivisorByDivisorKernels hits[offset] = firstHitValue; byte anyHit = firstHitValue; - int endIndex = offset + count; - if (hasCycle) + int remaining = count - 1; + if (remaining > 0) { - for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) + int exponentIndex = offset + 1; + int endIndex = offset + count; + + if (hasCycle) { - ulong exponent = exponents[exponentIndex]; - ulong delta = exponent - previousExponent; - previousExponent = exponent; - - if (delta >= cycleLength) + int blockEnd = exponentIndex + (remaining & ~7); + for (; exponentIndex < blockEnd; exponentIndex += 8) { - delta -= cycleLength; - if (delta >= cycleLength) - { - delta %= cycleLength; - } + ulong exponent0 = exponents[exponentIndex]; + ulong delta0 = exponent0 - previousExponent; + previousExponent = exponent0; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta0, cycleLength); + byte hit0 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent0) : (byte)0; + hits[exponentIndex] = hit0; + + ulong exponent1 = exponents[exponentIndex + 1]; + ulong delta1 = exponent1 - previousExponent; + previousExponent = exponent1; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta1, cycleLength); + byte hit1 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent1) : (byte)0; + hits[exponentIndex + 1] = hit1; + + ulong exponent2 = exponents[exponentIndex + 2]; + ulong delta2 = exponent2 - previousExponent; + previousExponent = exponent2; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta2, cycleLength); + byte hit2 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent2) : (byte)0; + hits[exponentIndex + 2] = hit2; + + ulong exponent3 = exponents[exponentIndex + 3]; + ulong delta3 = exponent3 - previousExponent; + previousExponent = exponent3; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta3, cycleLength); + byte hit3 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent3) : (byte)0; + hits[exponentIndex + 3] = hit3; + + ulong exponent4 = exponents[exponentIndex + 4]; + ulong delta4 = exponent4 - previousExponent; + previousExponent = exponent4; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta4, cycleLength); + byte hit4 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent4) : (byte)0; + hits[exponentIndex + 4] = hit4; + + ulong exponent5 = exponents[exponentIndex + 5]; + ulong delta5 = exponent5 - previousExponent; + previousExponent = exponent5; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta5, cycleLength); + byte hit5 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent5) : (byte)0; + hits[exponentIndex + 5] = hit5; + + ulong exponent6 = exponents[exponentIndex + 6]; + ulong delta6 = exponent6 - previousExponent; + previousExponent = exponent6; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta6, cycleLength); + byte hit6 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent6) : (byte)0; + hits[exponentIndex + 6] = hit6; + + ulong exponent7 = exponents[exponentIndex + 7]; + ulong delta7 = exponent7 - previousExponent; + previousExponent = exponent7; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta7, cycleLength); + byte hit7 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent7) : (byte)0; + hits[exponentIndex + 7] = hit7; + + byte blockHits = (byte)(hit0 | hit1 | hit2 | hit3 | hit4 | hit5 | hit6 | hit7); + anyHit |= blockHits; } - ulong nextRemainder = cycleRemainder + delta; - if (nextRemainder >= cycleLength) + for (; exponentIndex < endIndex; exponentIndex++) { - nextRemainder -= cycleLength; + ulong exponent = exponents[exponentIndex]; + ulong delta = exponent - previousExponent; + previousExponent = exponent; + cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta, cycleLength); + byte hit = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent) : (byte)0; + hits[exponentIndex] = hit; + anyHit |= hit; } - - cycleRemainder = nextRemainder; - - byte hit = 0; - if (cycleRemainder == 0UL) - { - hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; - } - - hits[exponentIndex] = hit; - anyHit |= hit; } - } - else - { - for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) + else { - ulong exponent = exponents[exponentIndex]; - byte hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; - hits[exponentIndex] = hit; - anyHit |= hit; + int blockEnd = exponentIndex + (remaining & ~7); + for (; exponentIndex < blockEnd; exponentIndex += 8) + { + ulong exponent0 = exponents[exponentIndex]; + byte hit0 = ComputeUnityHit(ref stepper, exponent0); + hits[exponentIndex] = hit0; + + ulong exponent1 = exponents[exponentIndex + 1]; + byte hit1 = ComputeUnityHit(ref stepper, exponent1); + hits[exponentIndex + 1] = hit1; + + ulong exponent2 = exponents[exponentIndex + 2]; + byte hit2 = ComputeUnityHit(ref stepper, exponent2); + hits[exponentIndex + 2] = hit2; + + ulong exponent3 = exponents[exponentIndex + 3]; + byte hit3 = ComputeUnityHit(ref stepper, exponent3); + hits[exponentIndex + 3] = hit3; + + ulong exponent4 = exponents[exponentIndex + 4]; + byte hit4 = ComputeUnityHit(ref stepper, exponent4); + hits[exponentIndex + 4] = hit4; + + ulong exponent5 = exponents[exponentIndex + 5]; + byte hit5 = ComputeUnityHit(ref stepper, exponent5); + hits[exponentIndex + 5] = hit5; + + ulong exponent6 = exponents[exponentIndex + 6]; + byte hit6 = ComputeUnityHit(ref stepper, exponent6); + hits[exponentIndex + 6] = hit6; + + ulong exponent7 = exponents[exponentIndex + 7]; + byte hit7 = ComputeUnityHit(ref stepper, exponent7); + hits[exponentIndex + 7] = hit7; + + byte blockHits = (byte)(hit0 | hit1 | hit2 | hit3 | hit4 | hit5 | hit6 | hit7); + anyHit |= blockHits; + } + + for (; exponentIndex < endIndex; exponentIndex++) + { + ulong exponent = exponents[exponentIndex]; + byte hit = ComputeUnityHit(ref stepper, exponent); + hits[exponentIndex] = hit; + anyHit |= hit; + } } } @@ -102,4 +189,31 @@ internal static class DivisorByDivisorKernels Atomic.Min(ref firstHit[0], globalIndex); } } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong AdvanceCycleRemainder(ulong currentRemainder, ulong delta, ulong cycleLength) + { + if (delta >= cycleLength) + { + delta -= cycleLength; + if (delta >= cycleLength) + { + delta %= cycleLength; + } + } + + ulong nextRemainder = currentRemainder + delta; + if (nextRemainder >= cycleLength) + { + nextRemainder -= cycleLength; + } + + return nextRemainder; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static byte ComputeUnityHit(ref ExponentRemainderStepperGpu stepper, ulong exponent) + { + return stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; + } } From 880b1cd96287017afb3ecb74c1ed219a5f53e0bc Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Wed, 29 Oct 2025 04:51:28 +0100 Subject: [PATCH 13/16] Revert "Unroll divisor-by-divisor GPU kernel loops" This reverts commit 5d14efb7397167795f9efe3424155d3a88c3fb3a. --- .../Gpu/Kernels/DivisorByDivisorKernels.cs | 184 ++++-------------- 1 file changed, 35 insertions(+), 149 deletions(-) diff --git a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs index e18874b..42f64b5 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/DivisorByDivisorKernels.cs @@ -1,4 +1,3 @@ -using System.Runtime.CompilerServices; using ILGPU; using ILGPU.Algorithms; using ILGPU.Runtime; @@ -51,136 +50,50 @@ internal static class DivisorByDivisorKernels hits[offset] = firstHitValue; byte anyHit = firstHitValue; - int remaining = count - 1; - if (remaining > 0) + int endIndex = offset + count; + if (hasCycle) { - int exponentIndex = offset + 1; - int endIndex = offset + count; - - if (hasCycle) + for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) { - int blockEnd = exponentIndex + (remaining & ~7); - for (; exponentIndex < blockEnd; exponentIndex += 8) + ulong exponent = exponents[exponentIndex]; + ulong delta = exponent - previousExponent; + previousExponent = exponent; + + if (delta >= cycleLength) { - ulong exponent0 = exponents[exponentIndex]; - ulong delta0 = exponent0 - previousExponent; - previousExponent = exponent0; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta0, cycleLength); - byte hit0 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent0) : (byte)0; - hits[exponentIndex] = hit0; - - ulong exponent1 = exponents[exponentIndex + 1]; - ulong delta1 = exponent1 - previousExponent; - previousExponent = exponent1; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta1, cycleLength); - byte hit1 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent1) : (byte)0; - hits[exponentIndex + 1] = hit1; - - ulong exponent2 = exponents[exponentIndex + 2]; - ulong delta2 = exponent2 - previousExponent; - previousExponent = exponent2; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta2, cycleLength); - byte hit2 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent2) : (byte)0; - hits[exponentIndex + 2] = hit2; - - ulong exponent3 = exponents[exponentIndex + 3]; - ulong delta3 = exponent3 - previousExponent; - previousExponent = exponent3; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta3, cycleLength); - byte hit3 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent3) : (byte)0; - hits[exponentIndex + 3] = hit3; - - ulong exponent4 = exponents[exponentIndex + 4]; - ulong delta4 = exponent4 - previousExponent; - previousExponent = exponent4; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta4, cycleLength); - byte hit4 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent4) : (byte)0; - hits[exponentIndex + 4] = hit4; - - ulong exponent5 = exponents[exponentIndex + 5]; - ulong delta5 = exponent5 - previousExponent; - previousExponent = exponent5; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta5, cycleLength); - byte hit5 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent5) : (byte)0; - hits[exponentIndex + 5] = hit5; - - ulong exponent6 = exponents[exponentIndex + 6]; - ulong delta6 = exponent6 - previousExponent; - previousExponent = exponent6; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta6, cycleLength); - byte hit6 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent6) : (byte)0; - hits[exponentIndex + 6] = hit6; - - ulong exponent7 = exponents[exponentIndex + 7]; - ulong delta7 = exponent7 - previousExponent; - previousExponent = exponent7; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta7, cycleLength); - byte hit7 = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent7) : (byte)0; - hits[exponentIndex + 7] = hit7; - - byte blockHits = (byte)(hit0 | hit1 | hit2 | hit3 | hit4 | hit5 | hit6 | hit7); - anyHit |= blockHits; + delta -= cycleLength; + if (delta >= cycleLength) + { + delta %= cycleLength; + } } - for (; exponentIndex < endIndex; exponentIndex++) + ulong nextRemainder = cycleRemainder + delta; + if (nextRemainder >= cycleLength) { - ulong exponent = exponents[exponentIndex]; - ulong delta = exponent - previousExponent; - previousExponent = exponent; - cycleRemainder = AdvanceCycleRemainder(cycleRemainder, delta, cycleLength); - byte hit = cycleRemainder == 0UL ? ComputeUnityHit(ref stepper, exponent) : (byte)0; - hits[exponentIndex] = hit; - anyHit |= hit; + nextRemainder -= cycleLength; } + + cycleRemainder = nextRemainder; + + byte hit = 0; + if (cycleRemainder == 0UL) + { + hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; + } + + hits[exponentIndex] = hit; + anyHit |= hit; } - else + } + else + { + for (int exponentIndex = offset + 1; exponentIndex < endIndex; exponentIndex++) { - int blockEnd = exponentIndex + (remaining & ~7); - for (; exponentIndex < blockEnd; exponentIndex += 8) - { - ulong exponent0 = exponents[exponentIndex]; - byte hit0 = ComputeUnityHit(ref stepper, exponent0); - hits[exponentIndex] = hit0; - - ulong exponent1 = exponents[exponentIndex + 1]; - byte hit1 = ComputeUnityHit(ref stepper, exponent1); - hits[exponentIndex + 1] = hit1; - - ulong exponent2 = exponents[exponentIndex + 2]; - byte hit2 = ComputeUnityHit(ref stepper, exponent2); - hits[exponentIndex + 2] = hit2; - - ulong exponent3 = exponents[exponentIndex + 3]; - byte hit3 = ComputeUnityHit(ref stepper, exponent3); - hits[exponentIndex + 3] = hit3; - - ulong exponent4 = exponents[exponentIndex + 4]; - byte hit4 = ComputeUnityHit(ref stepper, exponent4); - hits[exponentIndex + 4] = hit4; - - ulong exponent5 = exponents[exponentIndex + 5]; - byte hit5 = ComputeUnityHit(ref stepper, exponent5); - hits[exponentIndex + 5] = hit5; - - ulong exponent6 = exponents[exponentIndex + 6]; - byte hit6 = ComputeUnityHit(ref stepper, exponent6); - hits[exponentIndex + 6] = hit6; - - ulong exponent7 = exponents[exponentIndex + 7]; - byte hit7 = ComputeUnityHit(ref stepper, exponent7); - hits[exponentIndex + 7] = hit7; - - byte blockHits = (byte)(hit0 | hit1 | hit2 | hit3 | hit4 | hit5 | hit6 | hit7); - anyHit |= blockHits; - } - - for (; exponentIndex < endIndex; exponentIndex++) - { - ulong exponent = exponents[exponentIndex]; - byte hit = ComputeUnityHit(ref stepper, exponent); - hits[exponentIndex] = hit; - anyHit |= hit; - } + ulong exponent = exponents[exponentIndex]; + byte hit = stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; + hits[exponentIndex] = hit; + anyHit |= hit; } } @@ -189,31 +102,4 @@ internal static class DivisorByDivisorKernels Atomic.Min(ref firstHit[0], globalIndex); } } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong AdvanceCycleRemainder(ulong currentRemainder, ulong delta, ulong cycleLength) - { - if (delta >= cycleLength) - { - delta -= cycleLength; - if (delta >= cycleLength) - { - delta %= cycleLength; - } - } - - ulong nextRemainder = currentRemainder + delta; - if (nextRemainder >= cycleLength) - { - nextRemainder -= cycleLength; - } - - return nextRemainder; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static byte ComputeUnityHit(ref ExponentRemainderStepperGpu stepper, ulong exponent) - { - return stepper.ComputeNextIsUnityGpu(exponent) ? (byte)1 : (byte)0; - } } From 3147728d5bd8249eba454ffd5262793e8e30ffd4 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Wed, 29 Oct 2025 04:53:13 +0100 Subject: [PATCH 14/16] Benchmarks for heuristic combined methods added --- .../PrimeTesterBenchmarks.cs | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index 3a54d2a..f63e315 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -12,6 +12,7 @@ namespace EvenPerfectBitScanner.Benchmarks; [MemoryDiagnoser] public class PrimeTesterBenchmarks { + private static readonly HeuristicCombinedPrimeTester _heuristicCombinedTester = new(); private static readonly HeuristicPrimeTester _heuristicTester = new(); private ulong[] _candidates = []; @@ -71,6 +72,42 @@ public class PrimeTesterBenchmarks return primeCount; } + [Benchmark] + public int HeuristicCombinedGpu() + { + HeuristicCombinedPrimeTester tester = _heuristicCombinedTester; + ulong[] values = _candidates; + int primeCount = 0; + + for (int i = 0; i < values.Length; i++) + { + if (HeuristicCombinedPrimeTester.IsPrimeGpu(values[i])) + { + primeCount++; + } + } + + return primeCount; + } + + [Benchmark] + public int HeuristicCombinedCpu() + { + HeuristicCombinedPrimeTester tester = _heuristicCombinedTester; + ulong[] values = _candidates; + int primeCount = 0; + + for (int i = 0; i < values.Length; i++) + { + if (HeuristicCombinedPrimeTester.IsPrimeCpu(values[i])) + { + primeCount++; + } + } + + return primeCount; + } + private static readonly PrimeTester _nonHeuristicTester = new(); [Benchmark] From 4c79c4087e1acb0d8dc363f8784de14c0d87da52 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Wed, 29 Oct 2025 14:00:38 +0100 Subject: [PATCH 15/16] Various optimizations --- EvenPerfectBitScanner.Tests/ProgramTests.cs | 2 +- EvenPerfectBitScanner/Program.cs | 2 +- .../PrimeTester/PrimeTesterIsPrimeTests.cs | 2 +- .../Gpu/GpuPrimeWorkLimiter.cs | 29 +- ...MersenneNumberDivisorByDivisorGpuTester.cs | 8 +- .../Gpu/MersenneNumberResidueGpuTester.cs | 4 +- .../HeuristicCombinedPrimeTester.cs | 1334 ++++++++--------- PerfectNumbers.Core/HeuristicPrimeTester.cs | 14 +- PerfectNumbers.Core/PrimeTester.cs | 1013 ++++++------- 9 files changed, 1179 insertions(+), 1229 deletions(-) diff --git a/EvenPerfectBitScanner.Tests/ProgramTests.cs b/EvenPerfectBitScanner.Tests/ProgramTests.cs index 8b1a8fa..f08a676 100644 --- a/EvenPerfectBitScanner.Tests/ProgramTests.cs +++ b/EvenPerfectBitScanner.Tests/ProgramTests.cs @@ -55,7 +55,7 @@ public class ProgramTests [InlineData(81UL, false)] public void IsPrime_identifies_primes_correctly(ulong n, bool expected) { - new PrimeTester().IsPrime(n, CancellationToken.None).Should().Be(expected); + PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected); } [Fact] diff --git a/EvenPerfectBitScanner/Program.cs b/EvenPerfectBitScanner/Program.cs index 48ff15b..86b4c20 100644 --- a/EvenPerfectBitScanner/Program.cs +++ b/EvenPerfectBitScanner/Program.cs @@ -622,7 +622,7 @@ internal static class Program if (!_cliArguments.UseByDivisor) { if (!(GpuContextPool.ForceCpu - ? PrimeTesters.Value!.IsPrime(p, CancellationToken.None) + ? PrimeTester.IsPrime(p, CancellationToken.None) : PrimeTesters.Value!.IsPrimeGpu(p, CancellationToken.None))) { _lastCompositeP = true; diff --git a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs index c99d670..e86e1db 100644 --- a/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeTester/PrimeTesterIsPrimeTests.cs @@ -26,7 +26,7 @@ public class PrimeTesterIsPrimeTests { var tester = new PrimeTester(); - tester.IsPrime(n, CancellationToken.None).Should().Be(expected); + PrimeTester.IsPrime(n, CancellationToken.None).Should().Be(expected); } [Theory] diff --git a/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs b/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs index 97c08ba..1d39b7d 100644 --- a/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs +++ b/PerfectNumbers.Core/Gpu/GpuPrimeWorkLimiter.cs @@ -3,7 +3,6 @@ namespace PerfectNumbers.Core.Gpu; public static class GpuPrimeWorkLimiter { private static SemaphoreSlim _semaphore = new(1, int.MaxValue); - private static readonly object _lock = new(); private static int _currentLimit = 1; public static Lease Acquire() @@ -11,29 +10,21 @@ public static class GpuPrimeWorkLimiter // TODO: Inline the pooled limiter guard from the GpuLimiterThroughputBenchmarks so prime-sieve // GPU jobs stop allocating Lease instances on every acquisition. var sem = _semaphore; - sem.Wait(); + // sem.Wait(); return new Lease(sem); } public static void SetLimit(int value) { - if (value < 1) - { - value = 1; - } + if (value == _currentLimit) + { + return; + } - lock (_lock) - { - if (value == _currentLimit) - { - return; - } - - // TODO: Switch to a shared limiter implementation with GpuWorkLimiter so we can coordinate CPU/GPU prime work using - // the same semaphore pool without reallocating per adjustment. - _semaphore = new SemaphoreSlim(value, value); - _currentLimit = value; - } + // TODO: Switch to a shared limiter implementation with GpuWorkLimiter so we can coordinate CPU/GPU prime work using + // the same semaphore pool without reallocating per adjustment. + _semaphore = new SemaphoreSlim(value, value); + _currentLimit = value; } public sealed class Lease @@ -47,7 +38,7 @@ public static class GpuPrimeWorkLimiter public void Dispose() { - _sem.Release(); + // _sem.Release(); } } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index f480afd..2c9ca8b 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -94,7 +94,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv var gpuLease = RentAccelerator(); var accelerator = gpuLease.Accelerator; - Monitor.Enter(gpuLease.ExecutionLock); + // Monitor.Enter(gpuLease.ExecutionLock); var kernel = GetKernel(accelerator); var resources = RentBatchResources(accelerator, batchCapacity); @@ -121,7 +121,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv out processedCount); ReturnBatchResources(resources); - Monitor.Exit(gpuLease.ExecutionLock); + // Monitor.Exit(gpuLease.ExecutionLock); ReturnAccelerator(gpuLease); if (composite) @@ -560,7 +560,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } } - Monitor.Enter(_lease.ExecutionLock); + // Monitor.Enter(_lease.ExecutionLock); EnsureExecutionResourcesLocked(length); @@ -601,7 +601,7 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv ref byte hitRef = ref MemoryMarshal.GetReference(hitSlice); hitView.CopyToCPU(ref hitRef, length); - Monitor.Exit(_lease.ExecutionLock); + // Monitor.Exit(_lease.ExecutionLock); } public void Dispose() diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs index ca42ddd..759d575 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberResidueGpuTester.cs @@ -30,7 +30,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) var gpuLease = RentAccelerator(); var accelerator = gpuLease.Accelerator; - Monitor.Enter(gpuLease.ExecutionLock); + // Monitor.Enter(gpuLease.ExecutionLock); var stream = accelerator.CreateStream(); var kernel = GetPow2ModKernel(accelerator); @@ -108,7 +108,7 @@ public class MersenneNumberResidueGpuTester(bool useGpuOrder) ReturnResources(accelerator, resources); stream.Dispose(); - Monitor.Exit(gpuLease.ExecutionLock); + // Monitor.Exit(gpuLease.ExecutionLock); ReturnAccelerator(gpuLease); limiter.Dispose(); } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 1ae2880..b9ec773 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -13,682 +13,660 @@ namespace PerfectNumbers.Core; public sealed class HeuristicCombinedPrimeTester { - static HeuristicCombinedPrimeTester() - { - HeuristicPrimeSieves.EnsureInitialized(); - - CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB); - CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB); - - CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB); - CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB); - - CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB); - CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB); - - CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB); - CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB; - CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB; - CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB; - CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB; - } - - private const ulong Wheel210 = 210UL; - - private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. - - private static readonly uint[] CombinedDivisorsEnding1TwoAOneB; - private static readonly uint[] CombinedDivisorsEnding3TwoAOneB; - private static readonly uint[] CombinedDivisorsEnding7TwoAOneB; - private static readonly uint[] CombinedDivisorsEnding9TwoAOneB; - - private static readonly uint[] CombinedDivisorsEnding1OneAOneB; - private static readonly uint[] CombinedDivisorsEnding3OneAOneB; - private static readonly uint[] CombinedDivisorsEnding7OneAOneB; - private static readonly uint[] CombinedDivisorsEnding9OneAOneB; - - private static readonly uint[] CombinedDivisorsEnding1ThreeAOneB; - private static readonly uint[] CombinedDivisorsEnding3ThreeAOneB; - private static readonly uint[] CombinedDivisorsEnding7ThreeAOneB; - private static readonly uint[] CombinedDivisorsEnding9ThreeAOneB; - - private static readonly uint[] CombinedDivisorsEnding1ThreeATwoB; - private static readonly uint[] CombinedDivisorsEnding3ThreeATwoB; - private static readonly uint[] CombinedDivisorsEnding7ThreeATwoB; - private static readonly uint[] CombinedDivisorsEnding9ThreeATwoB; - - private static readonly uint[] CombinedDivisorsEnding1; - private static readonly uint[] CombinedDivisorsEnding3; - private static readonly uint[] CombinedDivisorsEnding7; - private static readonly uint[] CombinedDivisorsEnding9; - - private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); - - internal enum CombinedDivisorPattern : byte - { - TwoAOneB = 0, - OneAOneB = 1, - ThreeAOneB = 2, - ThreeATwoB = 3, - } - - internal enum HeuristicDivisorGroup : byte - { - None = 0, - GroupAConstant = 1, - GroupAWheel = 2, - GroupB = 3, - } - - internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue) - { - public readonly ulong Value = value; - public readonly HeuristicDivisorGroup Group = group; - public readonly byte Ending = ending; - public readonly byte PriorityIndex = priorityIndex; - public readonly ushort WheelResidue = wheelResidue; - } - - internal readonly struct HeuristicDivisorPreparation( - in HeuristicDivisorCandidate candidate, - in MontgomeryDivisorData divisorData, - ulong cycleLengthHint, - bool hasCycleLengthHint) - { - public readonly HeuristicDivisorCandidate Candidate = candidate; - public readonly MontgomeryDivisorData DivisorData = divisorData; - public readonly ulong CycleLengthHint = cycleLengthHint; - public readonly bool HasCycleLengthHint = hasCycleLengthHint; - public readonly bool RequiresCycleComputation = !hasCycleLengthHint; - } - - - public static int HeuristicGpuDivisorBatchSize { get; set; } = 4_096; - - public static int HeuristicDivisorInterleaveBatchSize { get; set; } = 64; - - [ThreadStatic] - private static HeuristicCombinedPrimeTester? _tester; - - public static HeuristicCombinedPrimeTester Exclusive => _tester ??= new(); - - public bool IsPrimeCpu(ulong n, CancellationToken ct) - { - ct.ThrowIfCancellationRequested(); - - byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); - - if (sqrtLimit < 3UL) - { - return EvaluateWithOpenNumericFallback(n); - } - - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); - } - - public bool IsPrimeCpu(ulong n) - { - return IsPrimeCpu(n, CancellationToken.None); - } - - public bool IsPrimeGpu(ulong n) - { - byte nMod10 = (byte)n.Mod10(); - ulong sqrtLimit = ComputeHeuristicSqrt(n); - return IsPrimeGpu(n, sqrtLimit, nMod10); - } - - public bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) - { - return HeuristicIsPrimeGpuCore(n, sqrtLimit, nMod10); - } - - private bool HeuristicIsPrimeGpuCore(ulong n, ulong sqrtLimit, byte nMod10) - { - // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? - if (sqrtLimit < 3UL) - { - return EvaluateWithOpenNumericFallback(n); - } - - if (!UseHeuristicGroupBTrialDivision) - { - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); - } - - bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); - return !compositeDetected; - } - - private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10) - { - ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); - if (combinedDivisors.IsEmpty) - { - return EvaluateWithOpenNumericFallback(n); - } - - int length = combinedDivisors.Length; - for (int i = 0; i < length; i++) - { - uint entry = combinedDivisors[i]; - if (entry > sqrtLimit) - { - break; - } - - if (n % entry == 0UL) - { - return false; - } - } - - return EvaluateWithOpenNumericFallback(n); - } - - private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) - { - int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); - var divisorPool = ThreadStaticPools.UlongPool; - var hitPool = ThreadStaticPools.BytePool; - - ulong[] divisorArray = divisorPool.Rent(batchCapacity); - byte[] hitFlags = hitPool.Rent(batchCapacity); - - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(); - var accelerator = gpu.Accelerator; - var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); - var scratch = state.RentScratch(batchCapacity, accelerator); - - bool compositeDetected = false; - int count = 0; - - lock (gpu.ExecutionLock) - { - var input = scratch.Input; - var output = scratch.Output; - - bool ProcessBatch(int length) - { - if (length <= 0) - { - return false; - } - - input.View.CopyFromCPU(ref divisorArray[0], length); - state.HeuristicTrialDivisionKernel(length, input.View, n, output.View); - accelerator.Synchronize(); - output.View.CopyToCPU(ref hitFlags[0], length); - - for (int i = 0; i < length; i++) - { - if (hitFlags[i] == 0) - { - continue; - } - - ulong divisor = divisorArray[i]; - if (divisor > 1UL && n % divisor == 0UL) - { - return true; - } - } - - return false; - } - - ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); - if (combinedDivisors.IsEmpty) - { - goto Cleanup; - } - - int length = combinedDivisors.Length; - for (int i = 0; i < length; i++) - { - uint entry = combinedDivisors[i]; - if (entry > sqrtLimit) - { - break; - } - - divisorArray[count] = entry; - count++; - - if (count == batchCapacity) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; - } - } - - if (!compositeDetected && count > 0) - { - if (ProcessBatch(count)) - { - compositeDetected = true; - goto Cleanup; - } - - count = 0; - } - } - -Cleanup: - state.ReturnScratch(scratch); - divisorPool.Return(divisorArray); - hitPool.Return(hitFlags); - gpu.Dispose(); - limiter.Dispose(); - return compositeDetected; - } - - private static bool EvaluateWithOpenNumericFallback(ulong n) - { - return Prime.Numbers.IsPrime(n); - } - - internal static ulong ComputeHeuristicSqrt(ulong n) - { - ulong sqrt = (ulong)Math.Sqrt(n); - UInt128 square = (UInt128)sqrt * sqrt; - - while (square > n) - { - sqrt--; - square = (UInt128)sqrt * sqrt; - } - - ulong next = sqrt + 1UL; - UInt128 nextSquare = (UInt128)next * next; - while (nextSquare <= n) - { - sqrt = next; - next++; - nextSquare = (UInt128)next * next; - } - - return sqrt; - } - - private static ReadOnlySpan GetCombinedDivisors(byte nMod10) - { - return nMod10 switch - { - 1 => CombinedDivisorsEnding1, - 3 => CombinedDivisorsEnding3, - 7 => CombinedDivisorsEnding7, - 9 => CombinedDivisorsEnding9, - _ => ReadOnlySpan.Empty, - }; - } - - internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) - { - return pattern switch - { - CombinedDivisorPattern.TwoAOneB => nMod10 switch - { - 1 => CombinedDivisorsEnding1TwoAOneB, - 3 => CombinedDivisorsEnding3TwoAOneB, - 7 => CombinedDivisorsEnding7TwoAOneB, - 9 => CombinedDivisorsEnding9TwoAOneB, - _ => ReadOnlySpan.Empty, - }, - CombinedDivisorPattern.OneAOneB => nMod10 switch - { - 1 => CombinedDivisorsEnding1OneAOneB, - 3 => CombinedDivisorsEnding3OneAOneB, - 7 => CombinedDivisorsEnding7OneAOneB, - 9 => CombinedDivisorsEnding9OneAOneB, - _ => ReadOnlySpan.Empty, - }, - CombinedDivisorPattern.ThreeAOneB => nMod10 switch - { - 1 => CombinedDivisorsEnding1ThreeAOneB, - 3 => CombinedDivisorsEnding3ThreeAOneB, - 7 => CombinedDivisorsEnding7ThreeAOneB, - 9 => CombinedDivisorsEnding9ThreeAOneB, - _ => ReadOnlySpan.Empty, - }, - CombinedDivisorPattern.ThreeATwoB => nMod10 switch - { - 1 => CombinedDivisorsEnding1ThreeATwoB, - 3 => CombinedDivisorsEnding3ThreeATwoB, - 7 => CombinedDivisorsEnding7ThreeATwoB, - 9 => CombinedDivisorsEnding9ThreeATwoB, - _ => ReadOnlySpan.Empty, - }, - _ => ReadOnlySpan.Empty, - }; - } - - private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) - { - ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; - ReadOnlySpan groupB = GetGroupBDivisors(nMod10); - - var combined = new List(groupA.Length + groupB.Length); - - int groupAIndex = 0; - while (groupAIndex < groupA.Length && groupA[groupAIndex] <= 19) - { - combined.Add((uint)groupA[groupAIndex]); - groupAIndex++; - } - - int groupBIndex = GetGroupBStartIndex(groupB); - int groupAInterleaveCount; - int groupBInterleaveCount; - switch (pattern) - { - case CombinedDivisorPattern.OneAOneB: - groupAInterleaveCount = 1; - groupBInterleaveCount = 1; - break; - case CombinedDivisorPattern.TwoAOneB: - groupAInterleaveCount = 2; - groupBInterleaveCount = 1; - break; - case CombinedDivisorPattern.ThreeAOneB: - groupAInterleaveCount = 3; - groupBInterleaveCount = 1; - break; - case CombinedDivisorPattern.ThreeATwoB: - groupAInterleaveCount = 3; - groupBInterleaveCount = 2; - break; - default: - groupAInterleaveCount = 2; - groupBInterleaveCount = 1; - break; - } - - while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) - { - int addedA = 0; - while (addedA < groupAInterleaveCount && groupAIndex < groupA.Length) - { - combined.Add((uint)groupA[groupAIndex]); - groupAIndex++; - addedA++; - } - - int addedB = 0; - while (addedB < groupBInterleaveCount && groupBIndex < groupB.Length) - { - combined.Add(groupB[groupBIndex]); - groupBIndex++; - addedB++; - } - } - - return combined.ToArray(); - } - - - - private static int GetGroupBStartIndex(ReadOnlySpan divisors) - { - int index = 0; - while (index < divisors.Length && IsGroupAPrefixValue(divisors[index])) - { - index++; - } - - return index; - } - - private static bool IsGroupAPrefixValue(uint value) => value is 3U or 7U or 11U or 13U; - - private static HeuristicDivisorGroup ResolveGroup(uint value) => value switch - { - 3U or 7U or 11U or 13U => HeuristicDivisorGroup.GroupAConstant, - _ when value % 10U == 3U => HeuristicDivisorGroup.GroupAWheel, - _ => HeuristicDivisorGroup.GroupB, - }; - - private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch - { - 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, - 3 => ReadOnlySpan.Empty, - 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, - 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, - _ => ReadOnlySpan.Empty, - }; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) - { - return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); - } - - internal struct HeuristicGroupBSequenceState - { - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor) - { - return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) - { - ulong divisor = candidate.Value; - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); - bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); - return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static ulong ResolveHeuristicCycleLength( - ulong exponent, - in HeuristicDivisorPreparation preparation, - out bool cycleFromHint, - out bool cycleComputed, - out bool primeOrderFailed) - { - if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL) - { - cycleFromHint = true; - cycleComputed = true; - primeOrderFailed = false; - return preparation.CycleLengthHint; - } - - ulong divisor = preparation.Candidate.Value; - MontgomeryDivisorData divisorData = preparation.DivisorData; - - bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( - divisor, - exponent, - divisorData, - out ulong computedCycle, - out bool primeOrderFailedLocal); - - if (trySuccess && computedCycle != 0UL) - { - cycleFromHint = false; - cycleComputed = true; - primeOrderFailed = primeOrderFailedLocal; - return computedCycle; - } - - primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL; - ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); - cycleFromHint = false; - cycleComputed = resolvedCycle != 0UL; - return resolvedCycle; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength) - { - ulong[] snapshot = HeuristicSmallCycleSnapshot; - if (divisor < (ulong)snapshot.Length) - { - cycleLength = snapshot[(int)divisor]; - if (cycleLength != 0UL) - { - return true; - } - } - - cycleLength = 0UL; - return false; - } - - internal ref struct HeuristicDivisorEnumerator - { - private readonly ulong sqrtLimit; - private readonly ReadOnlySpan combinedDivisors; - private int index; - - public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) - { - this.sqrtLimit = sqrtLimit; - _ = groupBBuffer; - combinedDivisors = GetCombinedDivisors(nMod10); - index = 0; - } - - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - while (index < combinedDivisors.Length) - { - uint entry = combinedDivisors[index++]; - if (entry > sqrtLimit) - { - index = combinedDivisors.Length; - break; - } - - HeuristicDivisorGroup group = ResolveGroup(entry); - ulong value = entry; - byte ending = (byte)(value % 10UL); - ushort residue = (ushort)(value % Wheel210); - candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); - return true; - } - - candidate = default; - return false; - } - } - - internal ref struct MersenneHeuristicDivisorEnumerator - { - private readonly GpuUInt128 step; - private readonly GpuUInt128 limit; - private GpuUInt128 current; - private MersenneDivisorResidueStepper residueStepper; - private bool active; - private ulong processedCount; - private ulong lastDivisor; - - public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor) - { - var stepLocal = new GpuUInt128(exponent); - stepLocal.ShiftLeft(1); - step = stepLocal; - - var limitLocal = new GpuUInt128(maxDivisor); - limit = limitLocal; - - var firstDivisor = stepLocal; - firstDivisor.Add(1UL); - - bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0; - - current = hasCandidates ? firstDivisor : GpuUInt128.Zero; - residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default; - active = hasCandidates; - processedCount = 0UL; - lastDivisor = 0UL; - } - - public bool TryGetNext(out HeuristicDivisorCandidate candidate) - { - while (active) - { - ulong value = current.Low; - processedCount++; - lastDivisor = value; - - bool admissible = residueStepper.IsAdmissible(); - - Advance(); - - if (admissible) - { - candidate = CreateCandidate(value); - return true; - } - } - - candidate = default; - return false; - } - - private void Advance() - { - GpuUInt128 next = current; - next.Add(step); - - if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0) - { - active = false; - return; - } - - current = next; - residueStepper.Advance(); - } - - public ulong ProcessedCount => processedCount; - - public ulong LastDivisor => lastDivisor; - - public bool Exhausted => !active; - - private static HeuristicDivisorCandidate CreateCandidate(ulong value) - { - byte ending = (byte)(value % 10UL); - HeuristicDivisorGroup group; - - if (value == 3UL || value == 7UL || value == 11UL || value == 13UL) - { - group = HeuristicDivisorGroup.GroupAConstant; - } - else if (ending == 3) - { - group = HeuristicDivisorGroup.GroupAWheel; - } - else - { - group = HeuristicDivisorGroup.GroupB; - } - - ushort residue = (ushort)(value % Wheel210); - return new HeuristicDivisorCandidate(value, group, ending, 0, residue); - } - } + static HeuristicCombinedPrimeTester() + { + HeuristicPrimeSieves.EnsureInitialized(); + + CombinedDivisorsEnding1TwoAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding3TwoAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding7TwoAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.TwoAOneB); + CombinedDivisorsEnding9TwoAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.TwoAOneB); + + CombinedDivisorsEnding1OneAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding3OneAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding7OneAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.OneAOneB); + CombinedDivisorsEnding9OneAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.OneAOneB); + + CombinedDivisorsEnding1ThreeAOneB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding3ThreeAOneB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding7ThreeAOneB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeAOneB); + CombinedDivisorsEnding9ThreeAOneB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeAOneB); + + CombinedDivisorsEnding1ThreeATwoB = BuildCombinedDivisors(1, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding3ThreeATwoB = BuildCombinedDivisors(3, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding7ThreeATwoB = BuildCombinedDivisors(7, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding9ThreeATwoB = BuildCombinedDivisors(9, CombinedDivisorPattern.ThreeATwoB); + CombinedDivisorsEnding1 = CombinedDivisorsEnding1OneAOneB; + CombinedDivisorsEnding3 = CombinedDivisorsEnding3OneAOneB; + CombinedDivisorsEnding7 = CombinedDivisorsEnding7OneAOneB; + CombinedDivisorsEnding9 = CombinedDivisorsEnding9OneAOneB; + } + + private const ulong Wheel210 = 210UL; + + private static readonly bool UseHeuristicGroupBTrialDivision = false; // Temporary fallback gate for Group B. + + private static readonly uint[] CombinedDivisorsEnding1TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding3TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding7TwoAOneB; + private static readonly uint[] CombinedDivisorsEnding9TwoAOneB; + + private static readonly uint[] CombinedDivisorsEnding1OneAOneB; + private static readonly uint[] CombinedDivisorsEnding3OneAOneB; + private static readonly uint[] CombinedDivisorsEnding7OneAOneB; + private static readonly uint[] CombinedDivisorsEnding9OneAOneB; + + private static readonly uint[] CombinedDivisorsEnding1ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding3ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding7ThreeAOneB; + private static readonly uint[] CombinedDivisorsEnding9ThreeAOneB; + + private static readonly uint[] CombinedDivisorsEnding1ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding3ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding7ThreeATwoB; + private static readonly uint[] CombinedDivisorsEnding9ThreeATwoB; + + private static readonly uint[] CombinedDivisorsEnding1; + private static readonly uint[] CombinedDivisorsEnding3; + private static readonly uint[] CombinedDivisorsEnding7; + private static readonly uint[] CombinedDivisorsEnding9; + + private static readonly ulong[] HeuristicSmallCycleSnapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); + + internal enum CombinedDivisorPattern : byte + { + TwoAOneB = 0, + OneAOneB = 1, + ThreeAOneB = 2, + ThreeATwoB = 3, + } + + internal enum HeuristicDivisorGroup : byte + { + None = 0, + GroupAConstant = 1, + GroupAWheel = 2, + GroupB = 3, + } + + internal readonly struct HeuristicDivisorCandidate(ulong value, HeuristicDivisorGroup group, byte ending, byte priorityIndex, ushort wheelResidue) + { + public readonly ulong Value = value; + public readonly HeuristicDivisorGroup Group = group; + public readonly byte Ending = ending; + public readonly byte PriorityIndex = priorityIndex; + public readonly ushort WheelResidue = wheelResidue; + } + + internal readonly struct HeuristicDivisorPreparation( + in HeuristicDivisorCandidate candidate, + in MontgomeryDivisorData divisorData, + ulong cycleLengthHint, + bool hasCycleLengthHint) + { + public readonly HeuristicDivisorCandidate Candidate = candidate; + public readonly MontgomeryDivisorData DivisorData = divisorData; + public readonly ulong CycleLengthHint = cycleLengthHint; + public readonly bool HasCycleLengthHint = hasCycleLengthHint; + public readonly bool RequiresCycleComputation = !hasCycleLengthHint; + } + + + internal static int HeuristicGpuDivisorBatchSize = 4_096; + internal static int HeuristicDivisorInterleaveBatchSize = 64; + + [ThreadStatic] + private static HeuristicCombinedPrimeTester? _tester; + + public static HeuristicCombinedPrimeTester Exclusive => _tester ??= new(); + + public static bool IsPrimeCpu(ulong n) + { + byte nMod10 = (byte)n.Mod10(); + ulong sqrtLimit = ComputeHeuristicSqrt(n); + + if (sqrtLimit < 3UL) + { + return EvaluateWithOpenNumericFallback(n); + } + + return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10); + } + + public static bool IsPrimeGpu(ulong n) + { + byte nMod10 = (byte)n.Mod10(); + ulong sqrtLimit = ComputeHeuristicSqrt(n); + return IsPrimeGpu(n, sqrtLimit, nMod10); + } + + public static bool IsPrimeGpu(ulong n, ulong sqrtLimit, byte nMod10) + { + // TODO: Is this condition ever met on the execution path in EvenPerfectBitScanner? + if (sqrtLimit < 3UL) + { + return EvaluateWithOpenNumericFallback(n); + } + + bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); + return !compositeDetected; + } + + private static bool HeuristicTrialDivisionCpu(ulong n, ulong sqrtLimit, byte nMod10) + { + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + if (combinedDivisors.IsEmpty) + { + return EvaluateWithOpenNumericFallback(n); + } + + int length = combinedDivisors.Length; + for (int i = 0; i < length; i++) + { + uint entry = combinedDivisors[i]; + if (entry > sqrtLimit) + { + break; + } + + if (n % entry == 0UL) + { + return false; + } + } + + return EvaluateWithOpenNumericFallback(n); + } + + private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) + { + int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); + var divisorPool = ThreadStaticPools.UlongPool; + var hitPool = ThreadStaticPools.BytePool; + + ulong[] divisorArray = divisorPool.Rent(batchCapacity); + byte[] hitFlags = hitPool.Rent(batchCapacity); + + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTester.PrimeTesterGpuContextPool.Rent(); + var accelerator = gpu.Accelerator; + var state = PrimeTester.GpuKernelState.GetOrCreate(accelerator); + var scratch = state.RentScratch(batchCapacity, accelerator); + + bool compositeDetected = false; + int count = 0; + + var input = scratch.Input; + var output = scratch.Output; + + bool ProcessBatch(int length) + { + if (length <= 0) + { + return false; + } + + input.View.CopyFromCPU(ref divisorArray[0], length); + state.HeuristicTrialDivisionKernel(length, input.View, n, output.View); + accelerator.Synchronize(); + output.View.CopyToCPU(ref hitFlags[0], length); + + for (int i = 0; i < length; i++) + { + if (hitFlags[i] == 0) + { + continue; + } + + return true; + // ulong divisor = divisorArray[i]; + // if (divisor > 1UL && n % divisor == 0UL) + // { + // return true; + // } + } + + return false; + } + + ReadOnlySpan combinedDivisors = GetCombinedDivisors(nMod10); + if (combinedDivisors.IsEmpty) + { + goto Cleanup; + } + + int length = combinedDivisors.Length; + for (int i = 0; i < length; i++) + { + uint entry = combinedDivisors[i]; + if (entry > sqrtLimit) + { + break; + } + + divisorArray[count] = entry; + count++; + + if (count == batchCapacity) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + + count = 0; + } + } + + if (!compositeDetected && count > 0) + { + if (ProcessBatch(count)) + { + compositeDetected = true; + goto Cleanup; + } + } + + Cleanup: + state.ReturnScratch(scratch); + divisorPool.Return(divisorArray); + hitPool.Return(hitFlags); + gpu.Dispose(); + limiter.Dispose(); + return compositeDetected; + } + + private static bool EvaluateWithOpenNumericFallback(ulong n) + { + return Prime.Numbers.IsPrime(n); + } + + internal static ulong ComputeHeuristicSqrt(ulong n) + { + ulong sqrt = (ulong)Math.Sqrt(n); + UInt128 square = (UInt128)sqrt * sqrt; + + while (square > n) + { + sqrt--; + square = (UInt128)sqrt * sqrt; + } + + ulong next = sqrt + 1UL; + UInt128 nextSquare = (UInt128)next * next; + while (nextSquare <= n) + { + sqrt = next; + next++; + nextSquare = (UInt128)next * next; + } + + return sqrt; + } + + private static ReadOnlySpan GetCombinedDivisors(byte nMod10) + { + return nMod10 switch + { + 1 => CombinedDivisorsEnding1, + 3 => CombinedDivisorsEnding3, + 7 => CombinedDivisorsEnding7, + 9 => CombinedDivisorsEnding9, + _ => ReadOnlySpan.Empty, + }; + } + + internal static ReadOnlySpan GetCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) + { + return pattern switch + { + CombinedDivisorPattern.TwoAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1TwoAOneB, + 3 => CombinedDivisorsEnding3TwoAOneB, + 7 => CombinedDivisorsEnding7TwoAOneB, + 9 => CombinedDivisorsEnding9TwoAOneB, + _ => ReadOnlySpan.Empty, + }, + CombinedDivisorPattern.OneAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1OneAOneB, + 3 => CombinedDivisorsEnding3OneAOneB, + 7 => CombinedDivisorsEnding7OneAOneB, + 9 => CombinedDivisorsEnding9OneAOneB, + _ => ReadOnlySpan.Empty, + }, + CombinedDivisorPattern.ThreeAOneB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeAOneB, + 3 => CombinedDivisorsEnding3ThreeAOneB, + 7 => CombinedDivisorsEnding7ThreeAOneB, + 9 => CombinedDivisorsEnding9ThreeAOneB, + _ => ReadOnlySpan.Empty, + }, + CombinedDivisorPattern.ThreeATwoB => nMod10 switch + { + 1 => CombinedDivisorsEnding1ThreeATwoB, + 3 => CombinedDivisorsEnding3ThreeATwoB, + 7 => CombinedDivisorsEnding7ThreeATwoB, + 9 => CombinedDivisorsEnding9ThreeATwoB, + _ => ReadOnlySpan.Empty, + }, + _ => ReadOnlySpan.Empty, + }; + } + + private static uint[] BuildCombinedDivisors(byte nMod10, CombinedDivisorPattern pattern) + { + ReadOnlySpan groupA = HeuristicPrimeSieves.GroupADivisors; + ReadOnlySpan groupB = GetGroupBDivisors(nMod10); + + var combined = new List(groupA.Length + groupB.Length); + + int groupAIndex = 0; + while (groupAIndex < groupA.Length && groupA[groupAIndex] <= 19) + { + combined.Add((uint)groupA[groupAIndex]); + groupAIndex++; + } + + int groupBIndex = GetGroupBStartIndex(groupB); + int groupAInterleaveCount; + int groupBInterleaveCount; + switch (pattern) + { + case CombinedDivisorPattern.OneAOneB: + groupAInterleaveCount = 1; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.TwoAOneB: + groupAInterleaveCount = 2; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.ThreeAOneB: + groupAInterleaveCount = 3; + groupBInterleaveCount = 1; + break; + case CombinedDivisorPattern.ThreeATwoB: + groupAInterleaveCount = 3; + groupBInterleaveCount = 2; + break; + default: + groupAInterleaveCount = 2; + groupBInterleaveCount = 1; + break; + } + + while (groupAIndex < groupA.Length || groupBIndex < groupB.Length) + { + int addedA = 0; + while (addedA < groupAInterleaveCount && groupAIndex < groupA.Length) + { + combined.Add((uint)groupA[groupAIndex]); + groupAIndex++; + addedA++; + } + + int addedB = 0; + while (addedB < groupBInterleaveCount && groupBIndex < groupB.Length) + { + combined.Add(groupB[groupBIndex]); + groupBIndex++; + addedB++; + } + } + + return combined.ToArray(); + } + + + + private static int GetGroupBStartIndex(ReadOnlySpan divisors) + { + int index = 0; + while (index < divisors.Length && IsGroupAPrefixValue(divisors[index])) + { + index++; + } + + return index; + } + + private static bool IsGroupAPrefixValue(uint value) => value is 3U or 7U or 11U or 13U; + + private static HeuristicDivisorGroup ResolveGroup(uint value) => value switch + { + 3U or 7U or 11U or 13U => HeuristicDivisorGroup.GroupAConstant, + _ when value % 10U == 3U => HeuristicDivisorGroup.GroupAWheel, + _ => HeuristicDivisorGroup.GroupB, + }; + + private static ReadOnlySpan GetGroupBDivisors(byte ending) => ending switch + { + 1 => DivisorGenerator.SmallPrimesLastOneWithoutLastThree, + 3 => ReadOnlySpan.Empty, + 7 => DivisorGenerator.SmallPrimesLastSevenWithoutLastThree, + 9 => DivisorGenerator.SmallPrimesLastNineWithoutLastThree, + _ => ReadOnlySpan.Empty, + }; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorEnumerator CreateHeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + return new HeuristicDivisorEnumerator(sqrtLimit, nMod10, groupBBuffer); + } + + internal struct HeuristicGroupBSequenceState + { + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static MersenneHeuristicDivisorEnumerator CreateMersenneDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + return new MersenneHeuristicDivisorEnumerator(exponent, maxDivisor); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static HeuristicDivisorPreparation PrepareHeuristicDivisor(in HeuristicDivisorCandidate candidate) + { + ulong divisor = candidate.Value; + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + bool hasCycleHint = TryGetCycleLengthHint(divisor, out ulong cycleLength); + return new HeuristicDivisorPreparation(candidate, divisorData, cycleLength, hasCycleHint); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static ulong ResolveHeuristicCycleLength( + ulong exponent, + in HeuristicDivisorPreparation preparation, + out bool cycleFromHint, + out bool cycleComputed, + out bool primeOrderFailed) + { + if (preparation.HasCycleLengthHint && preparation.CycleLengthHint != 0UL) + { + cycleFromHint = true; + cycleComputed = true; + primeOrderFailed = false; + return preparation.CycleLengthHint; + } + + ulong divisor = preparation.Candidate.Value; + MontgomeryDivisorData divisorData = preparation.DivisorData; + + bool trySuccess = MersenneDivisorCycles.TryCalculateCycleLengthForExponentCpu( + divisor, + exponent, + divisorData, + out ulong computedCycle, + out bool primeOrderFailedLocal); + + if (trySuccess && computedCycle != 0UL) + { + cycleFromHint = false; + cycleComputed = true; + primeOrderFailed = primeOrderFailedLocal; + return computedCycle; + } + + primeOrderFailed = primeOrderFailedLocal || !trySuccess || computedCycle == 0UL; + ulong resolvedCycle = MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: primeOrderFailed); + cycleFromHint = false; + cycleComputed = resolvedCycle != 0UL; + return resolvedCycle; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool TryGetCycleLengthHint(ulong divisor, out ulong cycleLength) + { + ulong[] snapshot = HeuristicSmallCycleSnapshot; + if (divisor < (ulong)snapshot.Length) + { + cycleLength = snapshot[(int)divisor]; + if (cycleLength != 0UL) + { + return true; + } + } + + cycleLength = 0UL; + return false; + } + + internal ref struct HeuristicDivisorEnumerator + { + private readonly ulong sqrtLimit; + private readonly ReadOnlySpan combinedDivisors; + private int index; + + public HeuristicDivisorEnumerator(ulong sqrtLimit, byte nMod10, Span groupBBuffer) + { + this.sqrtLimit = sqrtLimit; + _ = groupBBuffer; + combinedDivisors = GetCombinedDivisors(nMod10); + index = 0; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (index < combinedDivisors.Length) + { + uint entry = combinedDivisors[index++]; + if (entry > sqrtLimit) + { + index = combinedDivisors.Length; + break; + } + + HeuristicDivisorGroup group = ResolveGroup(entry); + ulong value = entry; + byte ending = (byte)(value % 10UL); + ushort residue = (ushort)(value % Wheel210); + candidate = new HeuristicDivisorCandidate(value, group, ending, 0, residue); + return true; + } + + candidate = default; + return false; + } + } + + internal ref struct MersenneHeuristicDivisorEnumerator + { + private readonly GpuUInt128 step; + private readonly GpuUInt128 limit; + private GpuUInt128 current; + private MersenneDivisorResidueStepper residueStepper; + private bool active; + private ulong processedCount; + private ulong lastDivisor; + + public MersenneHeuristicDivisorEnumerator(ulong exponent, ulong maxDivisor) + { + var stepLocal = new GpuUInt128(exponent); + stepLocal.ShiftLeft(1); + step = stepLocal; + + var limitLocal = new GpuUInt128(maxDivisor); + limit = limitLocal; + + var firstDivisor = stepLocal; + firstDivisor.Add(1UL); + + bool hasCandidates = !stepLocal.IsZero && firstDivisor.CompareTo(limitLocal) <= 0; + + current = hasCandidates ? firstDivisor : GpuUInt128.Zero; + residueStepper = hasCandidates ? new MersenneDivisorResidueStepper(exponent, stepLocal, firstDivisor) : default; + active = hasCandidates; + processedCount = 0UL; + lastDivisor = 0UL; + } + + public bool TryGetNext(out HeuristicDivisorCandidate candidate) + { + while (active) + { + ulong value = current.Low; + processedCount++; + lastDivisor = value; + + bool admissible = residueStepper.IsAdmissible(); + + Advance(); + + if (admissible) + { + candidate = CreateCandidate(value); + return true; + } + } + + candidate = default; + return false; + } + + private void Advance() + { + GpuUInt128 next = current; + next.Add(step); + + if (next.CompareTo(limit) > 0 || next.CompareTo(current) <= 0) + { + active = false; + return; + } + + current = next; + residueStepper.Advance(); + } + + public ulong ProcessedCount => processedCount; + + public ulong LastDivisor => lastDivisor; + + public bool Exhausted => !active; + + private static HeuristicDivisorCandidate CreateCandidate(ulong value) + { + byte ending = (byte)(value % 10UL); + HeuristicDivisorGroup group; + + if (value == 3UL || value == 7UL || value == 11UL || value == 13UL) + { + group = HeuristicDivisorGroup.GroupAConstant; + } + else if (ending == 3) + { + group = HeuristicDivisorGroup.GroupAWheel; + } + else + { + group = HeuristicDivisorGroup.GroupB; + } + + ushort residue = (ushort)(value % Wheel210); + return new HeuristicDivisorCandidate(value, group, ending, 0, residue); + } + } } diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index c3983a1..ce7e227 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -122,10 +122,10 @@ public sealed class HeuristicPrimeTester return EvaluateWithOpenNumericFallback(n); } - if (!UseHeuristicGroupBTrialDivision) - { - return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); - } + // if (!UseHeuristicGroupBTrialDivision) + // { + // return HeuristicTrialDivisionCpu(n, sqrtLimit, nMod10, includeGroupB: false); + // } bool compositeDetected = HeuristicTrialDivisionGpuDetectsDivisor(n, sqrtLimit, nMod10); return !compositeDetected; @@ -260,8 +260,8 @@ public sealed class HeuristicPrimeTester private bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) { - int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); - int interleaveBatchSize = Math.Max(1, HeuristicDivisorInterleaveBatchSize); + int batchCapacity = HeuristicGpuDivisorBatchSize; + int interleaveBatchSize = HeuristicDivisorInterleaveBatchSize; var divisorPool = ThreadStaticPools.UlongPool; var hitPool = ThreadStaticPools.BytePool; @@ -277,7 +277,7 @@ public sealed class HeuristicPrimeTester bool compositeDetected = false; int count = 0; - lock (gpu.ExecutionLock) + // lock (gpu.ExecutionLock) { var input = scratch.Input; var output = scratch.Output; diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 476eea9..7c4e3af 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -1,11 +1,8 @@ -using System; using System.Buffers; using System.Collections.Concurrent; using System.Numerics; using System.Runtime.CompilerServices; -using System.Threading; using ILGPU; -using ILGPU.Algorithms; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; @@ -13,516 +10,500 @@ namespace PerfectNumbers.Core; public sealed class PrimeTester { - public PrimeTester() - { - } - - [ThreadStatic] - private static PrimeTester? _tester; - - public static PrimeTester Exclusive => _tester ??= new(); - - public bool IsPrime(ulong n, CancellationToken ct) - { - if (n <= 1UL) - { - return false; - } - - if (n == 2UL) - { - throw new InvalidOperationException("PrimeTester.IsPrime encountered the sentinel input 2."); - } - - bool isOdd = (n & 1UL) != 0UL; - bool result = isOdd; - - bool requiresTrialDivision = result && n >= 7UL; - - if (requiresTrialDivision) - { - // EvenPerfectBitScanner streams exponents starting at 136,279,841, so the Mod10/GCD guard never fires on the - // production path. Leave the logic commented out as instrumentation for diagnostic builds. - // bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); - // result &= !sharesMaxExponentFactor; - - if (result) - { - uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; - ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; - - ulong nMod10 = n.Mod10(); - switch (nMod10) - { - case 1UL: - smallPrimeDivisors = PrimesGenerator.SmallPrimesLastOne; - smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastOne; - break; - case 3UL: - smallPrimeDivisors = DivisorGenerator.SmallPrimesLastThree; - smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastThree; - break; - case 7UL: - smallPrimeDivisors = PrimesGenerator.SmallPrimesLastSeven; - smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastSeven; - break; - case 9UL: - smallPrimeDivisors = DivisorGenerator.SmallPrimesLastNine; - smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastNine; - break; - } - - int smallPrimeDivisorsLength = smallPrimeDivisors.Length; - for (int i = 0; i < smallPrimeDivisorsLength; i++) - { - if (smallPrimeDivisorsMul[i] > n) - { - break; - } - - if (n % smallPrimeDivisors[i] == 0) - { - result = false; - break; - } - } - } - } - - return result; - } - - public static bool IsPrimeGpu(ulong n) - { - return Exclusive.IsPrimeGpu(n, CancellationToken.None); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public bool IsPrimeGpu(ulong n, CancellationToken ct) - { - bool forceCpu = GpuContextPool.ForceCpu; - bool belowGpuRange = n < 31UL; - - if (forceCpu || belowGpuRange) - { - return IsPrimeCpu(n, ct); - } - - Span one = stackalloc ulong[1]; - Span outFlags = stackalloc byte[1]; - one[0] = n; - outFlags[0] = 0; - - IsPrimeBatchGpu(one, outFlags); - - return outFlags[0] != 0; - } - - public static bool IsPrimeCpu(ulong n, CancellationToken ct) - { - // Preserve the legacy entry point for callers that bypass the thread-local caches. - return Exclusive.IsPrime(n, ct); - } - - - public static int GpuBatchSize { get; set; } = 262_144; - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static void IsPrimeBatchGpu(ReadOnlySpan values, Span results) - { - var limiter = GpuPrimeWorkLimiter.Acquire(); - var gpu = PrimeTesterGpuContextPool.Rent(); - var state = gpu.State; - var accelerator = gpu.Accelerator; - int totalLength = values.Length; - int batchSize = Math.Max(1, GpuBatchSize); - - lock (gpu.ExecutionLock) - { - var scratch = state.RentScratch(batchSize, accelerator); - var input = scratch.Input; - var output = scratch.Output; - ArrayPool pool = ThreadStaticPools.UlongPool; - ulong[] temp = pool.Rent(batchSize); - - int pos = 0; - while (pos < totalLength) - { - int remaining = totalLength - pos; - int count = remaining > batchSize ? batchSize : remaining; - - values.Slice(pos, count).CopyTo(temp); - input.View.CopyFromCPU(ref temp[0], count); - - state.Kernel( - count, - input.View, - state.DevicePrimesDefault.View, - state.DevicePrimesLastOne.View, - state.DevicePrimesLastSeven.View, - state.DevicePrimesLastThree.View, - state.DevicePrimesLastNine.View, - state.DevicePrimesPow2Default.View, - state.DevicePrimesPow2LastOne.View, - state.DevicePrimesPow2LastSeven.View, - state.DevicePrimesPow2LastThree.View, - state.DevicePrimesPow2LastNine.View, - output.View); - accelerator.Synchronize(); - output.View.CopyToCPU(ref results[pos], count); - - pos += count; - } - - pool.Return(temp, clearArray: false); - state.ReturnScratch(scratch); - } - - gpu.Dispose(); - limiter.Dispose(); - } - - internal static class PrimeTesterGpuContextPool - { - internal sealed class PooledContext - { - public Context Context { get; } - - public Accelerator Accelerator { get; } - - public object ExecutionLock { get; } = new(); - - private KernelState? _kernelState; - - public PooledContext() - { - Context = Context.CreateDefault(); - Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context); - _kernelState = GpuKernelState.GetOrCreate(Accelerator); - } - - public KernelState KernelState - { - get - { - var state = _kernelState; - if (state is null || state.IsDisposed) - { - state = GpuKernelState.GetOrCreate(Accelerator); - _kernelState = state; - } - - return state; - } - } - - public void Dispose() - { - ClearGpuCaches(Accelerator); - _kernelState = null; - Accelerator.Dispose(); - Context.Dispose(); - } - } - - private static readonly ConcurrentQueue Pool = new(); - private static readonly object CreationLock = new(); - - internal static PrimeTesterGpuContextLease Rent() - { - if (Pool.TryDequeue(out var ctx)) - { - return new PrimeTesterGpuContextLease(ctx); - } - - lock (CreationLock) - { - return new PrimeTesterGpuContextLease(new PooledContext()); - } - } - - internal static void DisposeAll() - { - while (Pool.TryDequeue(out var ctx)) - { - ctx.Dispose(); - } - } - - private static void Return(PooledContext ctx) - { - ctx.Accelerator.Synchronize(); - Pool.Enqueue(ctx); - } - - internal struct PrimeTesterGpuContextLease - { - private PooledContext? _ctx; - - internal PrimeTesterGpuContextLease(PooledContext ctx) - { - _ctx = ctx; - } - - public Accelerator Accelerator => _ctx!.Accelerator; - - public object ExecutionLock => _ctx!.ExecutionLock; - - public KernelState State => _ctx!.KernelState; - - public void Dispose() - { - Return(_ctx ?? throw new InvalidOperationException("GPU context lease is not initialized.")); - } - } - } - - // Per-accelerator GPU state for prime sieve (kernel + uploaded primes). - internal sealed class KernelState - { - public Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> Kernel { get; } - public Action, ulong, ArrayView> HeuristicTrialDivisionKernel { get; } - public MemoryBuffer1D DevicePrimesDefault { get; } - public MemoryBuffer1D DevicePrimesLastOne { get; } - public MemoryBuffer1D DevicePrimesLastSeven { get; } - public MemoryBuffer1D DevicePrimesLastThree { get; } - public MemoryBuffer1D DevicePrimesLastNine { get; } - public MemoryBuffer1D DevicePrimesPow2Default { get; } - public MemoryBuffer1D DevicePrimesPow2LastOne { get; } - public MemoryBuffer1D DevicePrimesPow2LastSeven { get; } - public MemoryBuffer1D DevicePrimesPow2LastThree { get; } - public MemoryBuffer1D DevicePrimesPow2LastNine { get; } - public bool IsDisposed { get; private set; } - private readonly System.Collections.Concurrent.ConcurrentBag _scratchPool = []; - // TODO: Replace this ConcurrentBag with the lock-free ring buffer variant validated in - // GpuModularArithmeticBenchmarks so renting scratch buffers stops contending on the bag's internal locks when - // thousands of GPU batches execute per second. - - public KernelState(Accelerator accelerator) - { - IsDisposed = false; - // Compile once per accelerator and upload primes once. - Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); - HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel, ulong, ArrayView>(PrimeTesterKernels.HeuristicTrialDivisionKernel); - - var primesDefault = DivisorGenerator.SmallPrimes; - DevicePrimesDefault = accelerator.Allocate1D(primesDefault.Length); - DevicePrimesDefault.View.CopyFromCPU(primesDefault); - - var primesLastOne = DivisorGenerator.SmallPrimesLastOne; - DevicePrimesLastOne = accelerator.Allocate1D(primesLastOne.Length); - DevicePrimesLastOne.View.CopyFromCPU(primesLastOne); - - var primesLastSeven = DivisorGenerator.SmallPrimesLastSeven; - DevicePrimesLastSeven = accelerator.Allocate1D(primesLastSeven.Length); - DevicePrimesLastSeven.View.CopyFromCPU(primesLastSeven); - - var primesLastThree = DivisorGenerator.SmallPrimesLastThree; - DevicePrimesLastThree = accelerator.Allocate1D(primesLastThree.Length); - DevicePrimesLastThree.View.CopyFromCPU(primesLastThree); - - var primesLastNine = DivisorGenerator.SmallPrimesLastNine; - DevicePrimesLastNine = accelerator.Allocate1D(primesLastNine.Length); - DevicePrimesLastNine.View.CopyFromCPU(primesLastNine); - - var primesPow2Default = DivisorGenerator.SmallPrimesPow2; - DevicePrimesPow2Default = accelerator.Allocate1D(primesPow2Default.Length); - DevicePrimesPow2Default.View.CopyFromCPU(primesPow2Default); - - var primesPow2LastOne = DivisorGenerator.SmallPrimesPow2LastOne; - DevicePrimesPow2LastOne = accelerator.Allocate1D(primesPow2LastOne.Length); - DevicePrimesPow2LastOne.View.CopyFromCPU(primesPow2LastOne); - - var primesPow2LastSeven = DivisorGenerator.SmallPrimesPow2LastSeven; - DevicePrimesPow2LastSeven = accelerator.Allocate1D(primesPow2LastSeven.Length); - DevicePrimesPow2LastSeven.View.CopyFromCPU(primesPow2LastSeven); - - var primesPow2LastThree = DivisorGenerator.SmallPrimesPow2LastThree; - DevicePrimesPow2LastThree = accelerator.Allocate1D(primesPow2LastThree.Length); - DevicePrimesPow2LastThree.View.CopyFromCPU(primesPow2LastThree); - - var primesPow2LastNine = DivisorGenerator.SmallPrimesPow2LastNine; - DevicePrimesPow2LastNine = accelerator.Allocate1D(primesPow2LastNine.Length); - DevicePrimesPow2LastNine.View.CopyFromCPU(primesPow2LastNine); - } - - internal sealed class ScratchBuffers - { - public MemoryBuffer1D Input { get; private set; } - public MemoryBuffer1D Output { get; private set; } - public int Capacity { get; private set; } - - public ScratchBuffers(Accelerator accel, int capacity) - { - Capacity = Math.Max(1, capacity); - Input = accel.Allocate1D(Capacity); - Output = accel.Allocate1D(Capacity); - } - - public void Dispose() - { - Output.Dispose(); - Input.Dispose(); - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public ScratchBuffers RentScratch(int minCapacity, Accelerator accel) - { - minCapacity = Math.Max(1, minCapacity); - while (_scratchPool.TryTake(out var sb)) - { - if (sb.Capacity >= minCapacity) - { - return sb; - } - - sb.Dispose(); - } - - return new ScratchBuffers(accel, minCapacity); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public void ReturnScratch(ScratchBuffers scratch) - { - _scratchPool.Add(scratch); - } - - public void Clear() - { - if (IsDisposed) - { - return; - } - - while (_scratchPool.TryTake(out var sb)) - { - sb.Dispose(); - } - - DevicePrimesDefault.Dispose(); - DevicePrimesLastOne.Dispose(); - DevicePrimesLastSeven.Dispose(); - DevicePrimesLastThree.Dispose(); - DevicePrimesLastNine.Dispose(); - DevicePrimesPow2Default.Dispose(); - DevicePrimesPow2LastOne.Dispose(); - DevicePrimesPow2LastSeven.Dispose(); - DevicePrimesPow2LastThree.Dispose(); - DevicePrimesPow2LastNine.Dispose(); - IsDisposed = true; - } - } - - internal static class GpuKernelState - { - // Map accelerator to cached state; use Lazy to serialize kernel creation - private static readonly System.Collections.Concurrent.ConcurrentDictionary> States = new(); - // TODO: Prewarm this per-accelerator cache during startup (and reuse a simple array keyed by accelerator index) - // once the kernel pool exposes deterministic ordering; the Lazy wrappers showed measurable overhead in the - // GpuModularArithmeticBenchmarks hot path. - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static KernelState GetOrCreate(Accelerator accelerator) - { - var lazy = States.GetOrAdd( - accelerator, - acc => new Lazy(() => new KernelState(acc), System.Threading.LazyThreadSafetyMode.ExecutionAndPublication)); - return lazy.Value; - } - - public static void Clear(Accelerator accelerator) - { - if (States.TryRemove(accelerator, out var lazy)) - { - if (lazy.IsValueCreated) - { - var state = lazy.Value; - state.Clear(); - } - } - } - } - - // Expose cache clearing for accelerator disposal coordination - public static void ClearGpuCaches(Accelerator accelerator) - { - GpuKernelState.Clear(accelerator); - } - - internal static void DisposeGpuContexts() - { - PrimeTesterGpuContextPool.DisposeAll(); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal static bool SharesFactorWithMaxExponent(ulong n) - { - // TODO: Replace this on-the-fly GCD probe with the cached factor table derived from - // ResidueComputationBenchmarks so divisor-cycle metadata can short-circuit the test - // instead of recomputing binary GCD for every candidate. - ulong m = (ulong)BitOperations.Log2(n); - return BinaryGcd(n, m) != 1UL; - } - - internal static void SharesFactorWithMaxExponentBatch(ReadOnlySpan values, Span results) - { - // TODO: Route this batch helper through the shared GPU kernel pool from - // GpuUInt128BinaryGcdBenchmarks so we reuse cached kernels, pinned host buffers, - // and divisor-cycle staging instead of allocating new device buffers per call. - var gpu = PrimeTesterGpuContextPool.Rent(); - var accelerator = gpu.Accelerator; - var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel); - - int length = values.Length; - var inputBuffer = accelerator.Allocate1D(length); - var resultBuffer = accelerator.Allocate1D(length); - - ArrayPool pool = ThreadStaticPools.UlongPool; - ulong[] temp = pool.Rent(length); - values.CopyTo(temp); - inputBuffer.View.CopyFromCPU(ref temp[0], length); - kernel(length, inputBuffer.View, resultBuffer.View); - accelerator.Synchronize(); - resultBuffer.View.CopyToCPU(ref results[0], length); - pool.Return(temp, clearArray: false); - resultBuffer.Dispose(); - inputBuffer.Dispose(); - gpu.Dispose(); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong BinaryGcd(ulong u, ulong v) - { - // TODO: Swap this handwritten binary GCD for the optimized helper measured in - // GpuUInt128BinaryGcdBenchmarks so CPU callers share the faster subtract-less - // ladder once the common implementation is promoted into PerfectNumbers.Core. - if (u == 0UL) - { - return v; - } - - if (v == 0UL) - { - return u; - } - - int shift = BitOperations.TrailingZeroCount(u | v); - u >>= BitOperations.TrailingZeroCount(u); - - do - { - v >>= BitOperations.TrailingZeroCount(v); - if (u > v) - { - (u, v) = (v, u); - } - - v -= u; - } - while (v != 0UL); - - return u << shift; - } - } + public PrimeTester() + { + } + + [ThreadStatic] + private static PrimeTester? _tester; + + public static PrimeTester Exclusive => _tester ??= new(); + + public static bool IsPrime(ulong n, CancellationToken ct) + { + if (n <= 1UL) + { + return false; + } + + if (n == 2UL) + { + throw new InvalidOperationException("PrimeTester.IsPrime encountered the sentinel input 2."); + } + + bool isOdd = (n & 1UL) != 0UL; + bool result = isOdd; + + bool requiresTrialDivision = result && n >= 7UL; + + if (requiresTrialDivision) + { + // EvenPerfectBitScanner streams exponents starting at 136,279,841, so the Mod10/GCD guard never fires on the + // production path. Leave the logic commented out as instrumentation for diagnostic builds. + // bool sharesMaxExponentFactor = n.Mod10() == 1UL && SharesFactorWithMaxExponent(n); + // result &= !sharesMaxExponentFactor; + + if (result) + { + uint[] smallPrimeDivisors = PrimesGenerator.SmallPrimes; + ulong[] smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2; + + ulong nMod10 = n.Mod10(); + switch (nMod10) + { + case 1UL: + smallPrimeDivisors = PrimesGenerator.SmallPrimesLastOne; + smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastOne; + break; + case 3UL: + smallPrimeDivisors = DivisorGenerator.SmallPrimesLastThree; + smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastThree; + break; + case 7UL: + smallPrimeDivisors = PrimesGenerator.SmallPrimesLastSeven; + smallPrimeDivisorsMul = PrimesGenerator.SmallPrimesPow2LastSeven; + break; + case 9UL: + smallPrimeDivisors = DivisorGenerator.SmallPrimesLastNine; + smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastNine; + break; + } + + int smallPrimeDivisorsLength = smallPrimeDivisors.Length; + for (int i = 0; i < smallPrimeDivisorsLength; i++) + { + if (smallPrimeDivisorsMul[i] > n) + { + break; + } + + if (n % smallPrimeDivisors[i] == 0) + { + result = false; + break; + } + } + } + } + + return result; + } + + public static bool IsPrimeGpu(ulong n) + { + return Exclusive.IsPrimeGpu(n, CancellationToken.None); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public bool IsPrimeGpu(ulong n, CancellationToken ct) + { + bool forceCpu = GpuContextPool.ForceCpu; + bool belowGpuRange = n < 31UL; + + if (forceCpu || belowGpuRange) + { + return IsPrimeCpu(n, ct); + } + + Span one = stackalloc ulong[1]; + Span outFlags = stackalloc byte[1]; + one[0] = n; + outFlags[0] = 0; + + IsPrimeBatchGpu(one, outFlags); + + return outFlags[0] != 0; + } + + public static bool IsPrimeCpu(ulong n, CancellationToken ct) + { + // Preserve the legacy entry point for callers that bypass the thread-local caches. + return IsPrime(n, ct); + } + + + public static int GpuBatchSize { get; set; } = 262_144; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void IsPrimeBatchGpu(ReadOnlySpan values, Span results) + { + var limiter = GpuPrimeWorkLimiter.Acquire(); + var gpu = PrimeTesterGpuContextPool.Rent(); + var state = gpu.State; + var accelerator = gpu.Accelerator; + int totalLength = values.Length; + int batchSize = Math.Max(1, GpuBatchSize); + + var scratch = state.RentScratch(batchSize, accelerator); + var input = scratch.Input; + var output = scratch.Output; + ArrayPool pool = ThreadStaticPools.UlongPool; + ulong[] temp = pool.Rent(batchSize); + + int pos = 0; + while (pos < totalLength) + { + int remaining = totalLength - pos; + int count = remaining > batchSize ? batchSize : remaining; + + values.Slice(pos, count).CopyTo(temp); + input.View.CopyFromCPU(ref temp[0], count); + + state.Kernel( + count, + input.View, + state.DevicePrimesDefault.View, + state.DevicePrimesLastOne.View, + state.DevicePrimesLastSeven.View, + state.DevicePrimesLastThree.View, + state.DevicePrimesLastNine.View, + state.DevicePrimesPow2Default.View, + state.DevicePrimesPow2LastOne.View, + state.DevicePrimesPow2LastSeven.View, + state.DevicePrimesPow2LastThree.View, + state.DevicePrimesPow2LastNine.View, + output.View); + accelerator.Synchronize(); + output.View.CopyToCPU(ref results[pos], count); + + pos += count; + } + + pool.Return(temp, clearArray: false); + state.ReturnScratch(scratch); + + gpu.Dispose(); + limiter.Dispose(); + } + + internal static class PrimeTesterGpuContextPool + { + internal sealed class PooledContext + { + public Context Context { get; } + + public Accelerator Accelerator { get; } + + private KernelState? _kernelState; + + public PooledContext() + { + Context = Context.CreateDefault(); + Accelerator = Context.GetPreferredDevice(false).CreateAccelerator(Context); + _kernelState = GpuKernelState.GetOrCreate(Accelerator); + } + + public KernelState KernelState + { + get + { + var state = _kernelState; + if (state is null || state.IsDisposed) + { + state = GpuKernelState.GetOrCreate(Accelerator); + _kernelState = state; + } + + return state; + } + } + + public void Dispose() + { + ClearGpuCaches(Accelerator); + _kernelState = null; + Accelerator.Dispose(); + Context.Dispose(); + } + } + + private static readonly ConcurrentQueue Pool = new(); + + internal static PrimeTesterGpuContextLease Rent() + { + if (Pool.TryDequeue(out var ctx)) + { + return new PrimeTesterGpuContextLease(ctx); + } + + return new PrimeTesterGpuContextLease(new PooledContext()); + } + + internal static void DisposeAll() + { + while (Pool.TryDequeue(out var ctx)) + { + ctx.Dispose(); + } + } + + private static void Return(PooledContext ctx) + { + ctx.Accelerator.Synchronize(); + Pool.Enqueue(ctx); + } + + internal readonly struct PrimeTesterGpuContextLease(PrimeTesterGpuContextPool.PooledContext ctx) + { + private readonly PooledContext? _ctx = ctx; + + public readonly Accelerator Accelerator => _ctx!.Accelerator; + + + public readonly KernelState State => _ctx!.KernelState; + + public readonly void Dispose() + { + Return(_ctx ?? throw new InvalidOperationException("GPU context lease is not initialized.")); + } + } + } + + // Per-accelerator GPU state for prime sieve (kernel + uploaded primes). + internal sealed class KernelState + { + public Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> Kernel; + public Action, ulong, ArrayView> HeuristicTrialDivisionKernel; + public MemoryBuffer1D DevicePrimesDefault; + public MemoryBuffer1D DevicePrimesLastOne; + public MemoryBuffer1D DevicePrimesLastSeven; + public MemoryBuffer1D DevicePrimesLastThree; + public MemoryBuffer1D DevicePrimesLastNine; + public MemoryBuffer1D DevicePrimesPow2Default; + public MemoryBuffer1D DevicePrimesPow2LastOne; + public MemoryBuffer1D DevicePrimesPow2LastSeven; + public MemoryBuffer1D DevicePrimesPow2LastThree; + public MemoryBuffer1D DevicePrimesPow2LastNine; + public bool IsDisposed; + private readonly ConcurrentBag _scratchPool = []; + // TODO: Replace this ConcurrentBag with the lock-free ring buffer variant validated in + // GpuModularArithmeticBenchmarks so renting scratch buffers stops contending on the bag's internal locks when + // thousands of GPU batches execute per second. + + public KernelState(Accelerator accelerator) + { + IsDisposed = false; + // Compile once per accelerator and upload primes once. + Kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel); + HeuristicTrialDivisionKernel = accelerator.LoadAutoGroupedStreamKernel, ulong, ArrayView>(PrimeTesterKernels.HeuristicTrialDivisionKernel); + + var primesDefault = DivisorGenerator.SmallPrimes; + DevicePrimesDefault = accelerator.Allocate1D(primesDefault.Length); + DevicePrimesDefault.View.CopyFromCPU(primesDefault); + + var primesLastOne = DivisorGenerator.SmallPrimesLastOne; + DevicePrimesLastOne = accelerator.Allocate1D(primesLastOne.Length); + DevicePrimesLastOne.View.CopyFromCPU(primesLastOne); + + var primesLastSeven = DivisorGenerator.SmallPrimesLastSeven; + DevicePrimesLastSeven = accelerator.Allocate1D(primesLastSeven.Length); + DevicePrimesLastSeven.View.CopyFromCPU(primesLastSeven); + + var primesLastThree = DivisorGenerator.SmallPrimesLastThree; + DevicePrimesLastThree = accelerator.Allocate1D(primesLastThree.Length); + DevicePrimesLastThree.View.CopyFromCPU(primesLastThree); + + var primesLastNine = DivisorGenerator.SmallPrimesLastNine; + DevicePrimesLastNine = accelerator.Allocate1D(primesLastNine.Length); + DevicePrimesLastNine.View.CopyFromCPU(primesLastNine); + + var primesPow2Default = DivisorGenerator.SmallPrimesPow2; + DevicePrimesPow2Default = accelerator.Allocate1D(primesPow2Default.Length); + DevicePrimesPow2Default.View.CopyFromCPU(primesPow2Default); + + var primesPow2LastOne = DivisorGenerator.SmallPrimesPow2LastOne; + DevicePrimesPow2LastOne = accelerator.Allocate1D(primesPow2LastOne.Length); + DevicePrimesPow2LastOne.View.CopyFromCPU(primesPow2LastOne); + + var primesPow2LastSeven = DivisorGenerator.SmallPrimesPow2LastSeven; + DevicePrimesPow2LastSeven = accelerator.Allocate1D(primesPow2LastSeven.Length); + DevicePrimesPow2LastSeven.View.CopyFromCPU(primesPow2LastSeven); + + var primesPow2LastThree = DivisorGenerator.SmallPrimesPow2LastThree; + DevicePrimesPow2LastThree = accelerator.Allocate1D(primesPow2LastThree.Length); + DevicePrimesPow2LastThree.View.CopyFromCPU(primesPow2LastThree); + + var primesPow2LastNine = DivisorGenerator.SmallPrimesPow2LastNine; + DevicePrimesPow2LastNine = accelerator.Allocate1D(primesPow2LastNine.Length); + DevicePrimesPow2LastNine.View.CopyFromCPU(primesPow2LastNine); + } + + internal sealed class ScratchBuffers + { + public MemoryBuffer1D Input { get; private set; } + public MemoryBuffer1D Output { get; private set; } + public int Capacity { get; private set; } + + public ScratchBuffers(Accelerator accel, int capacity) + { + Capacity = Math.Max(1, capacity); + Input = accel.Allocate1D(Capacity); + Output = accel.Allocate1D(Capacity); + } + + public void Dispose() + { + Output.Dispose(); + Input.Dispose(); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public ScratchBuffers RentScratch(int minCapacity, Accelerator accel) + { + while (_scratchPool.TryTake(out var sb)) + { + if (sb.Capacity >= minCapacity) + { + return sb; + } + + sb.Dispose(); + } + + return new ScratchBuffers(accel, minCapacity); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void ReturnScratch(ScratchBuffers scratch) + { + _scratchPool.Add(scratch); + } + + public void Clear() + { + if (IsDisposed) + { + return; + } + + while (_scratchPool.TryTake(out var sb)) + { + sb.Dispose(); + } + + DevicePrimesDefault.Dispose(); + DevicePrimesLastOne.Dispose(); + DevicePrimesLastSeven.Dispose(); + DevicePrimesLastThree.Dispose(); + DevicePrimesLastNine.Dispose(); + DevicePrimesPow2Default.Dispose(); + DevicePrimesPow2LastOne.Dispose(); + DevicePrimesPow2LastSeven.Dispose(); + DevicePrimesPow2LastThree.Dispose(); + DevicePrimesPow2LastNine.Dispose(); + IsDisposed = true; + } + } + + internal static class GpuKernelState + { + // Map accelerator to cached state; use Lazy to serialize kernel creation + private static readonly System.Collections.Concurrent.ConcurrentDictionary> States = new(); + // TODO: Prewarm this per-accelerator cache during startup (and reuse a simple array keyed by accelerator index) + // once the kernel pool exposes deterministic ordering; the Lazy wrappers showed measurable overhead in the + // GpuModularArithmeticBenchmarks hot path. + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static KernelState GetOrCreate(Accelerator accelerator) + { + var lazy = States.GetOrAdd( + accelerator, + acc => new Lazy(() => new KernelState(acc), LazyThreadSafetyMode.ExecutionAndPublication)); + return lazy.Value; + } + + public static void Clear(Accelerator accelerator) + { + if (States.TryRemove(accelerator, out var lazy)) + { + if (lazy.IsValueCreated) + { + var state = lazy.Value; + state.Clear(); + } + } + } + } + + // Expose cache clearing for accelerator disposal coordination + public static void ClearGpuCaches(Accelerator accelerator) + { + GpuKernelState.Clear(accelerator); + } + + internal static void DisposeGpuContexts() + { + PrimeTesterGpuContextPool.DisposeAll(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static bool SharesFactorWithMaxExponent(ulong n) + { + // TODO: Replace this on-the-fly GCD probe with the cached factor table derived from + // ResidueComputationBenchmarks so divisor-cycle metadata can short-circuit the test + // instead of recomputing binary GCD for every candidate. + ulong m = (ulong)BitOperations.Log2(n); + return BinaryGcd(n, m) != 1UL; + } + + internal static void SharesFactorWithMaxExponentBatch(ReadOnlySpan values, Span results) + { + // TODO: Route this batch helper through the shared GPU kernel pool from + // GpuUInt128BinaryGcdBenchmarks so we reuse cached kernels, pinned host buffers, + // and divisor-cycle staging instead of allocating new device buffers per call. + var gpu = PrimeTesterGpuContextPool.Rent(); + var accelerator = gpu.Accelerator; + var kernel = accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel); + + int length = values.Length; + var inputBuffer = accelerator.Allocate1D(length); + var resultBuffer = accelerator.Allocate1D(length); + + ArrayPool pool = ThreadStaticPools.UlongPool; + ulong[] temp = pool.Rent(length); + values.CopyTo(temp); + inputBuffer.View.CopyFromCPU(ref temp[0], length); + kernel(length, inputBuffer.View, resultBuffer.View); + accelerator.Synchronize(); + resultBuffer.View.CopyToCPU(ref results[0], length); + pool.Return(temp, clearArray: false); + resultBuffer.Dispose(); + inputBuffer.Dispose(); + gpu.Dispose(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong BinaryGcd(ulong u, ulong v) + { + // TODO: Swap this handwritten binary GCD for the optimized helper measured in + // GpuUInt128BinaryGcdBenchmarks so CPU callers share the faster subtract-less + // ladder once the common implementation is promoted into PerfectNumbers.Core. + if (u == 0UL) + { + return v; + } + + if (v == 0UL) + { + return u; + } + + int shift = BitOperations.TrailingZeroCount(u | v); + u >>= BitOperations.TrailingZeroCount(u); + + do + { + v >>= BitOperations.TrailingZeroCount(v); + if (u > v) + { + (u, v) = (v, u); + } + + v -= u; + } + while (v != 0UL); + + return u << shift; + } +} From b2865b54ceb569bf6cde2e75453b0157ab5633b1 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Wed, 29 Oct 2025 14:04:25 +0100 Subject: [PATCH 16/16] Tiny optimization --- PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index b9ec773..6175a28 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -180,7 +180,7 @@ public sealed class HeuristicCombinedPrimeTester private static bool HeuristicTrialDivisionGpuDetectsDivisor(ulong n, ulong sqrtLimit, byte nMod10) { - int batchCapacity = Math.Max(1, HeuristicGpuDivisorBatchSize); + int batchCapacity = HeuristicGpuDivisorBatchSize; var divisorPool = ThreadStaticPools.UlongPool; var hitPool = ThreadStaticPools.BytePool; @@ -201,11 +201,6 @@ public sealed class HeuristicCombinedPrimeTester bool ProcessBatch(int length) { - if (length <= 0) - { - return false; - } - input.View.CopyFromCPU(ref divisorArray[0], length); state.HeuristicTrialDivisionKernel(length, input.View, n, output.View); accelerator.Synchronize();