From 131b3265bcf9bcf4546594c6634ee3de3a440721 Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Thu, 11 Dec 2025 00:00:12 +0100 Subject: [PATCH] PrimeOrderCalculator.Cpu cleaned up - hybrid & GPU related methods moved to their corresponding files. --- .../Cpu/MersenneCpuDivisorScanSession.cs | 21 +- ...ator.cs => PrimeOrderCalculator.Config.cs} | 0 .../PrimeOrderCalculator.Cpu.cs | 1024 +---------------- .../PrimeOrderCalculator.Gpu.cs | 948 ++++++++++----- .../PrimeOrderCalculator.Hybrid.cs | 612 ++++++++++ 5 files changed, 1304 insertions(+), 1301 deletions(-) rename PerfectNumbers.Core/{PrimeOrderCalculator.cs => PrimeOrderCalculator.Config.cs} (100%) diff --git a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs index a04b62b..c1cc3e7 100644 --- a/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs +++ b/PerfectNumbers.Core/Cpu/MersenneCpuDivisorScanSession.cs @@ -1,3 +1,4 @@ +using System.Runtime.CompilerServices; using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core.Cpu; @@ -55,21 +56,21 @@ internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelera // Keep these remainder steppers in place so future updates continue reusing the previously computed residues. // They are critical for avoiding repeated full Montgomery exponentiation work when scanning divisors. - var exponentStepper = new ExponentRemainderStepperCpu(divisorData, divisorCycle); - var cycleStepper = new CycleRemainderStepper(divisorCycle); + if (cycleStepper.Initialize(primes[0]) == 0UL) + { + return true; + } - bool initialUnity = exponentStepper.InitializeCpuIsUnity(primes[0]); - ulong remainder = cycleStepper.Initialize(primes[0]); - if (remainder == 0UL && initialUnity) + var exponentStepper = new ExponentRemainderStepperCpu(divisorData, divisorCycle); + if (exponentStepper.InitializeCpuIsUnity(primes[0])) { return true; } for (int i = 1; i < length; i++) { - remainder = cycleStepper.ComputeNext(primes[i]); - if (remainder != 0UL) + if (cycleStepper.ComputeNext(primes[i]) != 0UL) { continue; } @@ -83,8 +84,6 @@ internal sealed class MersenneCpuDivisorScanSession(PrimeOrderCalculatorAccelera return false; } - public void Return() - { - ThreadStaticPools.ReturnMersenneCpuDivisorSession(this); - } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public void Return() => ThreadStaticPools.ReturnMersenneCpuDivisorSession(this); } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Config.cs similarity index 100% rename from PerfectNumbers.Core/PrimeOrderCalculator.cs rename to PerfectNumbers.Core/PrimeOrderCalculator.Config.cs diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index c0cdd1e..64ab716 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -3,10 +3,8 @@ using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; -using ILGPU; using PerfectNumbers.Core.Cpu; using PerfectNumbers.Core.Gpu; -using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core; @@ -189,43 +187,6 @@ internal static partial class PrimeOrderCalculator return candidateOrder; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong RunHeuristicPipelineHybrid( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong? previousOrder, - in PrimeOrderCalculatorConfig config, - in MontgomeryDivisorData divisorData, - ulong phi, - PartialFactorResult phiFactors) - { - if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, prime, phiFactors, divisorData)) - { - return phi; - } - - ulong candidateOrder = InitializeStartingOrderHybrid(gpu, prime, phi, divisorData); - candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); - - if (TryConfirmOrderHybrid(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) - { - return candidateOrder; - } - - if (config.StrictMode) - { - orderFactors?.Dispose(); - return CalculateByFactorizationCpu(prime, divisorData); - } - - if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) - { - return order; - } - - return candidateOrder; - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static UInt128 RunHeuristicPipelineWideCpu( in UInt128 prime, @@ -550,28 +511,6 @@ internal static partial class PrimeOrderCalculator return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData); } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - int length = factors.Count; - // The phi factorization on the scanner path always yields at least one entry for phi >= 2, - // preventing the zero-length case from occurring outside synthetic tests. - // if (length == 0) - // { - // return true; - // } - - ReadOnlySpan factorSpan = new(factors.Factors, 0, length); - Span stackBuffer = stackalloc ulong[length]; - - if (length <= 8) - { - return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData); - } - - return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, prime, divisorData); - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static bool EvaluateSpecialMaxCandidatesCpu(Span buffer, ReadOnlySpan factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData) { @@ -633,22 +572,6 @@ internal static partial class PrimeOrderCalculator return order; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong InitializeStartingOrderHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong phi, in MontgomeryDivisorData divisorData) - { - ulong order = phi; - if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) - { - ulong half = phi >> 1; - if (Pow2EqualsOneHybrid(gpu, half, prime, divisorData)) - { - order = half; - } - } - - return order; - } - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) { @@ -818,55 +741,6 @@ internal static partial class PrimeOrderCalculator return true; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmOrderHybrid( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong order, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - out PartialFactorResult? reusableFactorization) - { - reusableFactorization = null; - - if (order == 0UL) - { - return false; - } - - if (!Pow2EqualsOneHybrid(gpu, order, prime, divisorData)) - { - return false; - } - - PartialFactorResult factorization = PartialFactorHybrid(gpu, order, config); - if (factorization.Factors is null) - { - factorization.Dispose(); - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) - { - factorization.Dispose(); - return false; - } - - factorization.WithAdditionalPrime(factorization.Cofactor); - } - - if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization)) - { - reusableFactorization = factorization; - return false; - } - - factorization.Dispose(); - return true; - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static bool ValidateOrderAgainstFactors( ulong prime, @@ -1100,338 +974,6 @@ internal static partial class PrimeOrderCalculator } } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryHeuristicFinishHybrid( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong order, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - PartialFactorResult phiFactors, - PartialFactorResult? cachedOrderFactors, - out ulong result) - { - result = 0UL; - if (order <= 1UL) - { - cachedOrderFactors?.Dispose(); - return false; - } - - // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, config); - try - { - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= 1UL || !orderFactors.CofactorIsPrime) - { - return false; - } - - orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); - candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); - int candidateCount = candidates.Count; - if (candidateCount == 0) - { - ThreadStaticPools.ReturnUlongList(candidates); - return false; - } - - SortCandidates(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; - int powUsed = 0; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool powStepperInitialized = false; - - int index = 0; - Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; - ArrayPool pool = ThreadStaticPools.UlongPool; - while (index < candidateCount && powUsed < powBudget) - { - int remaining = candidateCount - index; - int budgetRemaining = powBudget - powUsed; - int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize)); - if (batchSize <= 0) - { - break; - } - - ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); - int startIndex = 0; - - if (!powStepperInitialized) - { - ulong candidate = batch[0]; - powUsed++; - - bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); - powStepperInitialized = true; - if (equalsOne && TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - startIndex = 1; - } - - for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) - { - ulong candidate = batch[i]; - powUsed++; - - bool equalsOne = powStepper.ComputeNextIsUnity(candidate); - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - index += batchSize; - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - - return false; - } - finally - { - orderFactors.Dispose(); - } - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryHeuristicFinishGpu( - PrimeOrderCalculatorAccelerator gpu, - ulong prime, - ulong order, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderCalculatorConfig config, - PartialFactorResult phiFactors, - PartialFactorResult? cachedOrderFactors, - out ulong result) - { - result = 0UL; - if (order <= 1UL) - { - cachedOrderFactors?.Dispose(); - return false; - } - - // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, order, config); - try - { - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= 1UL || !orderFactors.CofactorIsPrime) - { - return false; - } - - orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); - candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); - int candidateCount = candidates.Count; - if (candidateCount == 0) - { - ThreadStaticPools.ReturnUlongList(candidates); - return false; - } - - SortCandidates(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; - int powUsed = 0; - bool allowGpuBatch = true; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool powStepperInitialized = false; - - int index = 0; - Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; - ArrayPool pool = ThreadStaticPools.UlongPool; - while (index < candidateCount && powUsed < powBudget) - { - int remaining = candidateCount - index; - int budgetRemaining = powBudget - powUsed; - int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize)); - if (batchSize <= 0) - { - break; - } - - ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); - ulong[]? gpuPool = null; - Span pooledGpuRemainders = default; - bool gpuSuccess = false; - bool gpuStackRemainders = false; - GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; - - if (allowGpuBatch) - { - if (batchSize <= PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount) - { - Span localRemainders = stackGpuRemainders[..batchSize]; - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, localRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - gpuSuccess = true; - gpuStackRemainders = true; - } - } - else - { - gpuPool = pool.Rent(batchSize); - Span pooledRemainders = gpuPool.AsSpan(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, pooledRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - pooledGpuRemainders = pooledRemainders; - gpuSuccess = true; - } - else - { - pool.Return(gpuPool, clearArray: false); - gpuPool = null; - } - } - - if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) - { - allowGpuBatch = false; - } - } - - int startIndex = 0; - - if (!gpuSuccess && !powStepperInitialized) - { - ulong candidate = batch[0]; - powUsed++; - - bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); - powStepperInitialized = true; - if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - startIndex = 1; - } - - for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) - { - ulong candidate = batch[i]; - powUsed++; - - bool equalsOne; - if (gpuSuccess) - { - ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; - equalsOne = remainderValue == 1UL; - } - else - { - equalsOne = powStepper.ComputeNextIsUnity(candidate); - } - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - result = candidate; - return true; - } - - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - index += batchSize; - } - - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - // DebugLog("No candidate confirmed"); - return false; - } - finally - { - orderFactors.Dispose(); - } - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) { @@ -1706,85 +1248,6 @@ internal static partial class PrimeOrderCalculator return !violates; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmCandidateHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) - { - PartialFactorResult factorization = PartialFactorHybrid(gpu, candidate, config); - - if (factorization.Factors is null) - { - factorization.Dispose(); - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) - { - factorization.Dispose(); - return false; - } - - factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan factorsSpan = factorization.Factors; - ReadOnlySpan exponentsSpan = factorization.Exponents; - int length = factorization.Count; - - const int StackExponentCapacity = 32; - const int ExponentHardLimit = 256; - - Span stackBuffer = stackalloc ulong[StackExponentCapacity]; - ulong[]? heapCandidateArray = null; - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool violates = false; - - for (int i = 0; i < length; i++) - { - ulong primeFactor = factorsSpan[i]; - int exponent = exponentsSpan[i]; - - // Factorization entries for candidate confirmation are always positive here. - // if (exponent <= 0) - // { - // continue; - // } - - if (exponent > ExponentHardLimit) - { - throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}."); - } - - if (exponent <= StackExponentCapacity) - { - if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) - { - violates = true; - break; - } - - continue; - } - - heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit); - if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) - { - violates = true; - break; - } - } - - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - if (heapCandidateArray is not null) - { - ThreadStaticPools.UlongPool.Return(heapCandidateArray); - } - factorization.Dispose(); - - return !violates; - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static bool CheckCandidateViolation(Span buffer, ulong primeFactor, int exponent, ulong candidate, ulong prime, ref int powUsed, int powBudget, ref ExponentRemainderStepperCpu stepper) { @@ -1871,28 +1334,6 @@ internal static partial class PrimeOrderCalculator private static bool Pow2EqualsOneCpu(ulong exponent, in MontgomeryDivisorData divisorData) => exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static bool Pow2EqualsOneHybrid(PrimeOrderCalculatorAccelerator gpu, ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) - { - // if (exponent >= PerfectNumberConstants.MaxQForDivisorCycles || prime >= PerfectNumberConstants.MaxQForDivisorCycles) - // { - // // Atomic.Add(ref _pow2GpuHits, 1); - // // Console.WriteLine($"pow2 GPU Hits {_pow2GpuHits}"); - // ulong remainder = exponent.Pow2MontgomeryModWindowedConvertGpu(gpu, divisorData); - // // GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); - // // if (status == GpuPow2ModStatus.Success) - // // { - // return remainder == 1UL; - // // } - // } - - // Atomic.Add(ref _pow2CpuHits, 1); - // Console.WriteLine($"pow2 CPU Hits {_pow2CpuHits}"); - return exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; - } - - private static int _cpuCount; - private static readonly ulong[] _emptyUlongArray = []; private static readonly int[] _emptyIntArray = []; @@ -1912,14 +1353,11 @@ internal static partial class PrimeOrderCalculator // primeSlots.Clear(); // exponentSlots.Clear(); - ArgumentOutOfRangeException.ThrowIfEqual(config.SmallFactorLimit, 0U); - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(PrimeOrderConstants.GpuSmallPrimeFactorSlots); PopulateSmallPrimeFactorsCpu( value, - limit, + config.SmallFactorLimit, primeSlots, exponentSlots, out int factorCount, @@ -2061,392 +1499,6 @@ internal static partial class PrimeOrderCalculator return result; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } - - // stackalloc is faster than pooling - Span primeSlots = stackalloc ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; - Span exponentSlots = stackalloc int[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; - 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(); - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong remaining = value; - - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); - - PopulateSmallPrimeFactorsCpu( - value, - limit, - primeSlots, - exponentSlots, - out int factorCount, - out remaining); - - bool limitReached; - if (remaining > 1UL) - { - if (config.PollardRhoMilliseconds > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - Stack compositeStack = ThreadStaticPools.RentUlongStack(4); - compositeStack.Push(remaining); - - CollectFactorsHybrid(gpu, primeSlots, exponentSlots, primeSlotArray: ref primeSlotArray, ref exponentSlotArray, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); - - if (limitReached) - { - while (compositeStack.Count > 0) - { - pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false)); - } - } - - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); - } - else - { - pending.Add(new PendingEntry(remaining, knownComposite: false)); - } - } - - ulong cofactor = 1UL; - bool cofactorContainsComposite = false; - int pendingCount = pending.Count; - int index = 0; - for (; index < pendingCount; index++) - { - PendingEntry entry = pending[index]; - ulong composite = entry.Value; - - if (entry.KnownComposite) - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - continue; - } - - if (!entry.HasKnownPrimality) - { - bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); - - entry.WithPrimality(isPrime); - pending[index] = entry; - } - - if (entry.IsPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); - } - // composite is never smaller on the execution path - // else if (composite > 1UL) - else - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - } - } - - bool cofactorIsPrime; - if (cofactor <= 1UL) - { - cofactorIsPrime = false; - } - else if (cofactorContainsComposite) - { - cofactorIsPrime = false; - } - else - { - cofactorIsPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(cofactor) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor); - } - - PartialFactorResult result; - if (factorCount == 0) - { - if (cofactor == value) - { - result = PartialFactorResult.Rent(cofactor, false, 0, cofactorIsPrime); - goto ReturnResult; - } - - result = PartialFactorResult.Rent(cofactor, cofactor == 1UL, 0, cofactorIsPrime); - goto ReturnResult; - } - - // This will never happen in production code. We'll always get at least 1 factor - // if (factorCount == 0) - // { - // if (cofactor == value) - // { - // result = PartialFactorResult.Rent(null, value, false, 0); - // goto ReturnResult; - // } - - // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - // goto ReturnResult; - // } - - ArrayPool ulongPool = ThreadStaticPools.UlongPool; - ArrayPool intPool = ThreadStaticPools.IntPool; - - ulong[] factorsArray = ulongPool.Rent(factorCount); - int[] exponentsArray = intPool.Rent(factorCount); - - for (int i = 0; i < factorCount; i++) - { - factorsArray[i] = primeSlots[i]; - exponentsArray[i] = exponentSlots[i]; - - // This will never happen on the execution path from production code - // if (primeValue == 0UL || exponentValue == 0) - // { - // throw new Exception("Prime value or exponent equals zero"); - // continue; - // } - } - - Array.Sort(factorsArray, exponentsArray, 0, factorCount); - - // Check if it was factored completely - limitReached = cofactor == 1UL; - - result = PartialFactorResult.Rent(factorsArray, exponentsArray, cofactor, limitReached, factorCount, cofactorIsPrime); - - ReturnResult: - if (primeSlotArray.Length > 0) - { - ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false); - ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false); - } - - pending.Clear(); - ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } - - // stackalloc is faster than pooling - Span primeSlots = stackalloc ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; - Span exponentSlots = stackalloc int[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; - 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(); - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong remaining = value; - - bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out int factorCount, out remaining); - - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); - - if (!gpuFactored) - { - PopulateSmallPrimeFactorsCpu( - value, - limit, - primeSlots, - exponentSlots, - out factorCount, - out remaining); - } - - bool limitReached; - if (remaining > 1UL) - { - if (config.PollardRhoMilliseconds > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - Stack compositeStack = ThreadStaticPools.RentUlongStack(4); - compositeStack.Push(remaining); - - CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); - - if (limitReached) - { - while (compositeStack.Count > 0) - { - pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false)); - } - } - - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); - } - else - { - pending.Add(new PendingEntry(remaining, knownComposite: false)); - } - } - - ulong cofactor = 1UL; - bool cofactorContainsComposite = false; - int pendingCount = pending.Count; - int index = 0; - for (; index < pendingCount; index++) - { - PendingEntry entry = pending[index]; - ulong composite = entry.Value; - - if (entry.KnownComposite) - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - continue; - } - - if (!entry.HasKnownPrimality) - { - bool isPrime = PrimeTesterByLastDigit.IsPrimeCpu(composite); - - entry.WithPrimality(isPrime); - pending[index] = entry; - } - - if (entry.IsPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); - } - // composite is never smaller on the execution path - // else if (composite > 1UL) - else - { - cofactor = checked(cofactor * composite); - cofactorContainsComposite = true; - } - } - - bool cofactorIsPrime; - if (cofactor <= 1UL) - { - cofactorIsPrime = false; - } - else if (cofactorContainsComposite) - { - cofactorIsPrime = false; - } - else - { - cofactorIsPrime = PrimeTesterByLastDigit.IsPrimeCpu(cofactor); - } - - PartialFactorResult result; - if (factorCount == 0) - { - if (cofactor == value) - { - result = PartialFactorResult.Rent(cofactor, false, 0, cofactorIsPrime); - goto ReturnResult; - } - - result = PartialFactorResult.Rent(cofactor, cofactor == 1UL, 0, cofactorIsPrime); - goto ReturnResult; - } - - // This will never happen in production code. We'll always get at least 1 factor - // if (factorCount == 0) - // { - // if (cofactor == value) - // { - // result = PartialFactorResult.Rent(null, value, false, 0); - // goto ReturnResult; - // } - - // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - // goto ReturnResult; - // } - - ArrayPool ulongPool = ThreadStaticPools.UlongPool; - ArrayPool intPool = ThreadStaticPools.IntPool; - - ulong[] factorsArray = ulongPool.Rent(factorCount); - int[] exponentsArray = intPool.Rent(factorCount); - - for (int i = 0; i < factorCount; i++) - { - factorsArray[i] = primeSlots[i]; - exponentsArray[i] = exponentSlots[i]; - - // This will never happen on the execution path from production code - // if (primeValue == 0UL || exponentValue == 0) - // { - // throw new Exception("Prime value or exponent equals zero"); - // continue; - // } - } - - Array.Sort(factorsArray, exponentsArray, 0, factorCount); - - // Check if it was factored completely - limitReached = cofactor == 1UL; - - result = PartialFactorResult.Rent(factorsArray, exponentsArray, cofactor, limitReached, factorCount, cofactorIsPrime); - - ReturnResult: - if (primeSlotArray.Length > 0) - { - ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false); - ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false); - } - - pending.Clear(); - ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); - return result; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span primeSlots, Span exponentSlots, ref ulong[] primeSlotArray, ref int[] exponentSlotArray, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) - { - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); - - if (isPrime) - { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); - continue; - } - - pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; - if (pollardRhoDeadlineReached) - { - pending.Add(new PendingEntry(composite, knownComposite: true)); - continue; - } - - 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); - } - - pollardRhoDeadlineReached = false; - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static void CollectFactorsCpu(Span primeSlots, Span exponentSlots, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) { @@ -2482,80 +1534,6 @@ internal static partial class PrimeOrderCalculator pollardRhoDeadlineReached = false; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - 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(primeSlots, exponentSlots, 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.AggressiveOptimization)] - private static bool RunOnCpu() - { - if (PerfectNumberConstants.GpuRatio == 0) - { - return false; - } - - if (PerfectNumberConstants.GpuRatio == 1) - { - while (true) - { - if (Atomic.CompareExchange(ref _cpuCount, 0, 1) == 0) - { - return true; - } - - if (Atomic.CompareExchange(ref _cpuCount, 1, 0) == 1) - { - return false; - } - - Thread.Yield(); - } - } - - int cpuCount = Atomic.Add(ref _cpuCount, 1); - if (cpuCount == PerfectNumberConstants.GpuRatio) - { - Atomic.Add(ref _cpuCount, -PerfectNumberConstants.GpuRatio); - } - else if (cpuCount > PerfectNumberConstants.GpuRatio) - { - cpuCount -= PerfectNumberConstants.GpuRatio; - } - - return cpuCount != PerfectNumberConstants.GpuRatio; - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static void PopulateSmallPrimeFactorsCpu( ulong value, diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index bf3ffc3..9ef5553 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -1,5 +1,7 @@ using System.Buffers; +using System.Diagnostics; using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; using ILGPU; using ILGPU.Runtime; using PerfectNumbers.Core.Gpu; @@ -11,6 +13,49 @@ internal static partial class PrimeOrderCalculator { private const int GpuSmallPrimeFactorSlots = 64; + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static ulong CalculateByFactorizationGpu(PrimeOrderCalculatorAccelerator gpu, ulong prime, in MontgomeryDivisorData divisorData) + { + ulong phi = prime - 1UL; + Dictionary counts = gpu.Pow2ModEntriesToTestOnHost; + counts.Clear(); + + FactorCompletelyCpu(phi, counts, false); + if (counts.Count == 0) + { + return phi; + } + + int entryCount = counts.Count; + Span> entries = stackalloc KeyValuePair[entryCount]; + + int stored = 0; + foreach (var entry in counts) + { + entries[stored++] = entry; + } + + entries.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); + + ulong order = phi; + var acceleratorIndex = gpu.AcceleratorIndex; + var kernel = gpu.CheckFactorsKernel; + + gpu.EnsureCapacity(entryCount, 1); + var pow2ModEntriesToTestOnDeviceView = gpu.Pow2ModEntriesToTestOnDeviceView; + + AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); + pow2ModEntriesToTestOnDeviceView.CopyFromCPU(stream, entries); + + kernel.Launch(stream, 1, entryCount, phi, pow2ModEntriesToTestOnDeviceView, divisorData.Modulus, divisorData.NPrime, divisorData.MontgomeryOne, divisorData.MontgomeryTwo, divisorData.MontgomeryTwoSquared, gpu.OutputUlongView); + + gpu.OutputUlongView.CopyToCPU(stream, ref order, 1); + stream.Synchronize(); + + AcceleratorStreamPool.Return(acceleratorIndex, stream); + return order; + } + public static ulong CalculateGpu( PrimeOrderCalculatorAccelerator gpu, ulong prime, @@ -115,112 +160,78 @@ internal static partial class PrimeOrderCalculator return result; } - private static UInt128 RunHeuristicPipelineWideGpu( - PrimeOrderCalculatorAccelerator gpu, - in UInt128 prime, - in UInt128? previousOrder, - in MontgomeryDivisorData? divisorData, - in PrimeOrderCalculatorConfig config, - in UInt128 phi, - in PartialFactorResult128 phiFactors) - { - if (phiFactors.FullyFactored && TrySpecialMaxWideGpu(gpu, phi, prime, divisorData, phiFactors)) - { - return phi; - } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + 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); - UInt128 candidateOrder = InitializeStartingOrderWideGpu(gpu, prime, phi, divisorData); - candidateOrder = ExponentLoweringWideGpu(gpu, candidateOrder, prime, divisorData, phiFactors); + if (isPrime) + { + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + continue; + } - if (TryConfirmOrderWideGpu(gpu, prime, candidateOrder, divisorData, config)) - { - return candidateOrder; - } + long currentTimestamp = Stopwatch.GetTimestamp(); + if (currentTimestamp > deadlineTimestamp) + { + pollardRhoDeadlineReached = true; + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } - if (config.StrictMode) - { - return FinishStrictlyWideGpu(gpu, prime, divisorData); - } + // if (!TryPollardRhoGpu(gpu, composite, out ulong factor)) + if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) + { + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } - if (TryHeuristicFinishWideGpu(gpu, prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) - { - return order; - } + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + } - return candidateOrder; - } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool EvaluateSpecialMaxCandidatesGpu( + PrimeOrderCalculatorAccelerator gpu, + ReadOnlySpan factors, + ulong phi, + ulong prime, + in MontgomeryDivisorData divisorData) + { + int factorCount = factors.Length; + gpu.EnsureUlongInputOutputCapacity(factorCount); - private static UInt128 FinishStrictlyWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in MontgomeryDivisorData? divisorData) - { - UInt128 phi = prime - UInt128.One; - Dictionary counts = new(capacity: 8); - FactorCompletelyWide(phi, counts); - if (counts.Count == 0) - { - return phi; - } + int acceleratorIndex = gpu.AcceleratorIndex; + // GpuPrimeWorkLimiter.Acquire(); + ArrayView specialMaxFactorsView = gpu.InputView; + ArrayView specialMaxResultView = gpu.OutputUlongView2; + var kernelLauncher = gpu.SpecialMaxKernelLauncher; - List> entries = new(counts); - entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); + var stream = AcceleratorStreamPool.Rent(acceleratorIndex); + specialMaxFactorsView.SubView(0, factorCount).CopyFromCPU(stream, factors); - UInt128 order = phi; - int entryCount = entries.Count; - for (int i = 0; i < entryCount; i++) - { - UInt128 primeFactor = entries[i].Key; - int exponent = entries[i].Value; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != UInt128.Zero) - { - UInt128 candidate = order / primeFactor; - if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) == UInt128.One) - { - order = candidate; - } - else - { - break; - } - } - } - } + kernelLauncher( + stream, + 1, + phi, + specialMaxFactorsView, + factorCount, + divisorData.Modulus, + specialMaxResultView); - return order; - } - - private static bool TrySpecialMaxWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) - { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - for (int i = 0; i < length; i++) - { - UInt128 factor = factorSpan[i].Value; - UInt128 reduced = phi / factor; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - return false; - } - } - - return true; - } - - private static UInt128 InitializeStartingOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData? divisorData) - { - UInt128 order = phi; - UInt128 mod8 = prime & (UInt128)7UL; - if (mod8 == UInt128.One || mod8 == (UInt128)7UL) - { - UInt128 half = phi >> 1; - if (Pow2ModWideGpu(gpu, half, prime, divisorData) == UInt128.One) - { - order = half; - } - } - - return order; - } + ulong result = 0UL; + specialMaxResultView.CopyToCPU(stream, ref result, 1); + stream.Synchronize(); + + AcceleratorStreamPool.Return(acceleratorIndex, stream); + // GpuPrimeWorkLimiter.Release(); + return result != 0; + } private static UInt128 ExponentLoweringWideGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 order, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) { @@ -270,6 +281,355 @@ internal static partial class PrimeOrderCalculator return order; } + private static UInt128 FinishStrictlyWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in MontgomeryDivisorData? divisorData) + { + UInt128 phi = prime - UInt128.One; + Dictionary counts = new(capacity: 8); + FactorCompletelyWide(phi, counts); + if (counts.Count == 0) + { + return phi; + } + + List> entries = new(counts); + entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); + + UInt128 order = phi; + int entryCount = entries.Count; + for (int i = 0; i < entryCount; i++) + { + UInt128 primeFactor = entries[i].Key; + int exponent = entries[i].Value; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != UInt128.Zero) + { + UInt128 candidate = order / primeFactor; + if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) == UInt128.One) + { + order = candidate; + } + else + { + break; + } + } + } + } + + return order; + } + + private static UInt128 InitializeStartingOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData? divisorData) + { + UInt128 order = phi; + UInt128 mod8 = prime & (UInt128)7UL; + if (mod8 == UInt128.One || mod8 == (UInt128)7UL) + { + UInt128 half = phi >> 1; + if (Pow2ModWideGpu(gpu, half, prime, divisorData) == UInt128.One) + { + order = half; + } + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } + + // stackalloc is faster than pooling + Span primeSlots = stackalloc ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; + Span exponentSlots = stackalloc int[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; + 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(); + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong remaining = value; + + bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out int factorCount, out remaining); + + List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); + + if (!gpuFactored) + { + PopulateSmallPrimeFactorsCpu( + value, + limit, + primeSlots, + exponentSlots, + out factorCount, + out remaining); + } + + bool limitReached; + if (remaining > 1UL) + { + if (config.PollardRhoMilliseconds > 0) + { + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + Stack compositeStack = ThreadStaticPools.RentUlongStack(4); + compositeStack.Push(remaining); + + CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); + + if (limitReached) + { + while (compositeStack.Count > 0) + { + pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false)); + } + } + + compositeStack.Clear(); + ThreadStaticPools.ReturnUlongStack(compositeStack); + } + else + { + pending.Add(new PendingEntry(remaining, knownComposite: false)); + } + } + + ulong cofactor = 1UL; + bool cofactorContainsComposite = false; + int pendingCount = pending.Count; + int index = 0; + for (; index < pendingCount; index++) + { + PendingEntry entry = pending[index]; + ulong composite = entry.Value; + + if (entry.KnownComposite) + { + cofactor = checked(cofactor * composite); + cofactorContainsComposite = true; + continue; + } + + if (!entry.HasKnownPrimality) + { + bool isPrime = PrimeTesterByLastDigit.IsPrimeCpu(composite); + + entry.WithPrimality(isPrime); + pending[index] = entry; + } + + if (entry.IsPrime) + { + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + } + // composite is never smaller on the execution path + // else if (composite > 1UL) + else + { + cofactor = checked(cofactor * composite); + cofactorContainsComposite = true; + } + } + + bool cofactorIsPrime; + if (cofactor <= 1UL) + { + cofactorIsPrime = false; + } + else if (cofactorContainsComposite) + { + cofactorIsPrime = false; + } + else + { + cofactorIsPrime = PrimeTesterByLastDigit.IsPrimeCpu(cofactor); + } + + PartialFactorResult result; + if (factorCount == 0) + { + if (cofactor == value) + { + result = PartialFactorResult.Rent(cofactor, false, 0, cofactorIsPrime); + goto ReturnResult; + } + + result = PartialFactorResult.Rent(cofactor, cofactor == 1UL, 0, cofactorIsPrime); + goto ReturnResult; + } + + // This will never happen in production code. We'll always get at least 1 factor + // if (factorCount == 0) + // { + // if (cofactor == value) + // { + // result = PartialFactorResult.Rent(null, value, false, 0); + // goto ReturnResult; + // } + + // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + // goto ReturnResult; + // } + + ArrayPool ulongPool = ThreadStaticPools.UlongPool; + ArrayPool intPool = ThreadStaticPools.IntPool; + + ulong[] factorsArray = ulongPool.Rent(factorCount); + int[] exponentsArray = intPool.Rent(factorCount); + + for (int i = 0; i < factorCount; i++) + { + factorsArray[i] = primeSlots[i]; + exponentsArray[i] = exponentSlots[i]; + + // This will never happen on the execution path from production code + // if (primeValue == 0UL || exponentValue == 0) + // { + // throw new Exception("Prime value or exponent equals zero"); + // continue; + // } + } + + Array.Sort(factorsArray, exponentsArray, 0, factorCount); + + // Check if it was factored completely + limitReached = cofactor == 1UL; + + result = PartialFactorResult.Rent(factorsArray, exponentsArray, cofactor, limitReached, factorCount, cofactorIsPrime); + + ReturnResult: + if (primeSlotArray.Length > 0) + { + ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false); + ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false); + } + + pending.Clear(); + ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static UInt128 Pow2ModWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData? divisorData) + { + if (modulus == UInt128.One) + { + return UInt128.Zero; + } + + if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue) + { + ulong prime64 = (ulong)modulus; + ulong exponent64 = (ulong)exponent; + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent64, prime64, out ulong remainder, divisorData); + if (status == GpuPow2ModStatus.Success) + { + return remainder; + } + } + + GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent, modulus, out UInt128 wideRemainder); + if (wideStatus == GpuPow2ModStatus.Success) + { + return wideRemainder; + } + + return exponent.Pow2MontgomeryModWindowed(modulus); + } + + private static UInt128 RunHeuristicPipelineWideGpu( + PrimeOrderCalculatorAccelerator gpu, + in UInt128 prime, + in UInt128? previousOrder, + in MontgomeryDivisorData? divisorData, + in PrimeOrderCalculatorConfig config, + in UInt128 phi, + in PartialFactorResult128 phiFactors) + { + if (phiFactors.FullyFactored && TrySpecialMaxWideGpu(gpu, phi, prime, divisorData, phiFactors)) + { + return phi; + } + + UInt128 candidateOrder = InitializeStartingOrderWideGpu(gpu, prime, phi, divisorData); + candidateOrder = ExponentLoweringWideGpu(gpu, candidateOrder, prime, divisorData, phiFactors); + + if (TryConfirmOrderWideGpu(gpu, prime, candidateOrder, divisorData, config)) + { + return candidateOrder; + } + + if (config.StrictMode) + { + return FinishStrictlyWideGpu(gpu, prime, divisorData); + } + + if (TryHeuristicFinishWideGpu(gpu, prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) + { + return order; + } + + return candidateOrder; + } + + private static bool TryConfirmCandidateWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData? divisorData) + { + PartialFactorResult128 factorization = PartialFactorWide(candidate, config); + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= UInt128.One) + { + return false; + } + + if (!IsPrimeWide(factorization.Cofactor)) + { + return false; + } + + factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + } + + ReadOnlySpan span = factorization.Factors; + int length = factorization.Count; + for (int i = 0; i < length; i++) + { + UInt128 primeFactor = span[i].Value; + UInt128 reduced = candidate; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + if ((reduced % primeFactor) != UInt128.Zero) + { + break; + } + + reduced /= primeFactor; + if (powUsed >= powBudget && powBudget > 0) + { + return false; + } + + powUsed++; + if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) + { + return false; + } + } + } + + return true; + } + private static bool TryConfirmOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 order, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) { if (order == UInt128.Zero) @@ -327,6 +687,205 @@ internal static partial class PrimeOrderCalculator return true; } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TryHeuristicFinishGpu( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + ulong order, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + PartialFactorResult phiFactors, + PartialFactorResult? cachedOrderFactors, + out ulong result) + { + result = 0UL; + if (order <= 1UL) + { + cachedOrderFactors?.Dispose(); + return false; + } + + // Reuse the partial factorization from TryConfirmOrderCpu when available. + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, order, config); + try + { + if (orderFactors.Factors is null) + { + return false; + } + + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= 1UL || !orderFactors.CofactorIsPrime) + { + return false; + } + + orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; + List candidates = ThreadStaticPools.RentUlongList(capacity); + candidates.Clear(); + ulong[] factorArray = orderFactors.Factors!; + int[] exponentsArray = orderFactors.Exponents!; + + BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); + int candidateCount = candidates.Count; + if (candidateCount == 0) + { + ThreadStaticPools.ReturnUlongList(candidates); + return false; + } + + SortCandidates(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; + int powUsed = 0; + bool allowGpuBatch = true; + Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + bool powStepperInitialized = false; + + int index = 0; + Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; + ArrayPool pool = ThreadStaticPools.UlongPool; + while (index < candidateCount && powUsed < powBudget) + { + int remaining = candidateCount - index; + int budgetRemaining = powBudget - powUsed; + int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize)); + if (batchSize <= 0) + { + break; + } + + ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); + ulong[]? gpuPool = null; + Span pooledGpuRemainders = default; + bool gpuSuccess = false; + bool gpuStackRemainders = false; + GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; + + if (allowGpuBatch) + { + if (batchSize <= PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount) + { + Span localRemainders = stackGpuRemainders[..batchSize]; + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, localRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) + { + gpuSuccess = true; + gpuStackRemainders = true; + } + } + else + { + gpuPool = pool.Rent(batchSize); + Span pooledRemainders = gpuPool.AsSpan(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(gpu, batch, prime, pooledRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) + { + pooledGpuRemainders = pooledRemainders; + gpuSuccess = true; + } + else + { + pool.Return(gpuPool, clearArray: false); + gpuPool = null; + } + } + + if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) + { + allowGpuBatch = false; + } + } + + int startIndex = 0; + + if (!gpuSuccess && !powStepperInitialized) + { + ulong candidate = batch[0]; + powUsed++; + + bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); + powStepperInitialized = true; + if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + result = candidate; + return true; + } + + startIndex = 1; + } + + for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) + { + ulong candidate = batch[i]; + powUsed++; + + bool equalsOne; + if (gpuSuccess) + { + ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; + equalsOne = remainderValue == 1UL; + } + else + { + equalsOne = powStepper.ComputeNextIsUnity(candidate); + } + if (!equalsOne) + { + continue; + } + + if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + continue; + } + + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + result = candidate; + return true; + } + + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + index += batchSize; + } + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + // DebugLog("No candidate confirmed"); + return false; + } + finally + { + orderFactors.Dispose(); + } + } + private static bool TryHeuristicFinishWideGpu( PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, @@ -405,47 +964,33 @@ internal static partial class PrimeOrderCalculator return false; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong CalculateByFactorizationGpu(PrimeOrderCalculatorAccelerator gpu, ulong prime, in MontgomeryDivisorData divisorData) + private static bool TryPollardRhoGpu(PrimeOrderCalculatorAccelerator gpu, ulong n, out ulong factor) { - ulong phi = prime - 1UL; - Dictionary counts = gpu.Pow2ModEntriesToTestOnHost; - counts.Clear(); + Span randomStateSpan = stackalloc ulong[1]; + randomStateSpan[0] = ThreadStaticDeterministicRandomGpu.Exclusive.State; + + Span factoredSpan = stackalloc byte[1]; + Span factorSpan = stackalloc ulong[1]; - FactorCompletelyCpu(phi, counts, false); - if (counts.Count == 0) - { - return phi; - } - - int entryCount = counts.Count; - Span> entries = stackalloc KeyValuePair[entryCount]; - - int stored = 0; - foreach (var entry in counts) - { - entries[stored++] = entry; - } - - entries.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); - - ulong order = phi; - var acceleratorIndex = gpu.AcceleratorIndex; - var kernel = gpu.CheckFactorsKernel; - - gpu.EnsureCapacity(entryCount, 1); - var pow2ModEntriesToTestOnDeviceView = gpu.Pow2ModEntriesToTestOnDeviceView; + int acceleratorIndex = gpu.AcceleratorIndex; + var kernel = gpu.PollardRhoKernel; AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); - pow2ModEntriesToTestOnDeviceView.CopyFromCPU(stream, entries); + gpu.InputView.CopyFromCPU(stream, randomStateSpan); - kernel.Launch(stream, 1, entryCount, phi, pow2ModEntriesToTestOnDeviceView, divisorData.Modulus, divisorData.NPrime, divisorData.MontgomeryOne, divisorData.MontgomeryTwo, divisorData.MontgomeryTwoSquared, gpu.OutputUlongView); + kernel.Launch(stream, 1, n, 1, gpu.InputView, gpu.OutputByteView, gpu.OutputUlongView); - gpu.OutputUlongView.CopyToCPU(stream, ref order, 1); + gpu.OutputByteView.CopyToCPU(stream, factoredSpan); + gpu.OutputUlongView.CopyToCPU(stream, factorSpan); + gpu.InputView.CopyToCPU(stream, randomStateSpan); stream.Synchronize(); AcceleratorStreamPool.Return(acceleratorIndex, stream); - return order; + ThreadStaticDeterministicRandomGpu.Exclusive.SetState(randomStateSpan[0]); + + bool factored = factoredSpan[0] != 0; + factor = factored ? factorSpan[0] : 0UL; + return factored; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] @@ -515,152 +1060,21 @@ internal static partial class PrimeOrderCalculator return true; } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool EvaluateSpecialMaxCandidatesGpu( - PrimeOrderCalculatorAccelerator gpu, - ReadOnlySpan factors, - ulong phi, - ulong prime, - in MontgomeryDivisorData divisorData) - { - int factorCount = factors.Length; - gpu.EnsureUlongInputOutputCapacity(factorCount); - - int acceleratorIndex = gpu.AcceleratorIndex; - // GpuPrimeWorkLimiter.Acquire(); - ArrayView specialMaxFactorsView = gpu.InputView; - ArrayView specialMaxResultView = gpu.OutputUlongView2; - var kernelLauncher = gpu.SpecialMaxKernelLauncher; - - var stream = AcceleratorStreamPool.Rent(acceleratorIndex); - specialMaxFactorsView.SubView(0, factorCount).CopyFromCPU(stream, factors); - - kernelLauncher( - stream, - 1, - phi, - specialMaxFactorsView, - factorCount, - divisorData.Modulus, - specialMaxResultView); - - ulong result = 0UL; - specialMaxResultView.CopyToCPU(stream, ref result, 1); - stream.Synchronize(); - - AcceleratorStreamPool.Return(acceleratorIndex, stream); - // GpuPrimeWorkLimiter.Release(); - return result != 0; - } - - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 Pow2ModWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData? divisorData) + private static bool TrySpecialMaxWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) { - if (modulus == UInt128.One) - { - return UInt128.Zero; - } - - if (modulus <= (UInt128)ulong.MaxValue && exponent <= (UInt128)ulong.MaxValue) - { - ulong prime64 = (ulong)modulus; - ulong exponent64 = (ulong)exponent; - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent64, prime64, out ulong remainder, divisorData); - if (status == GpuPow2ModStatus.Success) - { - return remainder; - } - } - - GpuPow2ModStatus wideStatus = PrimeOrderGpuHeuristics.TryPow2Mod(gpu, exponent, modulus, out UInt128 wideRemainder); - if (wideStatus == GpuPow2ModStatus.Success) - { - return wideRemainder; - } - - return exponent.Pow2MontgomeryModWindowed(modulus); - } - - private static bool TryConfirmCandidateWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData? divisorData) - { - PartialFactorResult128 factorization = PartialFactorWide(candidate, config); - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= UInt128.One) - { - return false; - } - - if (!IsPrimeWide(factorization.Cofactor)) - { - return false; - } - - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); - } - - ReadOnlySpan span = factorization.Factors; - int length = factorization.Count; + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; for (int i = 0; i < length; i++) { - UInt128 primeFactor = span[i].Value; - UInt128 reduced = candidate; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) + UInt128 factor = factorSpan[i].Value; + UInt128 reduced = phi / factor; + if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) { - if ((reduced % primeFactor) != UInt128.Zero) - { - break; - } - - reduced /= primeFactor; - if (powUsed >= powBudget && powBudget > 0) - { - return false; - } - - powUsed++; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - return false; - } + return false; } } return true; } - - private static bool TryPollardRhoGpu(PrimeOrderCalculatorAccelerator gpu, ulong n, out ulong factor) - { - Span randomStateSpan = stackalloc ulong[1]; - randomStateSpan[0] = ThreadStaticDeterministicRandomGpu.Exclusive.State; - - Span factoredSpan = stackalloc byte[1]; - Span factorSpan = stackalloc ulong[1]; - - int acceleratorIndex = gpu.AcceleratorIndex; - var kernel = gpu.PollardRhoKernel; - - AcceleratorStream stream = AcceleratorStreamPool.Rent(acceleratorIndex); - gpu.InputView.CopyFromCPU(stream, randomStateSpan); - - kernel.Launch(stream, 1, n, 1, gpu.InputView, gpu.OutputByteView, gpu.OutputUlongView); - - gpu.OutputByteView.CopyToCPU(stream, factoredSpan); - gpu.OutputUlongView.CopyToCPU(stream, factorSpan); - gpu.InputView.CopyToCPU(stream, randomStateSpan); - stream.Synchronize(); - - AcceleratorStreamPool.Return(acceleratorIndex, stream); - ThreadStaticDeterministicRandomGpu.Exclusive.SetState(randomStateSpan[0]); - - bool factored = factoredSpan[0] != 0; - factor = factored ? factorSpan[0] : 0UL; - return factored; - } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs index 26c3528..0e04b77 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs @@ -1,3 +1,9 @@ +using System.Buffers; +using System.Diagnostics; +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; +using ILGPU; +using PerfectNumbers.Core.Gpu; using PerfectNumbers.Core.Gpu.Accelerators; namespace PerfectNumbers.Core; @@ -75,4 +81,610 @@ internal static partial class PrimeOrderCalculator return result; } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span primeSlots, Span exponentSlots, ref ulong[] primeSlotArray, ref int[] exponentSlotArray, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) + { + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); + + if (isPrime) + { + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + continue; + } + + pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; + if (pollardRhoDeadlineReached) + { + pending.Add(new PendingEntry(composite, knownComposite: true)); + continue; + } + + 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); + } + + pollardRhoDeadlineReached = false; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static ulong InitializeStartingOrderHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong phi, in MontgomeryDivisorData divisorData) + { + ulong order = phi; + if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) + { + ulong half = phi >> 1; + if (Pow2EqualsOneHybrid(gpu, half, prime, divisorData)) + { + order = half; + } + } + + return order; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } + + // stackalloc is faster than pooling + Span primeSlots = stackalloc ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; + Span exponentSlots = stackalloc int[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; + 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(); + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong remaining = value; + + List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); + + PopulateSmallPrimeFactorsCpu( + value, + limit, + primeSlots, + exponentSlots, + out int factorCount, + out remaining); + + bool limitReached; + if (remaining > 1UL) + { + if (config.PollardRhoMilliseconds > 0) + { + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + Stack compositeStack = ThreadStaticPools.RentUlongStack(4); + compositeStack.Push(remaining); + + CollectFactorsHybrid(gpu, primeSlots, exponentSlots, primeSlotArray: ref primeSlotArray, ref exponentSlotArray, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); + + if (limitReached) + { + while (compositeStack.Count > 0) + { + pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false)); + } + } + + compositeStack.Clear(); + ThreadStaticPools.ReturnUlongStack(compositeStack); + } + else + { + pending.Add(new PendingEntry(remaining, knownComposite: false)); + } + } + + ulong cofactor = 1UL; + bool cofactorContainsComposite = false; + int pendingCount = pending.Count; + int index = 0; + for (; index < pendingCount; index++) + { + PendingEntry entry = pending[index]; + ulong composite = entry.Value; + + if (entry.KnownComposite) + { + cofactor = checked(cofactor * composite); + cofactorContainsComposite = true; + continue; + } + + if (!entry.HasKnownPrimality) + { + bool isPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(composite) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, composite); + + entry.WithPrimality(isPrime); + pending[index] = entry; + } + + if (entry.IsPrime) + { + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + } + // composite is never smaller on the execution path + // else if (composite > 1UL) + else + { + cofactor = checked(cofactor * composite); + cofactorContainsComposite = true; + } + } + + bool cofactorIsPrime; + if (cofactor <= 1UL) + { + cofactorIsPrime = false; + } + else if (cofactorContainsComposite) + { + cofactorIsPrime = false; + } + else + { + cofactorIsPrime = RunOnCpu() ? PrimeTesterByLastDigit.IsPrimeCpu(cofactor) : HeuristicCombinedPrimeTester.IsPrimeGpu(gpu, cofactor); + } + + PartialFactorResult result; + if (factorCount == 0) + { + if (cofactor == value) + { + result = PartialFactorResult.Rent(cofactor, false, 0, cofactorIsPrime); + goto ReturnResult; + } + + result = PartialFactorResult.Rent(cofactor, cofactor == 1UL, 0, cofactorIsPrime); + goto ReturnResult; + } + + // This will never happen in production code. We'll always get at least 1 factor + // if (factorCount == 0) + // { + // if (cofactor == value) + // { + // result = PartialFactorResult.Rent(null, value, false, 0); + // goto ReturnResult; + // } + + // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + // goto ReturnResult; + // } + + ArrayPool ulongPool = ThreadStaticPools.UlongPool; + ArrayPool intPool = ThreadStaticPools.IntPool; + + ulong[] factorsArray = ulongPool.Rent(factorCount); + int[] exponentsArray = intPool.Rent(factorCount); + + for (int i = 0; i < factorCount; i++) + { + factorsArray[i] = primeSlots[i]; + exponentsArray[i] = exponentSlots[i]; + + // This will never happen on the execution path from production code + // if (primeValue == 0UL || exponentValue == 0) + // { + // throw new Exception("Prime value or exponent equals zero"); + // continue; + // } + } + + Array.Sort(factorsArray, exponentsArray, 0, factorCount); + + // Check if it was factored completely + limitReached = cofactor == 1UL; + + result = PartialFactorResult.Rent(factorsArray, exponentsArray, cofactor, limitReached, factorCount, cofactorIsPrime); + + ReturnResult: + if (primeSlotArray.Length > 0) + { + ThreadStaticPools.UlongPool.Return(primeSlotArray, clearArray: false); + ThreadStaticPools.IntPool.Return(exponentSlotArray!, clearArray: false); + } + + pending.Clear(); + ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static bool Pow2EqualsOneHybrid(PrimeOrderCalculatorAccelerator gpu, ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) + { + // if (exponent >= PerfectNumberConstants.MaxQForDivisorCycles || prime >= PerfectNumberConstants.MaxQForDivisorCycles) + // { + // // Atomic.Add(ref _pow2GpuHits, 1); + // // Console.WriteLine($"pow2 GPU Hits {_pow2GpuHits}"); + // ulong remainder = exponent.Pow2MontgomeryModWindowedConvertGpu(gpu, divisorData); + // // GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); + // // if (status == GpuPow2ModStatus.Success) + // // { + // return remainder == 1UL; + // // } + // } + + // Atomic.Add(ref _pow2CpuHits, 1); + // Console.WriteLine($"pow2 CPU Hits {_pow2CpuHits}"); + return exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static ulong RunHeuristicPipelineHybrid( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + ulong? previousOrder, + in PrimeOrderCalculatorConfig config, + in MontgomeryDivisorData divisorData, + ulong phi, + PartialFactorResult phiFactors) + { + if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, prime, phiFactors, divisorData)) + { + return phi; + } + + ulong candidateOrder = InitializeStartingOrderHybrid(gpu, prime, phi, divisorData); + candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); + + if (TryConfirmOrderHybrid(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) + { + return candidateOrder; + } + + if (config.StrictMode) + { + orderFactors?.Dispose(); + return CalculateByFactorizationCpu(prime, divisorData); + } + + if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) + { + return order; + } + + return candidateOrder; + } + + private static int _cpuCount; + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool RunOnCpu() + { + if (PerfectNumberConstants.GpuRatio == 0) + { + return false; + } + + if (PerfectNumberConstants.GpuRatio == 1) + { + while (true) + { + if (Atomic.CompareExchange(ref _cpuCount, 0, 1) == 0) + { + return true; + } + + if (Atomic.CompareExchange(ref _cpuCount, 1, 0) == 1) + { + return false; + } + + Thread.Yield(); + } + } + + int cpuCount = Atomic.Add(ref _cpuCount, 1); + if (cpuCount == PerfectNumberConstants.GpuRatio) + { + Atomic.Add(ref _cpuCount, -PerfectNumberConstants.GpuRatio); + } + else if (cpuCount > PerfectNumberConstants.GpuRatio) + { + cpuCount -= PerfectNumberConstants.GpuRatio; + } + + return cpuCount != PerfectNumberConstants.GpuRatio; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TryConfirmCandidateHybrid(PrimeOrderCalculatorAccelerator gpu, ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) + { + PartialFactorResult factorization = PartialFactorHybrid(gpu, candidate, config); + + if (factorization.Factors is null) + { + factorization.Dispose(); + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) + { + factorization.Dispose(); + return false; + } + + factorization.WithAdditionalPrime(factorization.Cofactor); + } + + ReadOnlySpan factorsSpan = factorization.Factors; + ReadOnlySpan exponentsSpan = factorization.Exponents; + int length = factorization.Count; + + const int StackExponentCapacity = 32; + const int ExponentHardLimit = 256; + + Span stackBuffer = stackalloc ulong[StackExponentCapacity]; + ulong[]? heapCandidateArray = null; + ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + bool violates = false; + + for (int i = 0; i < length; i++) + { + ulong primeFactor = factorsSpan[i]; + int exponent = exponentsSpan[i]; + + // Factorization entries for candidate confirmation are always positive here. + // if (exponent <= 0) + // { + // continue; + // } + + if (exponent > ExponentHardLimit) + { + throw new InvalidOperationException($"Candidate factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}."); + } + + if (exponent <= StackExponentCapacity) + { + if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) + { + violates = true; + break; + } + + continue; + } + + heapCandidateArray ??= ThreadStaticPools.UlongPool.Rent(ExponentHardLimit); + if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) + { + violates = true; + break; + } + } + + ThreadStaticPools.ReturnExponentStepperCpu(stepper); + if (heapCandidateArray is not null) + { + ThreadStaticPools.UlongPool.Return(heapCandidateArray); + } + factorization.Dispose(); + + return !violates; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TryConfirmOrderHybrid( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + ulong order, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + out PartialFactorResult? reusableFactorization) + { + reusableFactorization = null; + + if (order == 0UL) + { + return false; + } + + if (!Pow2EqualsOneHybrid(gpu, order, prime, divisorData)) + { + return false; + } + + PartialFactorResult factorization = PartialFactorHybrid(gpu, order, config); + if (factorization.Factors is null) + { + factorization.Dispose(); + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL || !factorization.CofactorIsPrime) + { + factorization.Dispose(); + return false; + } + + factorization.WithAdditionalPrime(factorization.Cofactor); + } + + if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization)) + { + reusableFactorization = factorization; + return false; + } + + factorization.Dispose(); + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TryHeuristicFinishHybrid( + PrimeOrderCalculatorAccelerator gpu, + ulong prime, + ulong order, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderCalculatorConfig config, + PartialFactorResult phiFactors, + PartialFactorResult? cachedOrderFactors, + out ulong result) + { + result = 0UL; + if (order <= 1UL) + { + cachedOrderFactors?.Dispose(); + return false; + } + + // Reuse the partial factorization from TryConfirmOrderCpu when available. + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, config); + try + { + if (orderFactors.Factors is null) + { + return false; + } + + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= 1UL || !orderFactors.CofactorIsPrime) + { + return false; + } + + orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; + List candidates = ThreadStaticPools.RentUlongList(capacity); + candidates.Clear(); + ulong[] factorArray = orderFactors.Factors!; + int[] exponentsArray = orderFactors.Exponents!; + + BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); + int candidateCount = candidates.Count; + if (candidateCount == 0) + { + ThreadStaticPools.ReturnUlongList(candidates); + return false; + } + + SortCandidates(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; + int powUsed = 0; + Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + bool powStepperInitialized = false; + + int index = 0; + Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; + ArrayPool pool = ThreadStaticPools.UlongPool; + while (index < candidateCount && powUsed < powBudget) + { + int remaining = candidateCount - index; + int budgetRemaining = powBudget - powUsed; + int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, PrimeOrderConstants.MaxGpuBatchSize)); + if (batchSize <= 0) + { + break; + } + + ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); + int startIndex = 0; + + if (!powStepperInitialized) + { + ulong candidate = batch[0]; + powUsed++; + + bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); + powStepperInitialized = true; + if (equalsOne && TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + result = candidate; + return true; + } + + startIndex = 1; + } + + for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) + { + ulong candidate = batch[i]; + powUsed++; + + bool equalsOne = powStepper.ComputeNextIsUnity(candidate); + if (!equalsOne) + { + continue; + } + + if (!TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + continue; + } + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + result = candidate; + return true; + } + + index += batchSize; + } + + candidates.Clear(); + ThreadStaticPools.ReturnUlongList(candidates); + ThreadStaticPools.ReturnExponentStepperCpu(powStepper); + + return false; + } + finally + { + orderFactors.Dispose(); + } + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + { + int length = factors.Count; + // The phi factorization on the scanner path always yields at least one entry for phi >= 2, + // preventing the zero-length case from occurring outside synthetic tests. + // if (length == 0) + // { + // return true; + // } + + ReadOnlySpan factorSpan = new(factors.Factors, 0, length); + Span stackBuffer = stackalloc ulong[length]; + + if (length <= 8) + { + return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData); + } + + return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, prime, divisorData); + } + }