From c311b27dccd3b1857fa62570f23e004b1b6e3408 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Thu, 27 Nov 2025 16:45:02 +0100 Subject: [PATCH] Prime testing using divisor tables depending on the last digit separated to its own class. PrimeTester switched back to standard endings. --- .../PrimeTesterBenchmarks.cs | 27 +- .../Gpu/Accelerators/AcceleratorStreamPool.cs | 8 +- .../PrimeOrderCalculatorAccelerator.cs | 21 +- .../Kernels/PrimeOrderSpecialMaxKernels.cs | 57 +--- .../Gpu/Kernels/PrimeTesterKernels.cs | 14 +- .../HeuristicCombinedPrimeTester.cs | 5 +- PerfectNumbers.Core/HeuristicPrimeTester.cs | 12 +- PerfectNumbers.Core/PerfectNumberConstants.cs | 6 +- .../PrimeOrderCalculator.Cpu.cs | 268 +++++++++------- .../PrimeOrderCalculator.Gpu.cs | 7 +- PerfectNumbers.Core/PrimeTester.cs | 291 ++++++++++++++++++ PerfectNumbers.Core/ULongExtensions.cs | 25 +- 12 files changed, 526 insertions(+), 215 deletions(-) diff --git a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs index 6c30aac..b27b469 100644 --- a/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/PrimeTesterBenchmarks.cs @@ -26,6 +26,8 @@ public class PrimeTesterBenchmarks return; } + PerfectNumberConstants.RollingAccelerators = 1; + PrimeOrderCalculatorAccelerator.WarmUp(); // HeuristicGroupABPrimeTesterAccelerator.WarmUp(); var gpu = PrimeOrderCalculatorAccelerator.Rent(1); @@ -53,7 +55,7 @@ public class PrimeTesterBenchmarks return primeCount; } - [Benchmark] + // [Benchmark] public int HeuristicGpu() { ulong[] values = _candidates; @@ -72,7 +74,7 @@ public class PrimeTesterBenchmarks return primeCount; } - [Benchmark] + // [Benchmark] public int HeuristicCombinedGpu() { ulong[] values = _candidates; @@ -90,7 +92,7 @@ public class PrimeTesterBenchmarks return primeCount; } - [Benchmark] + // [Benchmark] public int HeuristicCombinedCpu() { ulong[] values = _candidates; @@ -107,7 +109,7 @@ public class PrimeTesterBenchmarks return primeCount; } - [Benchmark] + // [Benchmark] public int NonHeuristicGpu() { ulong[] values = _candidates; @@ -143,6 +145,23 @@ public class PrimeTesterBenchmarks return primeCount; } + [Benchmark] + public int ByLastDigitCpu() + { + ulong[] values = _candidates; + int primeCount = 0; + + for (int i = 0; i < values.Length; i++) + { + if (PrimeTester.IsPrime(values[i])) + { + primeCount++; + } + } + + return primeCount; + } + // [Benchmark] public int OpenNumericPrimes() { diff --git a/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs b/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs index c6a5317..30e79ce 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/AcceleratorStreamPool.cs @@ -16,10 +16,10 @@ public static class AcceleratorStreamPool public static AcceleratorStream Rent(int acceleratorIndex) { - var streamLock = _locks[acceleratorIndex]; + // var streamLock = _locks[acceleratorIndex]; var queue = _streams[acceleratorIndex]; - streamLock.Wait(); + // streamLock.Wait(); return queue.TryDequeue(out var stream) ? stream @@ -28,10 +28,10 @@ public static class AcceleratorStreamPool public static void Return(int acceleratorIndex, AcceleratorStream stream) { - var streamLock = _locks[acceleratorIndex]; + // var streamLock = _locks[acceleratorIndex]; var queue = _streams[acceleratorIndex]; queue.Enqueue(stream); - streamLock.Release(); + // streamLock.Release(); } public static void WarmUp(int acceleratorIndex) diff --git a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs index 98454a4..1ba6153 100644 --- a/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs +++ b/PerfectNumbers.Core/Gpu/Accelerators/PrimeOrderCalculatorAccelerator.cs @@ -115,7 +115,7 @@ public sealed class PrimeOrderCalculatorAccelerator public MemoryBuffer1D OutputUlong; public readonly MemoryBuffer1D OutputUlong2; public ArrayView1D OutputUlongView; - public readonly ArrayView1D OutputUlongView2; + public readonly ArrayView OutputUlongView2; public MemoryBuffer1D OutputInt; public readonly MemoryBuffer1D OutputInt2; @@ -159,7 +159,8 @@ public sealed class PrimeOrderCalculatorAccelerator public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> ConvertToStandardKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelConvertToStandard)).CreateLauncherDelegate, ulong, ulong, ulong, ulong, ulong, ArrayView1D>>(); - public readonly Action, ulong, ulong, HeuristicCombinedGpuViews> HeuristicCombinedTrialDivisionKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, HeuristicCombinedGpuViews>(PrimeTesterKernels.HeuristicTrialCombinedDivisionKernel)).CreateLauncherDelegate, ulong, ulong, HeuristicCombinedGpuViews>>(); + // public readonly Action, ulong, ulong, HeuristicCombinedGpuViews> HeuristicCombinedTrialDivisionKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, HeuristicCombinedGpuViews>(PrimeTesterKernels.HeuristicTrialCombinedDivisionKernel)).CreateLauncherDelegate, ulong, ulong, HeuristicCombinedGpuViews>>(); + public readonly Action, ulong, ulong, ArrayView, ArrayView> HeuristicCombinedTrialDivisionKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ArrayView, ArrayView>(PrimeTesterKernels.HeuristicTrialCombinedDivisionKernel)).CreateLauncherDelegate, ulong, ulong, ArrayView, ArrayView>>(); public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> KeepMontgomeryKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ulong, ulong, ulong, ulong, ulong, ArrayView1D>(Pow2MontgomeryKernels.Pow2MontgomeryKernelKeepMontgomery)).CreateLauncherDelegate, ulong, ulong, ulong, ulong, ulong, ArrayView1D>>(); @@ -173,17 +174,17 @@ public sealed class PrimeOrderCalculatorAccelerator // public readonly Action, ArrayView> SharesFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView>(PrimeTesterKernels.SharesFactorKernel)).CreateLauncherDelegate, ArrayView>>(); - public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan)).CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); + // public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>(SmallPrimeFactorKernels.SmallPrimeFactorKernelScan)).CreateLauncherDelegate, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D>>(); - public readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> SmallPrimeSieveKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel)).CreateLauncherDelegate, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>(); + // public readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> SmallPrimeSieveKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>(PrimeTesterKernels.SmallPrimeSieveKernel)).CreateLauncherDelegate, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView>>(); - public readonly Action, int, ulong, ArrayView1D, ArrayView1D> SpecialMaxKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, int, ulong, ArrayView1D, ArrayView1D>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel)).CreateLauncherDelegate, int, ulong, ArrayView1D, ArrayView1D>>(); + public readonly Action, int, ulong, ArrayView> SpecialMaxKernelLauncher = KernelUtil.GetKernel(accelerator.LoadAutoGroupedStreamKernel, int, ulong, ArrayView>(PrimeOrderGpuHeuristics.EvaluateSpecialMaxCandidatesKernel)).CreateLauncherDelegate, int, ulong, ArrayView>>(); } public readonly Action, ArrayView1D, int, OrderKernelBuffers> CalculateOrderKernelLauncher; public readonly Kernel CheckFactorsKernel; public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> ConvertToStandardKernelLauncher; - public readonly Action, ulong, ulong, HeuristicCombinedGpuViews> HeuristicCombinedTrialDivisionKernelLauncher; + public readonly Action, ulong, ulong, ArrayView, ArrayView> HeuristicCombinedTrialDivisionKernelLauncher; public readonly Action, ulong, ulong, ulong, ulong, ulong, ArrayView1D> KeepMontgomeryKernelLauncher; public readonly Action, ArrayView1D, int, int, ulong, uint, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> PartialFactorKernelLauncher; public readonly Kernel PollardRhoKernel; @@ -192,7 +193,7 @@ public sealed class PrimeOrderCalculatorAccelerator public readonly Action, ArrayView> SharesFactorKernelLauncher; public readonly Action, ArrayView1D, int, ArrayView1D, ArrayView1D, ArrayView1D, ArrayView1D> SmallPrimeFactorKernelLauncher; public readonly Action, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView, ArrayView> SmallPrimeSieveKernelLauncher; - public readonly Action, int, ulong, ArrayView1D, ArrayView1D> SpecialMaxKernelLauncher; + public readonly Action, int, ulong, ArrayView> SpecialMaxKernelLauncher; public void EnsureCalculateOrderWideCapacity(int newSize) { @@ -442,14 +443,14 @@ public sealed class PrimeOrderCalculatorAccelerator CheckFactorsKernel = kernels.CheckFactorsKernel; ConvertToStandardKernelLauncher = kernels.ConvertToStandardKernelLauncher; HeuristicCombinedTrialDivisionKernelLauncher = kernels.HeuristicCombinedTrialDivisionKernelLauncher; - KeepMontgomeryKernelLauncher = kernels.KeepMontgomeryKernelLauncher; + // KeepMontgomeryKernelLauncher = kernels.KeepMontgomeryKernelLauncher; // PartialFactorKernelLauncher = kernels.PartialFactorKernelLauncher; // PollardRhoKernel = kernels.PollardRhoKernel; Pow2ModKernelLauncher = kernels.Pow2ModKernelLauncher; Pow2ModWideKernelLauncher = kernels.Pow2ModWideKernelLauncher; // SharesFactorKernelLauncher = kernels.SharesFactorKernelLauncher; - SmallPrimeFactorKernelLauncher = kernels.SmallPrimeFactorKernelLauncher; - SmallPrimeSieveKernelLauncher = kernels.SmallPrimeSieveKernelLauncher; + // SmallPrimeFactorKernelLauncher = kernels.SmallPrimeFactorKernelLauncher; + // SmallPrimeSieveKernelLauncher = kernels.SmallPrimeSieveKernelLauncher; SpecialMaxKernelLauncher = kernels.SpecialMaxKernelLauncher; // } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs index c785b4a..0e08f1b 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeOrderSpecialMaxKernels.cs @@ -8,72 +8,49 @@ internal static partial class PrimeOrderGpuHeuristics internal static void EvaluateSpecialMaxCandidatesKernel( Index1D index, ulong phi, - ArrayView1D factors, + ArrayView factors, int factorCount, ulong divisorModulus, - ArrayView1D candidates, - ArrayView1D resultOut) + ArrayView resultOut) { - if (index != 0) - { - return; - } - - int actual = 0; for (int i = 0; i < factorCount; i++) { - ulong factor = factors[i]; - if (factor <= 1UL) - { - continue; - } - - ulong reduced = phi / factor; - if (reduced == 0UL) - { - continue; - } - - candidates[actual] = reduced; - actual++; + factors[i] = phi / factors[i]; } - if (actual == 0) + for (int i = 1; i < factorCount; i++) { - resultOut[0] = 1; - return; - } - - for (int i = 1; i < actual; i++) - { - ulong current = candidates[i]; + ulong current = factors[i]; int insert = i - 1; - while (insert >= 0 && candidates[insert] > current) + while (insert >= 0 && factors[insert] > current) { - candidates[insert + 1] = candidates[insert]; + factors[insert + 1] = factors[insert]; insert--; } - candidates[insert + 1] = current; + factors[insert + 1] = current; } var divisorPartial = new GpuDivisorPartialData(divisorModulus); var stepper = new ExponentRemainderStepperGpu(divisorPartial); - if (stepper.InitializeIsUnityGpu(candidates[0])) + if (stepper.InitializeIsUnityGpu(factors[0])) { - Atomic.Exchange(ref resultOut[0], 0); + resultOut[0] = 0UL; + // Atomic.Exchange(ref resultOut[0], 0); return; } - for (int i = 1; i < actual; i++) + for (int i = 1; i < factorCount; i++) { - if (stepper.ComputeNextIsUnityGpu(candidates[i])) + if (stepper.ComputeNextIsUnityGpu(factors[i])) { - Atomic.Exchange(ref resultOut[0], 0); + resultOut[0] = 0UL; + // Atomic.Exchange(ref resultOut[0], 0); return; } } - Atomic.Exchange(ref resultOut[0], 1); + resultOut[0] = 1UL; + // Atomic.Exchange(ref resultOut[0], 1); } } diff --git a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs index 0ede5dd..bef19ca 100644 --- a/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs +++ b/PerfectNumbers.Core/Gpu/Kernels/PrimeTesterKernels.cs @@ -80,21 +80,17 @@ internal static class PrimeTesterKernels public static void HeuristicTrialCombinedDivisionKernel( Index1D index, - byte nMod10, + // byte nMod10, ArrayView resultFlag, ulong n, ulong maxDivisorSquare, - HeuristicCombinedGpuViews tables) + ArrayView divisors, + ArrayView divisorSquares) { - var (divisors, divisorSquares) = tables.SelectDivisorsAndSquares(nMod10); + // var (divisors, divisorSquares) = tables.SelectDivisorsAndSquares(nMod10); int threadIndex = index; - if (divisorSquares[threadIndex] > maxDivisorSquare) - { - return; - } - - if (n % divisors[threadIndex] == 0UL) + if (divisorSquares[threadIndex] <= maxDivisorSquare && (n % divisors[threadIndex] == 0UL)) { Atomic.Exchange(ref resultFlag[0], 1); } diff --git a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs index 7d0db10..aa04dcc 100644 --- a/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicCombinedPrimeTester.cs @@ -230,6 +230,7 @@ public sealed class HeuristicCombinedPrimeTester bool compositeDetected; Span compositeFlag = stackalloc int[1]; // int compositeFlag = 0; + var (divisors, divisorSquares) = gpu.DivisorTables.SelectDivisorsAndSquares(nMod10); var stream = AcceleratorStreamPool.Rent(acceleratorIndex); flagView1DView.MemSetToZero(stream); @@ -237,11 +238,11 @@ public sealed class HeuristicCombinedPrimeTester kernelLauncher( stream, 1, - nMod10, flagView1DView, n, maxDivisorSquare, - gpu.DivisorTables); + divisors, + divisorSquares); flagView1DView.CopyToCPU(stream, compositeFlag); stream.Synchronize(); diff --git a/PerfectNumbers.Core/HeuristicPrimeTester.cs b/PerfectNumbers.Core/HeuristicPrimeTester.cs index e7093b6..1d63c28 100644 --- a/PerfectNumbers.Core/HeuristicPrimeTester.cs +++ b/PerfectNumbers.Core/HeuristicPrimeTester.cs @@ -243,12 +243,13 @@ public sealed class HeuristicPrimeTester { // GpuPrimeWorkLimiter.Acquire(); int acceleratorIndex = gpu.AcceleratorIndex; - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - var flagView1DView = gpu.OutputIntView; - + var kernelLauncher = gpu.HeuristicCombinedTrialDivisionKernelLauncher; bool compositeDetected; int compositeFlag = 0; + var flagView1DView = gpu.OutputIntView; + var (divisors, divisorSquares) = gpu.DivisorTables.SelectDivisorsAndSquares(nMod10); + var stream = AcceleratorStreamPool.Rent(acceleratorIndex); flagView1DView.CopyFromCPU(stream, ref compositeFlag, 1); // Index1D index, @@ -258,15 +259,14 @@ public sealed class HeuristicPrimeTester // HeuristicGpuDivisorTableKind tableKind, // HeuristicGpuDivisorTables tables) - var kernelLauncher = gpu.HeuristicCombinedTrialDivisionKernelLauncher; kernelLauncher( stream, 1, - nMod10, flagView1DView, n, maxDivisorSquare, - gpu.DivisorTables); + divisors, + divisorSquares); flagView1DView.CopyToCPU(stream, ref compositeFlag, 1); stream.Synchronize(); diff --git a/PerfectNumbers.Core/PerfectNumberConstants.cs b/PerfectNumbers.Core/PerfectNumberConstants.cs index 903ed1c..8c9f271 100644 --- a/PerfectNumbers.Core/PerfectNumberConstants.cs +++ b/PerfectNumbers.Core/PerfectNumberConstants.cs @@ -8,14 +8,14 @@ public static class PerfectNumberConstants public const int DefaultSmallPrimeFactorSlotCount = 64; public const int DefaultSpecialMaxFactorCapacity = 1024; public const ulong ExtraDivisorCycleSearchLimit = 64UL; - public const int GpuFrequency = 256 - 1; + public const int GpuFrequency = 192 - 1; public const int MaxQForDivisorCycles = 4_000_000; public static readonly uint PrimesLimit = 1_000_000; //(ulong)Array.MaxLength;// 1_000_000; public const int MaxOddPowersCount = 128; public const int PooledArrayThreshold = 64; public const int Pow2WindowSize = 8; - public static readonly int RollingAccelerators = 298; //SharedGpuContext.Device.MaxNumThreadsPerGroup; - public const int ThreadsByAccelerator = 1; + public static int RollingAccelerators = 298; //SharedGpuContext.Device.MaxNumThreadsPerGroup; + public const int ThreadsByAccelerator = 16; // TODO: Load these limits from the benchmark-driven configuration so CPU and GPU scans stay aligned with the optimal // divisor-cycle datasets we generate offline. // TODO: Promote these magic numbers into a runtime profile derived from EvenPerfectBitScanner.Benchmarks so we can retune diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 1e59cc4..fb3c11f 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -94,31 +94,26 @@ internal static partial class PrimeOrderCalculator return phi; } - // GpuPrimeWorkLimiter.Acquire(); ulong candidateOrder = InitializeStartingOrderCpu(gpu, prime, phi, divisorData); candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); if (TryConfirmOrderCpu(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) { - // GpuPrimeWorkLimiter.Release(); return candidateOrder; } if (config.Mode == PrimeOrderMode.Strict) { orderFactors?.Dispose(); - // GpuPrimeWorkLimiter.Release(); // return CalculateByFactorizationCpu(gpu, prime, divisorData); return CalculateByFactorizationGpu(gpu, prime, divisorData); } if (TryHeuristicFinishCpu(gpu, prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) { - // GpuPrimeWorkLimiter.Release(); return order; } - // GpuPrimeWorkLimiter.Release(); return candidateOrder; } @@ -137,17 +132,16 @@ internal static partial class PrimeOrderCalculator Span stackBuffer = stackalloc ulong[length]; - if (length <= 7) + if (RunOnCpu()) { - return EvaluateSpecialMaxCandidates(stackBuffer, factorSpan, phi, prime, divisorData); + return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData); } return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, prime, divisorData); } - private static bool EvaluateSpecialMaxCandidates(Span buffer, ReadOnlySpan factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData) + private static bool EvaluateSpecialMaxCandidatesCpu(Span buffer, ReadOnlySpan factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData) { - int actual = 0; int factorCount = factors.Length; for (int i = 0; i < factorCount; i++) { @@ -159,36 +153,27 @@ internal static partial class PrimeOrderCalculator // continue; // } - ulong reduced = phi / factor; - if (reduced == 0UL) - { - continue; - } + // if (reduced == 0UL) + // { + // continue; + // } - buffer[actual] = reduced; - actual++; + buffer[i] = phi / factor; } - if (actual == 0) - { - return true; - } - - Span candidates = buffer[..actual]; - candidates.Sort(); + buffer.Sort(); ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - int candidateCount = candidates.Length; - - if (stepper.InitializeCpuIsUnity(candidates[0])) + + if (stepper.InitializeCpuIsUnity(buffer[0])) { ThreadStaticPools.ReturnExponentStepperCpu(stepper); return false; } - for (int i = 1; i < candidateCount; i++) + for (int i = 1; i < factorCount; i++) { - if (stepper.ComputeNextIsUnity(candidates[i])) + if (stepper.ComputeNextIsUnity(buffer[i])) { ThreadStaticPools.ReturnExponentStepperCpu(stepper); return false; @@ -1095,14 +1080,12 @@ internal static partial class PrimeOrderCalculator // Span exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount); Span primeSlots = stackalloc ulong[FactorSlotCount]; Span exponentSlots = stackalloc int[FactorSlotCount]; - ulong[]? primeSlotArray = null; - int[]? exponentSlotArray = null; + ulong[] primeSlotArray = []; + int[] exponentSlotArray = []; // We don't need to worry about leftovers, because we always use indexes within the calculated counts // primeSlots.Clear(); // exponentSlots.Clear(); - int factorCount = 0; - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; ulong remaining = value; @@ -1114,7 +1097,7 @@ internal static partial class PrimeOrderCalculator List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); Stack? compositeStack = null; - PrimeTester primeTester = PrimeTester.Exclusive; + // PrimeTester primeTester = PrimeTester.Exclusive; // HeuristicPrimeTester? tester = _tester ??= new(); // if (!gpuFactored) @@ -1126,7 +1109,7 @@ internal static partial class PrimeOrderCalculator ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, - out factorCount, + out int factorCount, out remaining); // } @@ -1146,72 +1129,21 @@ internal static partial class PrimeOrderCalculator bool pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; if (!pollardRhoDeadlineReached) { - // GpuPrimeWorkLimiter.Acquire(); - - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - // composite will never be smaller than 1 on the execution path - // if (composite <= 1UL) - // { - // continue; - // } - - // bool isPrime = primeTester.HeuristicIsPrimeGpu(composite); - // bool isPrime = PrimeTester.IsPrimeInternal(composite, CancellationToken.None); - - // Atomic.Add(ref _partialFactorHits, 1UL); - // Console.WriteLine($"Partial factor hits {Volatile.Read(ref _partialFactorHits)}"); - - // bool isPrime = PrimeTester.IsPrimeCpu(composite, CancellationToken.None); - // bool goToGpu = Volatile.Read(ref _goToGpu); - // Volatile.Write(ref _goToGpu, !goToGpu); - // bool isPrime = goToGpu - // ? HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite) - // : PrimeTester.IsPrime(composite); - // bool isPrime = composite > PerfectNumberConstants.MaxQForDivisorCycles - // ? HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite) - // : PrimeTester.IsPrime(composite); - // bool isPrime = PrimeTester.IsPrime(composite); - // bool isPrime = HeuristicCombinedPrimeTester.IsPrimeCpu(composite); - // bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); - bool runOnGpu = RunOnGpu(); - bool isPrime = !runOnGpu ? PrimeTester.IsPrime(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); - // bool isPrime = PrimeTester.IsPrime(composite); - // bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - - if (isPrime) - { - AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1); - continue; - } - - long currentTimestamp = Stopwatch.GetTimestamp(); - if (currentTimestamp > deadlineTimestamp) - { - pollardRhoDeadlineReached = true; - pending.Add(new PendingEntry(composite, knownComposite: true)); - continue; - } - - // if (!TryPollardRhoGpu(gpu, composite, out ulong factor)) - if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) - { - pending.Add(new PendingEntry(composite, knownComposite: true)); - continue; - } - - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } + CollectFactors(gpu, ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, pending, compositeStack, deadlineTimestamp, ref pollardRhoDeadlineReached); + // bool runOnCpu = !RunOnGpu(); + // if (runOnCpu) + // { + // CollectFactorsCpu(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, pending, compositeStack, deadlineTimestamp, ref pollardRhoDeadlineReached); + // } + // else + // { + // CollectFactorsGpu(gpu, ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, pending, compositeStack, deadlineTimestamp, ref pollardRhoDeadlineReached); + // } if (!pollardRhoDeadlineReached) { pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; } - - // GpuPrimeWorkLimiter.Release(); } if (pollardRhoDeadlineReached) @@ -1227,7 +1159,6 @@ internal static partial class PrimeOrderCalculator bool cofactorContainsComposite = false; int pendingCount = pending.Count; int index = 0; - // GpuPrimeWorkLimiter.Acquire(); for (; index < pendingCount; index++) { PendingEntry entry = pending[index]; @@ -1382,7 +1313,112 @@ internal static partial class PrimeOrderCalculator return result; } - private static bool RunOnGpu() + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void CollectFactors(PrimeOrderCalculatorAccelerator gpu, ref Span primeSlots, ref Span exponentSlots, ref ulong[] primeSlotArray, ref int[] exponentSlotArray, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, ref bool pollardRhoDeadlineReached) + { + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + bool isPrime = RunOnCpu() ? PrimeTester.IsPrime(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); + + if (isPrime) + { + AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1); + continue; + } + + long currentTimestamp = Stopwatch.GetTimestamp(); + if (currentTimestamp > deadlineTimestamp) + { + pollardRhoDeadlineReached = true; + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } + + // if (!TryPollardRhoGpu(gpu, composite, out ulong factor)) + if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) + { + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } + + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void CollectFactorsCpu(ref Span primeSlots, ref Span exponentSlots, ref ulong[] primeSlotArray, ref int[] exponentSlotArray, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, ref bool pollardRhoDeadlineReached) + { + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + bool isPrime = PrimeTester.IsPrime(composite); + + if (isPrime) + { + AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1); + continue; + } + + long currentTimestamp = Stopwatch.GetTimestamp(); + if (currentTimestamp > deadlineTimestamp) + { + pollardRhoDeadlineReached = true; + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } + + // if (!TryPollardRhoGpu(gpu, composite, out ulong factor)) + if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) + { + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } + + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void CollectFactorsGpu(PrimeOrderCalculatorAccelerator gpu, ref Span primeSlots, ref Span exponentSlots, ref ulong[] primeSlotArray, ref int[] exponentSlotArray, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, ref bool pollardRhoDeadlineReached) + { + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + bool isPrime = HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); + + if (isPrime) + { + AddFactorToCollector(ref primeSlots, ref exponentSlots, ref primeSlotArray, ref exponentSlotArray, ref factorCount, composite, 1); + continue; + } + + long currentTimestamp = Stopwatch.GetTimestamp(); + if (currentTimestamp > deadlineTimestamp) + { + pollardRhoDeadlineReached = true; + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } + + // if (!TryPollardRhoGpu(gpu, composite, out ulong factor)) + if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) + { + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } + + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + } + + private static bool RunOnCpu() { // return false; // return true; @@ -1412,7 +1448,7 @@ internal static partial class PrimeOrderCalculator cpuCount -= PerfectNumberConstants.GpuFrequency; } - return cpuCount == PerfectNumberConstants.GpuFrequency; + return cpuCount != PerfectNumberConstants.GpuFrequency; } private static void PopulateSmallPrimeFactorsCpu( @@ -1420,8 +1456,8 @@ internal static partial class PrimeOrderCalculator uint limit, ref Span primeTargets, ref Span exponentTargets, - ref ulong[]? primeArray, - ref int[]? exponentArray, + ref ulong[] primeArray, + ref int[] exponentArray, out int factorCount, out ulong remaining) { @@ -1592,12 +1628,6 @@ internal static partial class PrimeOrderCalculator long timestamp = 0L; // reused for deadline checks. while (true) { - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - return false; - } - ulong c = (DeterministicRandomCpu.NextUInt64() % (n - 1UL)) + 1UL; ulong x = (DeterministicRandomCpu.NextUInt64() % (n - 2UL)) + 2UL; ulong y = x; @@ -1612,8 +1642,7 @@ internal static partial class PrimeOrderCalculator } x = AdvancePolynomialCpu(x, c, n); - y = AdvancePolynomialCpu(y, c, n); - y = AdvancePolynomialCpu(y, c, n); + y = AdvancePolynomialTwiceCpu(y, c, n); ulong diff = x > y ? x - y : y - x; d = BinaryGcd(diff, n); } @@ -1633,6 +1662,15 @@ internal static partial class PrimeOrderCalculator return (ulong)(value % modulus); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong AdvancePolynomialTwiceCpu(ulong x, ulong c, ulong modulus) + { + UInt128 value = (UInt128)x * x + c; + var remainder = (ulong)(value % modulus); + value = (UInt128)remainder * remainder + c; + return (ulong)(value % modulus); + } + private static ulong BinaryGcd(ulong a, ulong b) { if (a == 0UL) @@ -1679,8 +1717,8 @@ internal static partial class PrimeOrderCalculator private static void AddFactorToCollector( ref Span primes, ref Span exponents, - ref ulong[]? primeArray, - ref int[]? exponentArray, + ref ulong[] primeArray, + ref int[] exponentArray, ref int count, ulong prime, int exponent) @@ -1715,8 +1753,8 @@ internal static partial class PrimeOrderCalculator private static void EnsureCollectorCapacity( ref Span primes, ref Span exponents, - ref ulong[]? primeArray, - ref int[]? exponentArray, + ref ulong[] primeArray, + ref int[] exponentArray, int count, int requiredCapacity) { @@ -1736,10 +1774,10 @@ internal static partial class PrimeOrderCalculator Span newPrimeSpan = newPrimeArray.AsSpan(); Span newExponentSpan = newExponentArray.AsSpan(); - primes.Slice(0, count).CopyTo(newPrimeSpan); - exponents.Slice(0, count).CopyTo(newExponentSpan); + primes[..count].CopyTo(newPrimeSpan); + exponents[..count].CopyTo(newExponentSpan); - if (primeArray is not null) + if (primeArray.Length > 0) { ThreadStaticPools.UlongPool.Return(primeArray, clearArray: false); ThreadStaticPools.IntPool.Return(exponentArray!, clearArray: false); diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index dd31835..753a6c6 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -41,7 +41,6 @@ internal static partial class PrimeOrderCalculator gpu.EnsureCapacity(entryCount, 1); var pow2ModEntriesToTestOnDeviceView = gpu.Pow2ModEntriesToTestOnDeviceView; - // GpuPrimeWorkLimiter.Acquire(); AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); pow2ModEntriesToTestOnDeviceView.CopyFromCPU(stream, entries); @@ -51,7 +50,6 @@ internal static partial class PrimeOrderCalculator stream.Synchronize(); AcceleratorStreamPool.Return(acceleratorIndex, stream); - // GpuPrimeWorkLimiter.Release(); return order; } @@ -133,8 +131,8 @@ internal static partial class PrimeOrderCalculator int acceleratorIndex = gpu.AcceleratorIndex; // GpuPrimeWorkLimiter.Acquire(); - ArrayView1D specialMaxFactorsView = gpu.InputView; - ArrayView1D specialMaxResultView = gpu.OutputUlongView2; + ArrayView specialMaxFactorsView = gpu.InputView; + ArrayView specialMaxResultView = gpu.OutputUlongView2; var kernelLauncher = gpu.SpecialMaxKernelLauncher; var stream = AcceleratorStreamPool.Rent(acceleratorIndex); @@ -147,7 +145,6 @@ internal static partial class PrimeOrderCalculator specialMaxFactorsView, factorCount, divisorData.Modulus, - gpu.OutputUlongView, specialMaxResultView); ulong result = 0UL; diff --git a/PerfectNumbers.Core/PrimeTester.cs b/PerfectNumbers.Core/PrimeTester.cs index 0ef6df5..28c30d7 100644 --- a/PerfectNumbers.Core/PrimeTester.cs +++ b/PerfectNumbers.Core/PrimeTester.cs @@ -6,6 +6,249 @@ using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core; +public sealed class PrimeTesterByLastDigit +{ + public PrimeTesterByLastDigit() + { + } + + [ThreadStatic] + private static PrimeTesterByLastDigit? _tester; + + public static PrimeTesterByLastDigit Exclusive => _tester ??= new(); + + public static bool IsPrimeCpu(ulong n) + { + 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; + ulong[] smallPrimeDivisorsMul; + + 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; + default: + smallPrimeDivisors = DivisorGenerator.SmallPrimesLastNine; + smallPrimeDivisorsMul = DivisorGenerator.SmallPrimesPow2LastNine; + break; + } + + int smallPrimeDivisorsLength = smallPrimeDivisors.Length; + for (int i = 0; i < smallPrimeDivisorsLength; i++) + { + ulong value = smallPrimeDivisorsMul[i]; + if (value > n) + { + break; + } + + value = n % smallPrimeDivisors[i]; + if (value == 0) + { + result = false; + break; + } + } + } + } + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool IsPrimeGpu(PrimeOrderCalculatorAccelerator gpu, ulong n) + { + byte flag = 0; + var inputView = gpu.InputView; + var outputView = gpu.OutputByteView; + + int acceleratorIndex = gpu.AcceleratorIndex; + AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); + inputView.CopyFromCPU(stream, ref n, 1); + + var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher; + + kernelLauncher( + stream, + 1, + inputView, + gpu.DevicePrimesLastOne, + gpu.DevicePrimesLastSeven, + gpu.DevicePrimesLastThree, + gpu.DevicePrimesLastNine, + gpu.DevicePrimesPow2LastOne, + gpu.DevicePrimesPow2LastSeven, + gpu.DevicePrimesPow2LastThree, + gpu.DevicePrimesPow2LastNine, + outputView); + + outputView.CopyToCPU(stream, ref flag, 1); + stream.Synchronize(); + AcceleratorStreamPool.Return(acceleratorIndex, stream); + + return flag != 0; + } + + public static int GpuBatchSize { get; set; } = 262_144; + + private static readonly object GpuWarmUpLock = new(); + private static int WarmedGpuLeaseCount; + + public static void WarmUpGpuKernels(int threadCount) + { + int target = threadCount >> 2; + if (target == 0) + { + target = threadCount; + } + + lock (GpuWarmUpLock) + { + if (target <= WarmedGpuLeaseCount) + { + return; + } + + PrimeOrderCalculatorAccelerator.WarmUp(); + WarmedGpuLeaseCount = target; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void IsPrimeBatchGpu(PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan values, Span results) + { + // GpuPrimeWorkLimiter.Acquire(); + int acceleratorIndex = gpu.AcceleratorIndex; + int totalLength = values.Length; + int batchSize = GpuBatchSize; + + var inputView = gpu.InputView; + var outputView = gpu.OutputByteView; + + int pos = 0; + AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); + var kernelLauncher = gpu.SmallPrimeSieveKernelLauncher; + + while (pos < totalLength) + { + int remaining = totalLength - pos; + int count = remaining > batchSize ? batchSize : remaining; + + var valueSlice = values.Slice(pos, count); + inputView.CopyFromCPU(stream, valueSlice); + + kernelLauncher( + stream, + count, + inputView, + gpu.DevicePrimesLastOne, + gpu.DevicePrimesLastSeven, + gpu.DevicePrimesLastThree, + gpu.DevicePrimesLastNine, + gpu.DevicePrimesPow2LastOne, + gpu.DevicePrimesPow2LastSeven, + gpu.DevicePrimesPow2LastThree, + gpu.DevicePrimesPow2LastNine, + outputView); + + var resultSlice = results.Slice(pos, count); + outputView.CopyToCPU(stream, resultSlice); + + pos += count; + } + + stream.Synchronize(); + AcceleratorStreamPool.Return(acceleratorIndex, stream); + // GpuPrimeWorkLimiter.Release(); + } + + // Expose cache clearing for accelerator disposal coordination + public static void ClearGpuCaches() + { + PrimeOrderCalculatorAccelerator.Clear(); + } + + internal static void DisposeGpuContexts() + { + PrimeOrderCalculatorAccelerator.DisposeAll(); + lock (GpuWarmUpLock) + { + WarmedGpuLeaseCount = 0; + } + } + + [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 n.BinaryGcd(m) != 1UL; + } + + internal static void SharesFactorWithMaxExponentBatch(PrimeOrderCalculatorAccelerator gpu, 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. + // Check in benchmarks, which implementation was the fastest, is compatible with GPU, + // and implement it. + + int length = values.Length; + + int acceleratorIndex = gpu.AcceleratorIndex; + gpu.EnsureCapacity(0, length); + var inputBufferView = gpu.InputView; + var resultBufferView = gpu.OutputByteView; + var kernelLauncher = gpu.SharesFactorKernelLauncher; + + AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); + inputBufferView.CopyFromCPU(stream, values); + + kernelLauncher(stream, length, inputBufferView, resultBufferView); + + resultBufferView.CopyToCPU(stream, in results); + stream.Synchronize(); + + AcceleratorStreamPool.Return(acceleratorIndex, stream); + } +} + public sealed class PrimeTester { public PrimeTester() @@ -34,6 +277,54 @@ public sealed class PrimeTester 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; + 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 IsPrimeLastDigit(ulong n) + { + 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 diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 86c4c35..85187bf 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -340,23 +340,14 @@ public static partial class ULongExtensions // Benchmarks (Mod5ULongBenchmarks) show the direct `% 5` is still cheaper (~0.26 ns vs 0.43 ns), so keep the modulo until a faster lookup is proven. // (Mod8/Mod10 stay masked because they win; Mod5 currently does not.) [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong Mod10(this ulong value) => (value & 1UL) == 0UL - ? (value % 5UL) switch - { - 0UL => 0UL, - 1UL => 6UL, - 2UL => 2UL, - 3UL => 8UL, - _ => 4UL, - } - : (value % 5UL) switch - { - 0UL => 5UL, - 1UL => 1UL, - 2UL => 7UL, - 3UL => 3UL, - _ => 9UL, - }; + public static ulong Mod10(this ulong value) => (value % 5UL) switch + { + 0UL => 0UL, + 1UL => 6UL, + 2UL => 2UL, + 3UL => 8UL, + _ => 4UL, + }; [MethodImpl(MethodImplOptions.AggressiveInlining)] public static ulong Mod128(this ulong value) => value & 127UL;