diff --git a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs index 3e875b6..c8c5bf6 100644 --- a/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs +++ b/PerfectNumbers.Core/ExponentRemainderStepperCpu.cs @@ -7,7 +7,7 @@ namespace PerfectNumbers.Core; /// Tracks successive Montgomery residues for a single divisor while scanning exponent candidates on the CPU path. /// Callers are responsible for seeding the stepper once per divisor and then advancing exponents in ascending order. /// -internal struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL) +public struct ExponentRemainderStepperCpu(in MontgomeryDivisorData divisor, ulong cycleLength = 0UL) { private readonly MontgomeryDivisorData _divisor = divisor; private readonly ulong _modulus = divisor.Modulus; diff --git a/PerfectNumbers.Core/PartialFactorPendingEntry.cs b/PerfectNumbers.Core/PartialFactorPendingEntry.cs new file mode 100644 index 0000000..25d1ce6 --- /dev/null +++ b/PerfectNumbers.Core/PartialFactorPendingEntry.cs @@ -0,0 +1,47 @@ +using System.Runtime.CompilerServices; + +namespace PerfectNumbers.Core; + +internal sealed class PartialFactorPendingEntry(ulong value, bool knownComposite) +{ + [ThreadStatic] + private static PartialFactorPendingEntry? _poolHead; + private PartialFactorPendingEntry? _next; + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static PartialFactorPendingEntry Rent(ulong value, bool knownComposite) + { + var poolHead = _poolHead; + if (poolHead == null) + { + return new(value, knownComposite); + } + + _poolHead = poolHead._next; + poolHead.Value = value; + poolHead.KnownComposite = knownComposite; + poolHead.HasKnownPrimality = knownComposite; + poolHead.IsPrime = false; + return poolHead; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public static void Return(PartialFactorPendingEntry entry) + { + entry._next = _poolHead; + _poolHead = entry; + } + + public ulong Value = value; + public bool KnownComposite = knownComposite; + public bool HasKnownPrimality = knownComposite; + public bool IsPrime = false; + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + public void WithPrimality(bool isPrime) + { + KnownComposite = KnownComposite || !isPrime; + HasKnownPrimality = true; + IsPrime = isPrime; + } +} diff --git a/PerfectNumbers.Core/PartialFactorResult.cs b/PerfectNumbers.Core/PartialFactorResult.cs index c19aefd..6ee5fb0 100644 --- a/PerfectNumbers.Core/PartialFactorResult.cs +++ b/PerfectNumbers.Core/PartialFactorResult.cs @@ -1,154 +1,109 @@ +using System.Collections.Generic; using System.Numerics; +using System.Runtime.CompilerServices; +using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core { public sealed class PartialFactorResult { + private const int ExponentHardLimit = 256; + [ThreadStatic] private static PartialFactorResult? s_poolHead; - public ulong[]? Factors; - public int[]? Exponents; - private PartialFactorResult? _next; - public ulong Cofactor; - public bool FullyFactored; - public int Count; - public bool CofactorIsPrime; - - private PartialFactorResult() - { - } - - public static PartialFactorResult Rent(in ulong[]? factors, in int[]? exponents, ulong cofactor, bool fullyFactored, int count, bool cofactorIsPrime) - { - PartialFactorResult? instance = s_poolHead; - if (instance is null) - { - instance = new PartialFactorResult(); - } - else - { - s_poolHead = instance._next; - } - - instance._next = null; - instance.Factors = factors; - instance.Exponents = exponents; - instance.Cofactor = cofactor; - instance.FullyFactored = fullyFactored; - instance.Count = count; - instance.CofactorIsPrime = cofactorIsPrime; - return instance; - } - - public static PartialFactorResult Rent(ulong cofactor, bool fullyFactored, int count, bool cofactorIsPrime) - { - PartialFactorResult? instance = s_poolHead; - if (instance is null) - { - instance = new PartialFactorResult(); - } - else - { - s_poolHead = instance._next; - } - - instance._next = null; - instance.Factors = null; - instance.Exponents = null; - instance.Cofactor = cofactor; - instance.FullyFactored = fullyFactored; - instance.Count = count; - instance.CofactorIsPrime = cofactorIsPrime; - return instance; - } - - public static readonly PartialFactorResult Empty = new() + public static readonly PartialFactorResult Empty = new(MontgomeryDivisorData.Empty) { Cofactor = 1UL, FullyFactored = true }; + internal readonly List Pending = new(PrimeOrderConstants.GpuSmallPrimeFactorSlots); + internal readonly FixedCapacityStack CompositeStack = new(PrimeOrderConstants.GpuSmallPrimeFactorSlots); + internal readonly ulong[] SpecialMaxBuffer = new ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; + + public readonly ulong[] Factors = new ulong[ExponentHardLimit]; + public readonly int[] Exponents = new int[ExponentHardLimit]; + public readonly ulong[] FactorCandidates = new ulong[ExponentHardLimit]; + public readonly List FactorCandidatesList = new(ExponentHardLimit); + public readonly ulong[] ValidateOrderFactorCandidates = new ulong[ExponentHardLimit]; + public readonly bool[] FactorCandidateEvaluations = new bool[ExponentHardLimit]; + internal readonly Dictionary FactorCounts = new(capacity: 8); + public ExponentRemainderStepperCpu ExponentRemainderStepper; + + private PartialFactorResult? _next; + public ulong Cofactor; + public bool CofactorIsPrime; + public int Count; + public bool FullyFactored; + public bool HasFactors; + + private PartialFactorResult(in MontgomeryDivisorData divisorData) + { + ExponentRemainderStepper = new(divisorData); + } + + public static PartialFactorResult Rent(in MontgomeryDivisorData divisorData) + { + if (s_poolHead is { } instance) + { + s_poolHead = instance._next; + instance._next = null; + instance.Cofactor = 0UL; + instance.FullyFactored = false; + instance.Count = 0; + instance.CofactorIsPrime = false; + instance.HasFactors = false; + if (!instance.ExponentRemainderStepper.MatchesDivisor(divisorData, 0UL)) + { + instance.ExponentRemainderStepper = new ExponentRemainderStepperCpu(divisorData, 0UL); + } + + return instance; + } + + return new PartialFactorResult(divisorData); + } + public void WithAdditionalPrime(ulong prime) { - ulong[]? sourceFactors = Factors; + ulong[] sourceFactors = Factors; int sourceCount = Count; - if (sourceFactors is null || sourceCount == 0) + if (sourceCount == 0) { InitializeWithPrime(prime); return; } - int newSize = (int)BitOperations.RoundUpToPowerOf2((uint)sourceCount + 1U); - ulong[] extendedFactors = sourceFactors.Length >= newSize - ? sourceFactors - : ResizeFactors(sourceFactors, sourceCount, newSize); + int[] sourceExponents = Exponents; - int[] sourceExponents = Exponents!; - int[] extendedExponents = sourceExponents.Length >= newSize - ? sourceExponents - : ResizeExponents(sourceCount, newSize, sourceExponents); + int newSize = sourceCount + 1; + sourceFactors[sourceCount] = prime; + sourceExponents[sourceCount] = 1; + Array.Sort(sourceFactors, sourceExponents, 0, newSize); - newSize = sourceCount + 1; - extendedFactors[sourceCount] = prime; - extendedExponents[sourceCount] = 1; - Array.Sort(extendedFactors, extendedExponents, 0, newSize); - - Factors = extendedFactors; - Exponents = extendedExponents; Cofactor = 1UL; FullyFactored = true; Count = newSize; CofactorIsPrime = false; - } - - private static int[] ResizeExponents(int sourceCount, int newSize, in int[] sourceExponents) - { - int[] extendedExponents = FixedCapacityPools.ExclusiveIntArray.Rent(newSize); - Array.Copy(sourceExponents, 0, extendedExponents, 0, sourceCount); - FixedCapacityPools.ExclusiveIntArray.Return(sourceExponents); - return extendedExponents; - } - - private static ulong[] ResizeFactors(in ulong[] sourceFactors, int sourceCount, int newSize) - { - ulong[] extendedFactors = FixedCapacityPools.ExclusiveUlongArray.Rent(newSize); - Array.Copy(sourceFactors, 0, extendedFactors, 0, sourceCount); - FixedCapacityPools.ExclusiveUlongArray.Return(sourceFactors); - return extendedFactors; + HasFactors = true; } private void InitializeWithPrime(ulong prime) { - ulong[] factors = FixedCapacityPools.ExclusiveUlongArray.Rent(1); - int[] exponents = FixedCapacityPools.ExclusiveIntArray.Rent(1); - factors[0] = prime; - exponents[0] = 1; - Factors = factors; - Exponents = exponents; + Factors[0] = prime; + Exponents[0] = 1; Cofactor = 1UL; FullyFactored = true; Count = 1; CofactorIsPrime = false; + HasFactors = true; } public void Dispose() { if (this != Empty) { - - ulong[]? factors = Factors; - - if (factors is not null) - { - int[] exponents = Exponents!; - Factors = null; - Exponents = null; - - FixedCapacityPools.ExclusiveUlongArray.Return(factors); - FixedCapacityPools.ExclusiveIntArray.Return(exponents); - } - _next = s_poolHead; s_poolHead = this; } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 95be940..ad9b3ea 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -1,59 +1,13 @@ -using System.Buffers; using System.Diagnostics; using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using PerfectNumbers.Core.Cpu; -using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - internal sealed class PendingEntry(ulong value, bool knownComposite) - { - [ThreadStatic] - private static PendingEntry? _poolHead; - private PendingEntry? _next; - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - public static PendingEntry Rent(ulong value, bool knownComposite) - { - var poolHead = _poolHead; - if (poolHead == null) - { - return new(value, knownComposite); - } - - _poolHead = poolHead._next; - poolHead.Value = value; - poolHead.KnownComposite = knownComposite; - poolHead.HasKnownPrimality = knownComposite; - poolHead.IsPrime = false; - return poolHead; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - public static void Return(PendingEntry entry) - { - entry._next = _poolHead; - _poolHead = entry; - } - - public ulong Value = value; - public bool KnownComposite = knownComposite; - public bool HasKnownPrimality = knownComposite; - public bool IsPrime = false; - - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - public void WithPrimality(bool isPrime) - { - KnownComposite = KnownComposite || !isPrime; - HasKnownPrimality = true; - IsPrime = isPrime; - } - } - public static ulong CalculateCpu( ulong prime, ulong? previousOrder, @@ -68,12 +22,12 @@ internal static partial class PrimeOrderCalculator ulong phi = prime - 1UL; - PartialFactorResult phiFactors = PartialFactorCpu(phi, config); + PartialFactorResult phiFactors = PartialFactorCpu(phi, divisorData, config); ulong result; - if (phiFactors.Factors is null) + if (!phiFactors.HasFactors) { - result = CalculateByFactorizationCpu(prime, divisorData); + result = CalculateByFactorizationCpu(prime, divisorData, phiFactors); phiFactors.Dispose(); return result; @@ -117,14 +71,14 @@ internal static partial class PrimeOrderCalculator else { divisorData = MontgomeryDivisorData.Empty; - result = CalculateWideInternalCpu(prime, previousOrder, divisorData, config); + result = CalculateWideInternalCpu(prime, previousOrder, config); } return result; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 CalculateWideInternalCpu(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) + private static UInt128 CalculateWideInternalCpu(in UInt128 prime, in UInt128? previousOrder, in PrimeOrderCalculatorConfig config) { if (prime <= UInt128.One) { @@ -138,17 +92,19 @@ internal static partial class PrimeOrderCalculator UInt128 phi = prime - UInt128.One; PartialFactorResult128 phiFactors = PartialFactorWide(phi, config); - UInt128 result; - if (phiFactors.Factors is null) - { - result = FinishStrictlyWideCpu(prime, divisorData); - } - else - { - result = RunHeuristicPipelineWideCpu(prime, previousOrder, divisorData, config, phi, phiFactors); - } - return result; + UInt128 result; + if (!phiFactors.HasFactors) + { + result = FinishStrictlyWideCpu(prime, phiFactors); + } + else + { + result = RunHeuristicPipelineWideCpu(prime, previousOrder, config, phi, phiFactors); + } + + phiFactors.Dispose(); + return result; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] @@ -160,13 +116,13 @@ internal static partial class PrimeOrderCalculator ulong phi, PartialFactorResult phiFactors) { - if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) + if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, phiFactors)) { return phi; } ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); - candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); + candidateOrder = ExponentLoweringCpu(candidateOrder, phiFactors, divisorData); if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) { @@ -176,7 +132,7 @@ internal static partial class PrimeOrderCalculator if (config.StrictMode) { orderFactors?.Dispose(); - return CalculateByFactorizationCpu(prime, divisorData); + return CalculateByFactorizationCpu(prime, divisorData, phiFactors); } if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, orderFactors, out ulong order)) @@ -191,61 +147,60 @@ internal static partial class PrimeOrderCalculator private static UInt128 RunHeuristicPipelineWideCpu( in UInt128 prime, in UInt128? previousOrder, - in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config, in UInt128 phi, in PartialFactorResult128 phiFactors) { - if (phiFactors.FullyFactored && TrySpecialMaxWideCpu(phi, prime, divisorData, phiFactors)) + if (phiFactors.FullyFactored && TrySpecialMaxWideCpu(phi, prime, phiFactors)) { return phi; } - UInt128 candidateOrder = InitializeStartingOrderWideCpu(prime, phi, divisorData); - candidateOrder = ExponentLoweringWideCpu(candidateOrder, prime, divisorData, phiFactors); + UInt128 candidateOrder = InitializeStartingOrderWideCpu(prime, phi); + candidateOrder = ExponentLoweringWideCpu(candidateOrder, prime, phiFactors); - if (TryConfirmOrderWideCpu(prime, candidateOrder, divisorData, config)) + if (TryConfirmOrderWideCpu(prime, candidateOrder, config, out PartialFactorResult128? orderFactors)) { return candidateOrder; } if (config.StrictMode) { - return FinishStrictlyWideCpu(prime, divisorData); + orderFactors?.Dispose(); + return FinishStrictlyWideCpu(prime, phiFactors); } - if (TryHeuristicFinishWideCpu(prime, candidateOrder, previousOrder, divisorData, config, out UInt128 order)) + if (TryHeuristicFinishWideCpu(prime, candidateOrder, previousOrder, config, orderFactors, out UInt128 order)) { + orderFactors?.Dispose(); return order; } + orderFactors?.Dispose(); return candidateOrder; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 FinishStrictlyWideCpu(in UInt128 prime, in MontgomeryDivisorData? divisorData) + private static UInt128 FinishStrictlyWideCpu(in UInt128 prime, PartialFactorResult128 phiFactors) { UInt128 phi = prime - UInt128.One; - FixedCapacityStack> uInt128IntDictionaryPool = ThreadStaticPools.UInt128IntDictionaryPool; - Dictionary counts = uInt128IntDictionaryPool.Rent(capacity: 8); + Dictionary counts = phiFactors.FactorCounts; + counts.Clear(); FactorCompletelyWide(phi, counts); int entryCount = counts.Count; if (entryCount == 0) { - uInt128IntDictionaryPool.Return(counts); return phi; } - FixedCapacityStack>> keyValuePairUInt128IntegerPool = ThreadStaticPools.KeyValuePairUInt128IntegerPool; - var entries = keyValuePairUInt128IntegerPool.Rent(entryCount); - entries.AddRange(counts); - uInt128IntDictionaryPool.Return(counts); + Span> entries = [.. counts]; entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); UInt128 order = phi; for (int i = 0; i < entryCount; i++) { + // (UInt128 primeFactor, int exponent) = entries[i]; UInt128 primeFactor = entries[i].Key; int exponent = entries[i].Value; for (int iteration = 0; iteration < exponent; iteration++) @@ -253,7 +208,7 @@ internal static partial class PrimeOrderCalculator if (order.ReduceCycleRemainder(primeFactor) != UInt128.Zero) { UInt128 candidate = order / primeFactor; - if (Pow2ModWideCpu(candidate, prime, divisorData) == UInt128.One) + if (Pow2ModWideCpu(candidate, prime) == UInt128.One) { order = candidate; } @@ -265,12 +220,11 @@ internal static partial class PrimeOrderCalculator } } - keyValuePairUInt128IntegerPool.Return(entries); return order; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TrySpecialMaxWideCpu(in UInt128 phi, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) + private static bool TrySpecialMaxWideCpu(in UInt128 phi, in UInt128 prime, in PartialFactorResult128 factors) { ReadOnlySpan factorSpan = factors.Factors; int length = factors.Count; @@ -278,7 +232,7 @@ internal static partial class PrimeOrderCalculator { UInt128 factor = factorSpan[i].Value; UInt128 reduced = phi / factor; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) + if (Pow2ModWideCpu(reduced, prime) == UInt128.One) { return false; } @@ -288,14 +242,14 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 InitializeStartingOrderWideCpu(in UInt128 prime, in UInt128 phi, in MontgomeryDivisorData? divisorData) + private static UInt128 InitializeStartingOrderWideCpu(in UInt128 prime, in UInt128 phi) { UInt128 order = phi; UInt128 mod8 = prime & (UInt128)7UL; if (mod8 == UInt128.One || mod8 == (UInt128)7UL) { UInt128 half = phi >> 1; - if (Pow2ModWideCpu(half, prime, divisorData) == UInt128.One) + if (Pow2ModWideCpu(half, prime) == UInt128.One) { order = half; } @@ -305,9 +259,8 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 ExponentLoweringWideCpu(UInt128 order, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) + private static UInt128 ExponentLoweringWideCpu(UInt128 order, in UInt128 prime, in PartialFactorResult128 factors) { - ArrayPool pool = ThreadStaticPools.FactorEntry128Pool; ReadOnlySpan factorSpan = factors.Factors; int length = factors.Count; // TODO: This should never trigger from production code - check @@ -315,14 +268,15 @@ internal static partial class PrimeOrderCalculator // { // return order; // } - FactorEntry128[] tempArray = pool.Rent(length + 1); - Span buffer = tempArray.AsSpan(0, length); + + int newLength = length + 1; + Span buffer = new(factors.TempFactors, 0, newLength); factorSpan.CopyTo(buffer); if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor)) { buffer[length] = new FactorEntry128(factors.Cofactor, 1); - length++; + length = newLength; } buffer[..length].Sort(static (a, b) => a.Value.CompareTo(b.Value)); @@ -336,7 +290,7 @@ internal static partial class PrimeOrderCalculator if (order.ReduceCycleRemainder(primeFactor) == UInt128.Zero) { UInt128 reduced = order / primeFactor; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) + if (Pow2ModWideCpu(reduced, prime) == UInt128.One) { order = reduced; continue; @@ -347,26 +301,27 @@ internal static partial class PrimeOrderCalculator } } - pool.Return(tempArray, clearArray: false); return order; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmOrderWideCpu(in UInt128 prime, in UInt128 order, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) + private static bool TryConfirmOrderWideCpu(in UInt128 prime, in UInt128 order, in PrimeOrderCalculatorConfig config, out PartialFactorResult128? reusableFactorization) { + reusableFactorization = null; if (order == UInt128.Zero) { return false; } - if (Pow2ModWideCpu(order, prime, divisorData) != UInt128.One) + if (Pow2ModWideCpu(order, prime) != UInt128.One) { return false; } PartialFactorResult128 factorization = PartialFactorWide(order, config); - if (factorization.Factors is null) + if (factorization.Count == 0) { + factorization.Dispose(); return false; } @@ -374,19 +329,21 @@ internal static partial class PrimeOrderCalculator { if (factorization.Cofactor <= UInt128.One) { + factorization.Dispose(); return false; } if (!IsPrimeWide(factorization.Cofactor)) { + factorization.Dispose(); return false; } - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.WithAdditionalPrime(factorization.Cofactor); } - ReadOnlySpan span = factorization.Factors; int length = factorization.Count; + ReadOnlySpan span = new(factorization.Factors, 0, length); for (int i = 0; i < length; i++) { UInt128 primeFactor = span[i].Value; @@ -399,13 +356,15 @@ internal static partial class PrimeOrderCalculator } reduced /= primeFactor; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) + if (Pow2ModWideCpu(reduced, prime) == UInt128.One) { + reusableFactorization = factorization; return false; } } } + factorization.Dispose(); return true; } @@ -414,8 +373,8 @@ internal static partial class PrimeOrderCalculator in UInt128 prime, in UInt128 order, in UInt128? previousOrder, - in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config, + PartialFactorResult128? cachedOrderFactors, out UInt128 result) { result = UInt128.Zero; @@ -424,80 +383,100 @@ internal static partial class PrimeOrderCalculator return false; } - PartialFactorResult128 orderFactors = PartialFactorWide(order, config); - if (orderFactors.Factors is null) + PartialFactorResult128 orderFactors; + if (cachedOrderFactors is not null) { + orderFactors = cachedOrderFactors; + } + else + { + orderFactors = PartialFactorWide(order, config); + } + + if (orderFactors.Count == 0) + { + if (cachedOrderFactors is null) + { + orderFactors.Dispose(); + } + return false; } - if (!orderFactors.FullyFactored) + try { - if (orderFactors.Cofactor <= UInt128.One) + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= UInt128.One) + { + return false; + } + + if (!IsPrimeWide(orderFactors.Cofactor)) + { + return false; + } + + orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; + List candidates = orderFactors.FactorCandidatesList; + candidates.Clear(); + FactorEntry128[] factorArray = orderFactors.Factors!; + + BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity); + + int candidateCount = candidates.Count; + if (candidateCount == 0) { return false; } - if (!IsPrimeWide(orderFactors.Cofactor)) + SortCandidatesWide(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; + int powUsed = 0; + + for (int i = 0; i < candidateCount; i++) { - return false; + if (powUsed >= powBudget && powBudget > 0) + { + break; + } + + UInt128 candidate = candidates[i]; + powUsed++; + + if (Pow2ModWideCpu(candidate, prime) != UInt128.One) + { + continue; + } + + if (!TryConfirmCandidateWideCpu(prime, candidate, config, ref powUsed, powBudget)) + { + continue; + } + + result = candidate; + return true; } - orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - FixedCapacityStack> uInt128ListPool = ThreadStaticPools.UInt128ListPool; - List candidates = uInt128ListPool.Rent(capacity); - FactorEntry128[] factorArray = orderFactors.Factors!; - - BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity); - - int candidateCount = candidates.Count; - if (candidateCount == 0) - { - uInt128ListPool.Return(candidates); return false; } - - SortCandidatesWide(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; - int powUsed = 0; - - for (int i = 0; i < candidateCount; i++) + finally { - if (powUsed >= powBudget && powBudget > 0) + if (cachedOrderFactors is null) { - break; + orderFactors.Dispose(); } - - UInt128 candidate = candidates[i]; - powUsed++; - - if (Pow2ModWideCpu(candidate, prime, divisorData) != UInt128.One) - { - continue; - } - - if (!TryConfirmCandidateWideCpu(prime, candidate, config, ref powUsed, powBudget, divisorData)) - { - continue; - } - - result = candidate; - uInt128ListPool.Return(candidates); - return true; } - - uInt128ListPool.Return(candidates); - return false; } // private static ulong specialMaxHits; [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + private static bool TrySpecialMaxCpu(ulong phi, PartialFactorResult factors) { - 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) @@ -505,16 +484,16 @@ internal static partial class PrimeOrderCalculator // return true; // } - ReadOnlySpan factorSpan = new(factors.Factors, 0, length); - Span stackBuffer = stackalloc ulong[length]; - - return EvaluateSpecialMaxCandidatesCpu(stackBuffer, factorSpan, phi, prime, divisorData); + return EvaluateSpecialMaxCandidatesCpu(factors, phi); } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool EvaluateSpecialMaxCandidatesCpu(Span buffer, ReadOnlySpan factors, ulong phi, ulong prime, in MontgomeryDivisorData divisorData) + private static bool EvaluateSpecialMaxCandidatesCpu(PartialFactorResult partialFactors, ulong phi) { - int factorCount = factors.Length; + int factorCount = partialFactors.Count; + ReadOnlySpan factors = new(partialFactors.Factors, 0, factorCount); + // Span buffer = new(partialFactors.SpecialMaxBuffer, 0, factorCount); + Span buffer = stackalloc ulong[factorCount]; for (int i = 0; i < factorCount; i++) { ulong factor = factors[i]; @@ -535,11 +514,11 @@ internal static partial class PrimeOrderCalculator buffer.Sort(); - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + ExponentRemainderStepperCpu stepper = partialFactors.ExponentRemainderStepper; + stepper.Reset(); if (stepper.InitializeCpuIsUnity(buffer[0])) { - ThreadStaticPools.ReturnExponentStepperCpu(stepper); return false; } @@ -547,12 +526,10 @@ internal static partial class PrimeOrderCalculator { if (stepper.ComputeNextIsUnity(buffer[i])) { - ThreadStaticPools.ReturnExponentStepperCpu(stepper); return false; } } - ThreadStaticPools.ReturnExponentStepperCpu(stepper); return true; } @@ -573,76 +550,45 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + private static ulong ExponentLoweringCpu(ulong order, PartialFactorResult partialFactors, in MontgomeryDivisorData divisorData) { - ReadOnlySpan factorSpan = factors.Factors; - ReadOnlySpan exponentsSpan = factors.Exponents; + int length = partialFactors.Count; + ulong[] factors = partialFactors.Factors; + int[] exponents = partialFactors.Exponents; + Span factorsSpan = new(factors, 0, length + 1); + Span exponentsSpan = new(exponents, 0, length + 1); - int length = factors.Count; int newLength = length + 1; - ulong[] tempFactorsArray = FixedCapacityPools.ExclusiveUlongArray.Rent(newLength); - int[] tempExponentsArray = FixedCapacityPools.ExclusiveIntArray.Rent(newLength); - Span ulongBuffer = tempFactorsArray; - factorSpan.CopyTo(ulongBuffer); - Span intBuffer = tempExponentsArray; - exponentsSpan.CopyTo(intBuffer); - - if (!factors.FullyFactored && factors.Cofactor > 1UL && factors.CofactorIsPrime) + if (!partialFactors.FullyFactored && partialFactors.Cofactor > 1UL && partialFactors.CofactorIsPrime) { - ulongBuffer[length] = factors.Cofactor; - intBuffer[length] = 1; + factorsSpan[length] = partialFactors.Cofactor; + exponentsSpan[length] = 1; length = newLength; } - Array.Sort(tempFactorsArray, tempExponentsArray, 0, length); + Array.Sort(factors, exponents, 0, length); - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + ExponentRemainderStepperCpu stepper = partialFactors.ExponentRemainderStepper; - const int StackExponentCapacity = 16; - const int ExponentHardLimit = 256; - - Span stackCandidates = stackalloc ulong[StackExponentCapacity]; - Span stackEvaluations = stackalloc bool[StackExponentCapacity]; - - ulong[]? heapCandidateArray = null; - bool[]? heapEvaluationArray = null; + // Span stackCandidates = new(partialFactors.FactorCandidates); + // Span stackEvaluations = new(partialFactors.FactorCandidateEvaluations); + Span stackCandidates = stackalloc ulong[32]; + Span stackEvaluations = stackalloc bool[32]; for (int i = 0; i < length; i++) { - ulong primeFactor = ulongBuffer[i]; - int exponent = intBuffer[i]; + ulong primeFactor = factorsSpan[i]; + int exponent = exponentsSpan[i]; // Factor exponents produced by Pollard-Rho and trial division are always positive in the scanner flow. // if (exponent <= 0) // { // continue; // } - if (exponent > ExponentHardLimit) - { - throw new InvalidOperationException($"Prime factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}."); - } - - if (exponent <= StackExponentCapacity) - { - ProcessExponentLoweringPrime(stackCandidates[..exponent], stackEvaluations[..exponent], ref order, primeFactor, exponent, ref stepper); - continue; - } - - heapCandidateArray ??= FixedCapacityPools.ExclusiveUlongArray.Rent(ExponentHardLimit); - heapEvaluationArray ??= ThreadStaticPools.BoolPool.Rent(ExponentHardLimit); - ProcessExponentLoweringPrime(heapCandidateArray.AsSpan(0, exponent), heapEvaluationArray.AsSpan(0, exponent), ref order, primeFactor, exponent, ref stepper); + ProcessExponentLoweringPrime(stackCandidates[..exponent], stackEvaluations[..exponent], ref order, primeFactor, exponent, ref stepper); } - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - - if (heapCandidateArray is not null) - { - FixedCapacityPools.ExclusiveUlongArray.Return(heapCandidateArray); - ThreadStaticPools.BoolPool.Return(heapEvaluationArray!); - } - - FixedCapacityPools.ExclusiveUlongArray.Return(tempFactorsArray); return order; } @@ -674,26 +620,23 @@ internal static partial class PrimeOrderCalculator return; } - Span candidates = candidateBuffer[..actual]; - Span evaluations = evaluationBuffer[..actual]; - stepper.Reset(); int last = actual - 1; - evaluations[last] = stepper.InitializeCpuIsUnity(candidates[last]); + evaluationBuffer[last] = stepper.InitializeCpuIsUnity(candidateBuffer[last]); for (int j = last - 1; j >= 0; j--) { - evaluations[j] = stepper.ComputeNextIsUnity(candidates[j]); + evaluationBuffer[j] = stepper.ComputeNextIsUnity(candidateBuffer[j]); } for (int j = 0; j < actual; j++) { - if (!evaluations[j]) + if (!evaluationBuffer[j]) { break; } - order = candidates[j]; + order = candidateBuffer[j]; } } @@ -713,8 +656,8 @@ internal static partial class PrimeOrderCalculator return false; } - PartialFactorResult factorization = PartialFactorCpu(order, config); - if (factorization.Factors is null) + PartialFactorResult factorization = PartialFactorCpu(order, divisorData, config); + if (factorization.Count == 0) { factorization.Dispose(); return false; @@ -731,7 +674,7 @@ internal static partial class PrimeOrderCalculator factorization.WithAdditionalPrime(factorization.Cofactor); } - if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization)) + if (!ValidateOrderAgainstFactors(order, factorization)) { reusableFactorization = factorization; return false; @@ -743,25 +686,17 @@ internal static partial class PrimeOrderCalculator [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static bool ValidateOrderAgainstFactors( - ulong prime, ulong order, - in MontgomeryDivisorData divisorData, PartialFactorResult factorization) { - ReadOnlySpan factorsSpan = factorization.Factors!; - ReadOnlySpan exponentsSpan = factorization.Exponents!; - int length = factorization.Count; + ReadOnlySpan factorsSpan = new(factorization.Factors, 0, length); + ReadOnlySpan exponentsSpan = new (factorization.Exponents, 0, length); + // Span stackBuffer = new(factorization.ValidateOrderFactorCandidates); + Span stackBuffer = stackalloc ulong[32]; + ExponentRemainderStepperCpu stepper = factorization.ExponentRemainderStepper; - 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]; @@ -773,37 +708,13 @@ internal static partial class PrimeOrderCalculator // continue; // } - if (exponent > ExponentHardLimit) - { - throw new InvalidOperationException($"Factor exponent {exponent} exceeds the supported limit of {ExponentHardLimit}."); - } - - if (exponent <= StackExponentCapacity) - { - if (ValidateOrderForFactor(stackBuffer[..exponent], primeFactor, exponent, order, ref stepper)) - { - violates = true; - break; - } - - continue; - } - - heapCandidateArray ??= FixedCapacityPools.ExclusiveUlongArray.Rent(ExponentHardLimit); - if (ValidateOrderForFactor(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, order, ref stepper)) + if (ValidateOrderForFactor(stackBuffer[..exponent], primeFactor, exponent, order, ref stepper)) { violates = true; break; } } - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - - if (heapCandidateArray is not null) - { - FixedCapacityPools.ExclusiveUlongArray.Return(heapCandidateArray); - } - return !violates; } @@ -874,10 +785,10 @@ internal static partial class PrimeOrderCalculator } // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorCpu(order, config); + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorCpu(order, divisorData, config); try { - if (orderFactors.Factors is null) + if (orderFactors.Count == 0) { return false; } @@ -898,16 +809,13 @@ internal static partial class PrimeOrderCalculator } int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); + List candidates = orderFactors.FactorCandidatesList; candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); + BuildCandidates(order, orderFactors, orderFactors.Count, candidates, capacity); int candidateCount = candidates.Count; if (candidateCount == 0) { - ThreadStaticPools.ReturnUlongList(candidates); return false; } @@ -916,7 +824,8 @@ internal static partial class PrimeOrderCalculator int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; int powUsed = 0; Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + ExponentRemainderStepperCpu powStepper = orderFactors.ExponentRemainderStepper; + powStepper.Reset(); ulong candidate = candidateSpan[0]; powUsed++; @@ -924,9 +833,6 @@ internal static partial class PrimeOrderCalculator bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); if (equalsOne && TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) { - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); result = candidate; return true; } @@ -954,18 +860,11 @@ internal static partial class PrimeOrderCalculator continue; } - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); result = candidate; return true; } - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - return false; } finally @@ -1096,13 +995,15 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static void BuildCandidates(ulong order, ulong[] factors, int[] exponents, int count, List candidates, int limit) + private static void BuildCandidates(ulong order, PartialFactorResult orderFactors, int count, List candidates, int limit) { if (count == 0) { return; } + ulong[] factors = orderFactors.Factors; + int[] exponents = orderFactors.Exponents; Array.Sort(factors, exponents, 0, count); BuildCandidatesRecursive(order, factors, exponents, 0, 1UL, candidates, limit); } @@ -1171,9 +1072,9 @@ internal static partial class PrimeOrderCalculator [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) { - PartialFactorResult factorization = PartialFactorCpu(candidate, config); + PartialFactorResult factorization = PartialFactorCpu(candidate, divisorData, config); - if (factorization.Factors is null) + if (factorization.Count == 0) { factorization.Dispose(); return false; @@ -1190,16 +1091,13 @@ internal static partial class PrimeOrderCalculator factorization.WithAdditionalPrime(factorization.Cofactor); } - ReadOnlySpan factorsSpan = factorization.Factors; - ReadOnlySpan exponentsSpan = factorization.Exponents; int length = factorization.Count; + ReadOnlySpan factorsSpan = new(factorization.Factors, 0, length); + ReadOnlySpan exponentsSpan = new(factorization.Exponents, 0, length); - const int StackExponentCapacity = 32; - const int ExponentHardLimit = 256; - - Span stackBuffer = stackalloc ulong[StackExponentCapacity]; - ulong[]? heapCandidateArray = null; - ExponentRemainderStepperCpu stepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + // Span stackBuffer = new(factorization.SpecialMaxBuffer); + Span stackBuffer = stackalloc ulong[32]; + ExponentRemainderStepperCpu stepper = factorization.ExponentRemainderStepper; bool violates = false; for (int i = 0; i < length; i++) @@ -1213,37 +1111,13 @@ internal static partial class PrimeOrderCalculator // 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 ??= FixedCapacityPools.ExclusiveUlongArray.Rent(ExponentHardLimit); - if (CheckCandidateViolation(heapCandidateArray.AsSpan(0, exponent), primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) + if (CheckCandidateViolation(stackBuffer[..exponent], primeFactor, exponent, candidate, prime, ref powUsed, powBudget, ref stepper)) { violates = true; break; } } - ThreadStaticPools.ReturnExponentStepperCpu(stepper); - - if (heapCandidateArray is not null) - { - FixedCapacityPools.ExclusiveUlongArray.Return(heapCandidateArray); - } - factorization.Dispose(); return !violates; } @@ -1335,22 +1209,23 @@ internal static partial class PrimeOrderCalculator => exponent.Pow2MontgomeryModWindowedConvertToStandardCpu(divisorData) == 1UL; [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorCpu(ulong value, in PrimeOrderCalculatorConfig config) + private static PartialFactorResult PartialFactorCpu(ulong value, in MontgomeryDivisorData divisorData, 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]; + var result = PartialFactorResult.Rent(divisorData); + Span primeSlots = new(result.Factors); + Span exponentSlots = new(result.Exponents); // We don't need to worry about leftovers, because we always use indexes within the calculated counts // primeSlots.Clear(); // exponentSlots.Clear(); - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(PrimeOrderConstants.GpuSmallPrimeFactorSlots); + List pending = result.Pending; + pending.Clear(); PopulateSmallPrimeFactorsCpu( value, @@ -1366,7 +1241,8 @@ internal static partial class PrimeOrderCalculator if (config.PollardRhoMilliseconds > 0) { long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - FixedCapacityStack compositeStack = ThreadStaticPools.RentUlongStack(PrimeOrderConstants.GpuSmallPrimeFactorSlots); + FixedCapacityStack compositeStack = result.CompositeStack; + compositeStack.Clear(); compositeStack.Push(remaining); CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); @@ -1375,38 +1251,29 @@ internal static partial class PrimeOrderCalculator { while (compositeStack.Count > 0) { - pending.Add(PendingEntry.Rent(compositeStack.Pop(), knownComposite: false)); + pending.Add(PartialFactorPendingEntry.Rent(compositeStack.Pop(), knownComposite: false)); } } - - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); } else { - pending.Add(PendingEntry.Rent(remaining, knownComposite: false)); + pending.Add(PartialFactorPendingEntry.Rent(remaining, knownComposite: false)); } } ulong cofactor = 1UL; bool cofactorContainsComposite = false; int pendingCount = pending.Count; - int index = 0; - if (pending.Capacity > PrimeOrderConstants.GpuSmallPrimeFactorSlots) + for (int index = 0; index < pendingCount; index++) { - Console.WriteLine($"Initial capacity is too small. Required capacity: {pending.Capacity}"); - } - - for (; index < pendingCount; index++) - { - PendingEntry entry = pending[index]; + PartialFactorPendingEntry entry = pending[index]; ulong composite = entry.Value; if (entry.KnownComposite) { cofactor = checked(cofactor * composite); cofactorContainsComposite = true; - PendingEntry.Return(entry); + PartialFactorPendingEntry.Return(entry); continue; } @@ -1419,7 +1286,7 @@ internal static partial class PrimeOrderCalculator if (entry.IsPrime) { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); } // composite is never smaller on the execution path // else if (composite > 1UL) @@ -1429,7 +1296,7 @@ internal static partial class PrimeOrderCalculator cofactorContainsComposite = true; } - PendingEntry.Return(entry); + PartialFactorPendingEntry.Return(entry); } bool cofactorIsPrime; @@ -1446,19 +1313,6 @@ internal static partial class PrimeOrderCalculator 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) // { @@ -1472,37 +1326,18 @@ internal static partial class PrimeOrderCalculator // goto ReturnResult; // } - ulong[] factorsArray = FixedCapacityPools.ExclusiveUlongArray.Rent(factorCount); - int[] exponentsArray = FixedCapacityPools.ExclusiveIntArray.Rent(factorCount); + Array.Sort(result.Factors, result.Exponents, 0, 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 it factored completely - limitReached = cofactor == 1UL; - - result = PartialFactorResult.Rent(factorsArray, exponentsArray, cofactor, limitReached, factorCount, cofactorIsPrime); - - ReturnResult: - pending.Clear(); - ThreadStaticPools.ReturnPrimeOrderPendingEntryList(pending); + result.Cofactor = cofactor; + result.FullyFactored = cofactor == 1UL; + result.Count = factorCount; + result.HasFactors = factorCount != 0; + result.CofactorIsPrime = cofactorIsPrime; return result; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void CollectFactorsCpu(Span primeSlots, Span exponentSlots, ref int factorCount, List pending, FixedCapacityStack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) + private static void CollectFactorsCpu(Span primeSlots, Span exponentSlots, ref int factorCount, List pending, FixedCapacityStack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) { while (compositeStack.Count > 0) { @@ -1511,20 +1346,20 @@ internal static partial class PrimeOrderCalculator if (isPrime) { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); continue; } pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; if (pollardRhoDeadlineReached) { - pending.Add(PendingEntry.Rent(composite, knownComposite: true)); + pending.Add(PartialFactorPendingEntry.Rent(composite, knownComposite: true)); continue; } if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) { - pending.Add(PendingEntry.Rent(composite, knownComposite: true)); + pending.Add(PartialFactorPendingEntry.Rent(composite, knownComposite: true)); continue; } @@ -1590,10 +1425,6 @@ internal static partial class PrimeOrderCalculator continue; } - EnsureCollectorCapacity( - primeTargets, - factorCount + 1); - primeTargets[factorCount] = primeValue; exponentTargets[factorCount] = exponent; factorCount++; @@ -1678,6 +1509,7 @@ internal static partial class PrimeOrderCalculator long startTimestamp = Stopwatch.GetTimestamp(); long deadline = ConvertMillisecondsToStopwatchTicks(milliseconds); deadline += startTimestamp; + // Return maximum value in case of overflow. if (deadline < startTimestamp) { return long.MaxValue; @@ -1801,62 +1633,65 @@ internal static partial class PrimeOrderCalculator Span primes, Span exponents, ref int count, - ulong prime, - int exponent) + ulong prime) { - if (exponent <= 0) - { - return; - } + // exponent is always 1 on the execution path. + // if (exponent <= 0) + // { + // return; + // } for (int i = 0; i < count; i++) { if (primes[i] == prime) { - exponents[i] += exponent; + exponents[i]++; return; } } - EnsureCollectorCapacity( - primes, - count + 1); - primes[count] = prime; - exponents[count] = exponent; + exponents[count] = 1; count++; } - [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - private static void EnsureCollectorCapacity( - Span primes, - int requiredCapacity) - { - ArgumentOutOfRangeException.ThrowIfGreaterThan(requiredCapacity, primes.Length); - } - [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData) + private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData, PartialFactorResult phiFactors) { ulong phi = prime - 1UL; - Dictionary counts = ThreadStaticPools.RentFactorCountDictionary(); + Dictionary counts = phiFactors.FactorCounts; + counts.Clear(); FactorCompletelyCpu(phi, counts, false); if (counts.Count == 0) { - ThreadStaticPools.ReturnFactorCountDictionary(counts); return phi; } - KeyValuePair[] entries = [.. counts]; - Array.Sort(entries, static (a, b) => a.Key.CompareTo(b.Key)); + // int entryCount = counts.Count; + // Span primes = new(phiFactors.Factors, 0, entryCount); + // Span exponents = new(phiFactors.Exponents, 0, entryCount); + // int index = 0; + // foreach (KeyValuePair kvp in counts) + // { + // primes[index] = kvp.Key; + // exponents[index] = kvp.Value; + // index++; + // } + // Array.Sort(phiFactors.Factors, phiFactors.Exponents, 0, entryCount); + + Span> entries = [.. counts]; + + entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); ulong order = phi; int entryCount = entries.Length; for (int i = 0; i < entryCount; i++) { + // (ulong primeFactor, int exponent) = entries[i]; ulong primeFactor = entries[i].Key; int exponent = entries[i].Value; + for (int iteration = 0; iteration < exponent; iteration++) { ulong remainder = order.ReduceCycleRemainder(primeFactor); @@ -1879,7 +1714,6 @@ internal static partial class PrimeOrderCalculator } } - ThreadStaticPools.ReturnFactorCountDictionary(counts); return order; } @@ -1974,7 +1808,7 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static UInt128 Pow2ModWideCpu(in UInt128 exponent, in UInt128 modulus, in MontgomeryDivisorData? divisorData) + private static UInt128 Pow2ModWideCpu(in UInt128 exponent, in UInt128 modulus) { if (modulus == UInt128.One) { @@ -1985,60 +1819,64 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TryConfirmCandidateWideCpu(in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget, in MontgomeryDivisorData? divisorData) + private static bool TryConfirmCandidateWideCpu(in UInt128 prime, in UInt128 candidate, in PrimeOrderCalculatorConfig config, ref int powUsed, int powBudget) { PartialFactorResult128 factorization = PartialFactorWide(candidate, config); - if (factorization.Factors is null) + if (factorization.Count == 0) { + factorization.Dispose(); return false; } - if (!factorization.FullyFactored) + try { - if (factorization.Cofactor <= UInt128.One) + if (!factorization.FullyFactored) { - 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.ReduceCycleRemainder(primeFactor) != UInt128.Zero) - { - break; - } - - reduced /= primeFactor; - if (powUsed >= powBudget && powBudget > 0) + if (factorization.Cofactor <= UInt128.One) { return false; } - powUsed++; - if (Pow2ModWideCpu(reduced, prime, divisorData) == UInt128.One) + if (!IsPrimeWide(factorization.Cofactor)) { return false; } - } - } - return true; + factorization.WithAdditionalPrime(factorization.Cofactor); + } + + int length = factorization.Count; + ReadOnlySpan span = new(factorization.Factors, 0, length); + 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.ReduceCycleRemainder(primeFactor) != UInt128.Zero) + { + break; + } + + reduced /= primeFactor; + if (powUsed >= powBudget && powBudget > 0) + { + return false; + } + + powUsed++; + if (Pow2ModWideCpu(reduced, prime) == UInt128.One) + { + return false; + } + } + } + + return true; + } + finally + { + factorization.Dispose(); + } } - } - - - diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs index c2349b6..7f86123 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -76,12 +76,12 @@ internal static partial class PrimeOrderCalculator return gpuOrder; } - PartialFactorResult phiFactors = PartialFactorGpu(gpu, phi, config); + PartialFactorResult phiFactors = PartialFactorGpu(gpu, divisorData, phi, config); ulong result; - if (phiFactors.Factors is null) + if (phiFactors.Count == 0) { - result = CalculateByFactorizationCpu(prime, divisorData); + result = CalculateByFactorizationCpu(prime, divisorData, phiFactors); phiFactors.Dispose(); return result; @@ -147,21 +147,28 @@ internal static partial class PrimeOrderCalculator UInt128 phi = prime - UInt128.One; PartialFactorResult128 phiFactors = PartialFactorWide(phi, config); - UInt128 result; - if (phiFactors.Factors is null) + try { - result = FinishStrictlyWideGpu(gpu, prime, divisorData); - } - else - { - result = RunHeuristicPipelineWideGpu(gpu, prime, previousOrder, divisorData, config, phi, phiFactors); - } + UInt128 result; + if (phiFactors.Count == 0) + { + result = FinishStrictlyWideGpu(gpu, prime, divisorData); + } + else + { + result = RunHeuristicPipelineWideGpu(gpu, prime, previousOrder, divisorData, config, phi, phiFactors); + } - return result; + return result; + } + finally + { + phiFactors.Dispose(); + } } [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) + private static void CollectFactorsGpu(PrimeOrderCalculatorAccelerator gpu, ref Span primeSlots, ref Span exponentSlots, ref int factorCount, List pending, Stack compositeStack, long deadlineTimestamp, ref bool pollardRhoDeadlineReached) { while (compositeStack.Count > 0) { @@ -170,7 +177,7 @@ internal static partial class PrimeOrderCalculator if (isPrime) { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); continue; } @@ -178,14 +185,14 @@ internal static partial class PrimeOrderCalculator if (currentTimestamp > deadlineTimestamp) { pollardRhoDeadlineReached = true; - pending.Add(new PendingEntry(composite, knownComposite: true)); + pending.Add(new PartialFactorPendingEntry(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)); + pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true)); continue; } @@ -200,7 +207,6 @@ internal static partial class PrimeOrderCalculator PrimeOrderCalculatorAccelerator gpu, ReadOnlySpan factors, ulong phi, - ulong prime, in MontgomeryDivisorData divisorData) { int factorCount = factors.Length; @@ -235,7 +241,6 @@ internal static partial class PrimeOrderCalculator private static UInt128 ExponentLoweringWideGpu(PrimeOrderCalculatorAccelerator gpu, UInt128 order, in UInt128 prime, in MontgomeryDivisorData? divisorData, in PartialFactorResult128 factors) { - ArrayPool pool = ThreadStaticPools.FactorEntry128Pool; ReadOnlySpan factorSpan = factors.Factors; int length = factors.Count; // TODO: This should never trigger from production code - check @@ -243,9 +248,7 @@ internal static partial class PrimeOrderCalculator // { // return order; // } - FactorEntry128[] tempArray = pool.Rent(length + 1); - - Span buffer = tempArray.AsSpan(0, length); + Span buffer = factors.TempFactors.AsSpan(0, length + 1); factorSpan.CopyTo(buffer); if (!factors.FullyFactored && factors.Cofactor > UInt128.One && IsPrimeWide(factors.Cofactor)) @@ -263,7 +266,7 @@ internal static partial class PrimeOrderCalculator int exponent = buffer[i].Exponent; for (int iteration = 0; iteration < exponent; iteration++) { - if ((order % primeFactor) == UInt128.Zero) + if (order.ReduceCycleRemainder(primeFactor) == UInt128.Zero) { UInt128 reduced = order / primeFactor; if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) @@ -277,7 +280,6 @@ internal static partial class PrimeOrderCalculator } } - pool.Return(tempArray, clearArray: false); return order; } @@ -302,7 +304,7 @@ internal static partial class PrimeOrderCalculator int exponent = entries[i].Value; for (int iteration = 0; iteration < exponent; iteration++) { - if ((order % primeFactor) != UInt128.Zero) + if (order.ReduceCycleRemainder(primeFactor) != UInt128.Zero) { UInt128 candidate = order / primeFactor; if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) == UInt128.One) @@ -337,28 +339,26 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) + private static PartialFactorResult PartialFactorGpu(PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisorData, 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 = []; + PartialFactorResult result = PartialFactorResult.Rent(divisorData); + Span primeSlots = new(result.Factors); + Span exponentSlots = new (result.Exponents); // 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; + uint limit = config.SmallFactorLimit; + bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out int factorCount, out ulong remaining); - bool gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(gpu, value, limit, primeSlots, exponentSlots, out int factorCount, out remaining); - - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); + List pending = result.Pending; + pending.Clear(); if (!gpuFactored) { @@ -371,31 +371,28 @@ internal static partial class PrimeOrderCalculator out remaining); } - bool limitReached; if (remaining > 1UL) { if (config.PollardRhoMilliseconds > 0) { long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - FixedCapacityStack compositeStack = ThreadStaticPools.RentUlongStack(4); + FixedCapacityStack compositeStack = result.CompositeStack; + compositeStack.Clear(); compositeStack.Push(remaining); - CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); + CollectFactorsCpu(primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out bool limitReached); if (limitReached) { while (compositeStack.Count > 0) { - pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false)); + pending.Add(PartialFactorPendingEntry.Rent(compositeStack.Pop(), knownComposite: false)); } } - - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); } else { - pending.Add(new PendingEntry(remaining, knownComposite: false)); + pending.Add(PartialFactorPendingEntry.Rent(remaining, knownComposite: false)); } } @@ -405,13 +402,14 @@ internal static partial class PrimeOrderCalculator int index = 0; for (; index < pendingCount; index++) { - PendingEntry entry = pending[index]; + PartialFactorPendingEntry entry = pending[index]; ulong composite = entry.Value; if (entry.KnownComposite) { cofactor = checked(cofactor * composite); cofactorContainsComposite = true; + PartialFactorPendingEntry.Return(entry); continue; } @@ -425,7 +423,7 @@ internal static partial class PrimeOrderCalculator if (entry.IsPrime) { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); } // composite is never smaller on the execution path // else if (composite > 1UL) @@ -434,6 +432,8 @@ internal static partial class PrimeOrderCalculator cofactor = checked(cofactor * composite); cofactorContainsComposite = true; } + + PartialFactorPendingEntry.Return(entry); } bool cofactorIsPrime; @@ -450,19 +450,6 @@ internal static partial class PrimeOrderCalculator 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) // { @@ -476,41 +463,13 @@ internal static partial class PrimeOrderCalculator // 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); + Array.Sort(result.Factors, result.Exponents, 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); + result.Cofactor = cofactor; + result.FullyFactored = cofactor == 1UL; + result.Count = factorCount; + result.CofactorIsPrime = cofactorIsPrime; return result; } @@ -580,8 +539,9 @@ internal static partial class PrimeOrderCalculator 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) + if (factorization.Count == 0) { + factorization.Dispose(); return false; } @@ -589,45 +549,54 @@ internal static partial class PrimeOrderCalculator { if (factorization.Cofactor <= UInt128.One) { + factorization.Dispose(); return false; } if (!IsPrimeWide(factorization.Cofactor)) { + factorization.Dispose(); return false; } - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.WithAdditionalPrime(factorization.Cofactor); } ReadOnlySpan span = factorization.Factors; int length = factorization.Count; - for (int i = 0; i < length; i++) + try { - UInt128 primeFactor = span[i].Value; - UInt128 reduced = candidate; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) + for (int i = 0; i < length; i++) { - if ((reduced % primeFactor) != UInt128.Zero) + UInt128 primeFactor = span[i].Value; + UInt128 reduced = candidate; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) { - break; - } + if (reduced.ReduceCycleRemainder(primeFactor) != UInt128.Zero) + { + break; + } - reduced /= primeFactor; - if (powUsed >= powBudget && powBudget > 0) - { - return false; - } + reduced /= primeFactor; + if (powUsed >= powBudget && powBudget > 0) + { + return false; + } - powUsed++; - if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) - { - return false; + powUsed++; + if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) + { + return false; + } } } - } - return true; + return true; + } + finally + { + factorization.Dispose(); + } } private static bool TryConfirmOrderWideGpu(PrimeOrderCalculatorAccelerator gpu, in UInt128 prime, in UInt128 order, in MontgomeryDivisorData? divisorData, in PrimeOrderCalculatorConfig config) @@ -643,8 +612,9 @@ internal static partial class PrimeOrderCalculator } PartialFactorResult128 factorization = PartialFactorWide(order, config); - if (factorization.Factors is null) + if (factorization.Count == 0) { + factorization.Dispose(); return false; } @@ -652,15 +622,17 @@ internal static partial class PrimeOrderCalculator { if (factorization.Cofactor <= UInt128.One) { + factorization.Dispose(); return false; } if (!IsPrimeWide(factorization.Cofactor)) { + factorization.Dispose(); return false; } - factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.WithAdditionalPrime(factorization.Cofactor); } ReadOnlySpan span = factorization.Factors; @@ -679,11 +651,13 @@ internal static partial class PrimeOrderCalculator reduced /= primeFactor; if (Pow2ModWideGpu(gpu, reduced, prime, divisorData) == UInt128.One) { + factorization.Dispose(); return false; } } } + factorization.Dispose(); return true; } @@ -695,7 +669,6 @@ internal static partial class PrimeOrderCalculator ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, - PartialFactorResult phiFactors, PartialFactorResult? cachedOrderFactors, out ulong result) { @@ -707,10 +680,10 @@ internal static partial class PrimeOrderCalculator } // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, order, config); + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorGpu(gpu, divisorData, order, config); try { - if (orderFactors.Factors is null) + if (orderFactors.Count == 0) { return false; } @@ -726,16 +699,13 @@ internal static partial class PrimeOrderCalculator } int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); + List candidates = orderFactors.FactorCandidatesList; candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); + BuildCandidates(order, orderFactors, orderFactors.Count, candidates, capacity); int candidateCount = candidates.Count; if (candidateCount == 0) { - ThreadStaticPools.ReturnUlongList(candidates); return false; } @@ -745,11 +715,12 @@ internal static partial class PrimeOrderCalculator int powUsed = 0; bool allowGpuBatch = true; Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); + ExponentRemainderStepperCpu powStepper = orderFactors.ExponentRemainderStepper; + powStepper.Reset(); bool powStepperInitialized = false; int index = 0; - Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; + Span stackGpuRemainders = new(orderFactors.SpecialMaxBuffer, 0, PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount); ArrayPool pool = ThreadStaticPools.UlongPool; while (index < candidateCount && powUsed < powBudget) { @@ -819,9 +790,6 @@ internal static partial class PrimeOrderCalculator pool.Return(gpuPool, clearArray: false); } - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); result = candidate; return true; } @@ -859,9 +827,6 @@ internal static partial class PrimeOrderCalculator pool.Return(gpuPool, clearArray: false); } - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); result = candidate; return true; } @@ -874,9 +839,6 @@ internal static partial class PrimeOrderCalculator index += batchSize; } - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); // DebugLog("No candidate confirmed"); return false; } @@ -902,66 +864,75 @@ internal static partial class PrimeOrderCalculator } PartialFactorResult128 orderFactors = PartialFactorWide(order, config); - if (orderFactors.Factors is null) + if (orderFactors.Count == 0) { + orderFactors.Dispose(); return false; } - if (!orderFactors.FullyFactored) + try { - if (orderFactors.Cofactor <= UInt128.One) + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= UInt128.One) + { + return false; + } + + if (!IsPrimeWide(orderFactors.Cofactor)) + { + return false; + } + + orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; + List candidates = orderFactors.FactorCandidatesList; + candidates.Clear(); + FactorEntry128[] factorArray = orderFactors.Factors!; + BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity); + if (candidates.Count == 0) { return false; } - if (!IsPrimeWide(orderFactors.Cofactor)) + SortCandidatesWide(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; + int powUsed = 0; + int candidateCount = candidates.Count; + + for (int i = 0; i < candidateCount; i++) { - return false; + if (powUsed >= powBudget && powBudget > 0) + { + break; + } + + UInt128 candidate = candidates[i]; + powUsed++; + + if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) != UInt128.One) + { + continue; + } + + if (!TryConfirmCandidateWideGpu(gpu, prime, candidate, config, ref powUsed, powBudget, divisorData)) + { + continue; + } + + result = candidate; + return true; } - orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = new(capacity); - FactorEntry128[] factorArray = orderFactors.Factors!; - BuildCandidatesWide(order, factorArray, orderFactors.Count, candidates, capacity); - if (candidates.Count == 0) - { return false; } - - SortCandidatesWide(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; - int powUsed = 0; - int candidateCount = candidates.Count; - - for (int i = 0; i < candidateCount; i++) + finally { - if (powUsed >= powBudget && powBudget > 0) - { - break; - } - - UInt128 candidate = candidates[i]; - powUsed++; - - if (Pow2ModWideGpu(gpu, candidate, prime, divisorData) != UInt128.One) - { - continue; - } - - if (!TryConfirmCandidateWideGpu(gpu, prime, candidate, config, ref powUsed, powBudget, divisorData)) - { - continue; - } - - result = candidate; - return true; + orderFactors.Dispose(); } - - return false; } private static bool TryPollardRhoGpu(PrimeOrderCalculatorAccelerator gpu, ulong n, out ulong factor) @@ -1077,4 +1048,3 @@ internal static partial class PrimeOrderCalculator return true; } } - diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs index 87f2dae..979d12d 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Hybrid.cs @@ -1,8 +1,6 @@ -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; @@ -25,12 +23,12 @@ internal static partial class PrimeOrderCalculator ulong phi = prime - 1UL; - PartialFactorResult phiFactors = PartialFactorHybrid(gpu, phi, config); + PartialFactorResult phiFactors = PartialFactorHybrid(gpu, divisorData, phi, config); ulong result; - if (phiFactors.Factors is null) + if (phiFactors.Count == 0) { - result = CalculateByFactorizationCpu(prime, divisorData); + result = CalculateByFactorizationCpu(prime, divisorData, phiFactors); phiFactors.Dispose(); return result; @@ -75,15 +73,14 @@ internal static partial class PrimeOrderCalculator } else { - divisorData = MontgomeryDivisorData.Empty; - result = CalculateWideInternalCpu(prime, previousOrder, divisorData, config); + result = CalculateWideInternalCpu(prime, previousOrder, config); } return result; } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span primeSlots, Span exponentSlots, ref int factorCount, List pending, FixedCapacityStack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) + private static void CollectFactorsHybrid(PrimeOrderCalculatorAccelerator gpu, Span primeSlots, Span exponentSlots, ref int factorCount, List pending, FixedCapacityStack compositeStack, long deadlineTimestamp, out bool pollardRhoDeadlineReached) { while (compositeStack.Count > 0) { @@ -92,20 +89,20 @@ internal static partial class PrimeOrderCalculator if (isPrime) { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); continue; } pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; if (pollardRhoDeadlineReached) { - pending.Add(new PendingEntry(composite, knownComposite: true)); + pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true)); continue; } if (!TryPollardRhoCpu(composite, deadlineTimestamp, out ulong factor)) { - pending.Add(new PendingEntry(composite, knownComposite: true)); + pending.Add(new PartialFactorPendingEntry(composite, knownComposite: true)); continue; } @@ -134,7 +131,7 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, ulong value, in PrimeOrderCalculatorConfig config) + private static PartialFactorResult PartialFactorHybrid(PrimeOrderCalculatorAccelerator gpu, in MontgomeryDivisorData divisorData, ulong value, in PrimeOrderCalculatorConfig config) { if (value <= 1UL) { @@ -142,18 +139,17 @@ internal static partial class PrimeOrderCalculator } // stackalloc is faster than pooling - Span primeSlots = stackalloc ulong[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; - Span exponentSlots = stackalloc int[PrimeOrderConstants.GpuSmallPrimeFactorSlots]; - ulong[] primeSlotArray = []; - int[] exponentSlotArray = []; + PartialFactorResult result = PartialFactorResult.Rent(divisorData); + Span primeSlots = new (result.Factors); + Span exponentSlots = new (result.Exponents); // 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; + uint limit = config.SmallFactorLimit; - List pending = ThreadStaticPools.RentPrimeOrderPendingEntryList(2); + List pending = result.Pending; + pending.Clear(); PopulateSmallPrimeFactorsCpu( value, @@ -161,33 +157,30 @@ internal static partial class PrimeOrderCalculator primeSlots, exponentSlots, out int factorCount, - out remaining); + out ulong remaining); - bool limitReached; if (remaining > 1UL) { if (config.PollardRhoMilliseconds > 0) { long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - FixedCapacityStack compositeStack = ThreadStaticPools.RentUlongStack(4); + FixedCapacityStack compositeStack = result.CompositeStack; + compositeStack.Clear(); compositeStack.Push(remaining); - CollectFactorsHybrid(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out limitReached); + CollectFactorsHybrid(gpu, primeSlots, exponentSlots, ref factorCount, pending, compositeStack, deadlineTimestamp, out bool limitReached); if (limitReached) { while (compositeStack.Count > 0) { - pending.Add(new PendingEntry(compositeStack.Pop(), knownComposite: false)); + pending.Add(PartialFactorPendingEntry.Rent(compositeStack.Pop(), knownComposite: false)); } } - - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); } else { - pending.Add(new PendingEntry(remaining, knownComposite: false)); + pending.Add(PartialFactorPendingEntry.Rent(remaining, knownComposite: false)); } } @@ -197,13 +190,14 @@ internal static partial class PrimeOrderCalculator int index = 0; for (; index < pendingCount; index++) { - PendingEntry entry = pending[index]; + PartialFactorPendingEntry entry = pending[index]; ulong composite = entry.Value; if (entry.KnownComposite) { cofactor = checked(cofactor * composite); cofactorContainsComposite = true; + PartialFactorPendingEntry.Return(entry); continue; } @@ -217,7 +211,7 @@ internal static partial class PrimeOrderCalculator if (entry.IsPrime) { - AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite, 1); + AddFactorToCollector(primeSlots, exponentSlots, ref factorCount, composite); } // composite is never smaller on the execution path // else if (composite > 1UL) @@ -226,6 +220,8 @@ internal static partial class PrimeOrderCalculator cofactor = checked(cofactor * composite); cofactorContainsComposite = true; } + + PartialFactorPendingEntry.Return(entry); } bool cofactorIsPrime; @@ -242,20 +238,16 @@ internal static partial class PrimeOrderCalculator 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) + // { + // result.Cofactor = cofactor; + // result.FullyFactored = cofactor != value && cofactor == 1UL; + // result.Count = 0; + // result.CofactorIsPrime = cofactorIsPrime; + // goto ReturnResult; + // } + // if (factorCount == 0) // { // if (cofactor == value) @@ -268,41 +260,13 @@ internal static partial class PrimeOrderCalculator // 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); + Array.Sort(result.Factors, result.Exponents, 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); + result.Cofactor = cofactor; + result.FullyFactored = cofactor == 1UL; + result.Count = factorCount; + result.CofactorIsPrime = cofactorIsPrime; return result; } @@ -336,13 +300,13 @@ internal static partial class PrimeOrderCalculator ulong phi, PartialFactorResult phiFactors) { - if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, prime, phiFactors, divisorData)) + if (phiFactors.FullyFactored && TrySpecialMaxHybrid(gpu, phi, phiFactors, divisorData)) { return phi; } ulong candidateOrder = InitializeStartingOrderHybrid(gpu, prime, phi, divisorData); - candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); + candidateOrder = ExponentLoweringCpu(candidateOrder, phiFactors, divisorData); if (TryConfirmOrderHybrid(gpu, prime, candidateOrder, divisorData, config, out PartialFactorResult? orderFactors)) { @@ -352,10 +316,10 @@ internal static partial class PrimeOrderCalculator if (config.StrictMode) { orderFactors?.Dispose(); - return CalculateByFactorizationCpu(prime, divisorData); + return CalculateByFactorizationCpu(prime, divisorData, phiFactors); } - if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) + if (TryHeuristicFinishHybrid(gpu, prime, candidateOrder, previousOrder, divisorData, config, orderFactors, out ulong order)) { return order; } @@ -407,9 +371,9 @@ internal static partial class PrimeOrderCalculator [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); + PartialFactorResult factorization = PartialFactorHybrid(gpu, divisorData, candidate, config); - if (factorization.Factors is null) + if (factorization.Count == 0) { factorization.Dispose(); return false; @@ -430,12 +394,8 @@ internal static partial class PrimeOrderCalculator 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); + Span stackBuffer = new(factorization.SpecialMaxBuffer); + ExponentRemainderStepperCpu stepper = factorization.ExponentRemainderStepper; bool violates = false; for (int i = 0; i < length; i++) @@ -449,35 +409,13 @@ internal static partial class PrimeOrderCalculator // 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)) + if (CheckCandidateViolation(stackBuffer[..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; @@ -504,8 +442,8 @@ internal static partial class PrimeOrderCalculator return false; } - PartialFactorResult factorization = PartialFactorHybrid(gpu, order, config); - if (factorization.Factors is null) + PartialFactorResult factorization = PartialFactorHybrid(gpu, divisorData, order, config); + if (factorization.Count == 0) { factorization.Dispose(); return false; @@ -522,7 +460,7 @@ internal static partial class PrimeOrderCalculator factorization.WithAdditionalPrime(factorization.Cofactor); } - if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization)) + if (!ValidateOrderAgainstFactors(order, factorization)) { reusableFactorization = factorization; return false; @@ -540,7 +478,6 @@ internal static partial class PrimeOrderCalculator ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderCalculatorConfig config, - PartialFactorResult phiFactors, PartialFactorResult? cachedOrderFactors, out ulong result) { @@ -552,10 +489,10 @@ internal static partial class PrimeOrderCalculator } // Reuse the partial factorization from TryConfirmOrderCpu when available. - PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, order, config); + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactorHybrid(gpu, divisorData, order, config); try { - if (orderFactors.Factors is null) + if (orderFactors.Count == 0) { return false; } @@ -571,16 +508,13 @@ internal static partial class PrimeOrderCalculator } int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecksCapacity; - List candidates = ThreadStaticPools.RentUlongList(capacity); + List candidates = orderFactors.FactorCandidatesList; candidates.Clear(); - ulong[] factorArray = orderFactors.Factors!; - int[] exponentsArray = orderFactors.Exponents!; - BuildCandidates(order, factorArray, exponentsArray, orderFactors.Count, candidates, capacity); + BuildCandidates(order, orderFactors, orderFactors.Count, candidates, capacity); int candidateCount = candidates.Count; if (candidateCount == 0) { - ThreadStaticPools.ReturnUlongList(candidates); return false; } @@ -589,74 +523,46 @@ internal static partial class PrimeOrderCalculator int powBudget = config.MaxPowChecks <= 0 ? candidateCount : config.MaxPowChecks; int powUsed = 0; Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - ExponentRemainderStepperCpu powStepper = ThreadStaticPools.RentExponentStepperCpu(divisorData); - bool powStepperInitialized = false; + ExponentRemainderStepperCpu powStepper = orderFactors.ExponentRemainderStepper; + powStepper.Reset(); - int index = 0; - Span stackGpuRemainders = stackalloc ulong[PerfectNumberConstants.DefaultSmallPrimeFactorSlotCount]; - ArrayPool pool = ThreadStaticPools.UlongPool; + ulong candidate = candidateSpan[0]; + powUsed++; + + bool equalsOne = powStepper.InitializeCpuIsUnity(candidate); + if (equalsOne && TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + result = candidate; + return true; + } + + int index = 1; 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) + if (budgetRemaining <= 0) { break; } - ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); - int startIndex = 0; + candidate = candidateSpan[index++]; + powUsed++; - if (!powStepperInitialized) + equalsOne = powStepper.ComputeNextIsUnity(candidate); + if (!equalsOne) { - 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; + continue; } - for (int i = startIndex; i < batchSize && powUsed < powBudget; i++) + if (!TryConfirmCandidateHybrid(gpu, prime, candidate, divisorData, config, ref powUsed, powBudget)) { - 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; + continue; } - index += batchSize; + result = candidate; + return true; } - candidates.Clear(); - ThreadStaticPools.ReturnUlongList(candidates); - ThreadStaticPools.ReturnExponentStepperCpu(powStepper); - return false; } finally @@ -666,7 +572,7 @@ internal static partial class PrimeOrderCalculator } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + private static bool TrySpecialMaxHybrid(PrimeOrderCalculatorAccelerator gpu, ulong phi, 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, @@ -676,15 +582,13 @@ internal static partial class PrimeOrderCalculator // 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 EvaluateSpecialMaxCandidatesCpu(factors, phi); } - return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, prime, divisorData); + ReadOnlySpan factorSpan = new(factors.Factors, 0, length); + return EvaluateSpecialMaxCandidatesGpu(gpu, factorSpan, phi, divisorData); } } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs index 3681aec..23d037a 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -1,64 +1,90 @@ using System.Buffers; using System.Diagnostics; using System.Numerics; +using System.Runtime.CompilerServices; namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private static void SortCandidatesWide(UInt128 prime, in UInt128? previousOrder, List candidates) { UInt128 previous = previousOrder ?? UInt128.Zero; - int previousGroup = previousOrder.HasValue ? GetGroupWide(previousOrder.Value, prime) : 1; bool hasPrevious = previousOrder.HasValue; + int previousGroup = hasPrevious ? GetGroupWide(previous, prime) : 1; - candidates.Sort((x, y) => - { - CandidateKey128 keyX = BuildKeyWide(x, prime, previous, previousGroup, hasPrevious); - CandidateKey128 keyY = BuildKeyWide(y, prime, previous, previousGroup, hasPrevious); - - int primary = keyX.Primary.CompareTo(keyY.Primary); - if (primary != 0) - { - return primary; - } - - int secondary = CompareComponents(keyX.SecondaryDescending, keyX.Secondary, keyY.SecondaryDescending, keyY.Secondary); - if (secondary != 0) - { - return secondary; - } - - return CompareComponents(keyX.TertiaryDescending, keyX.Tertiary, keyY.TertiaryDescending, keyY.Tertiary); - }); + candidates.Sort((left, right) => CompareCandidatesWide(prime, left, right, previous, hasPrevious, previousGroup)); } - private static CandidateKey128 BuildKeyWide( - in UInt128 value, - in UInt128 prime, - in UInt128 previous, - int previousGroup, - bool hasPrevious) - { - int group = GetGroupWide(value, prime); - if (group == 0) - { - return new CandidateKey128(int.MaxValue, false, UInt128.Zero, false, UInt128.Zero); - } + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static int CompareCandidatesWide(UInt128 prime, UInt128 left, UInt128 right, UInt128 previous, bool hasPrevious, int previousGroup) + { + int leftGroup = GetGroupWide(left, prime); + int rightGroup = GetGroupWide(right, prime); + if (leftGroup == 0) + { + return rightGroup == 0 ? 0 : 1; + } - bool isGe = !hasPrevious || value >= previous; - int primary = ComputePrimary(group, isGe, previousGroup); + if (rightGroup == 0) + { + return -1; + } - if (group == 3) - { - return new CandidateKey128(primary, true, value, true, value); - } + bool leftIsGe = !hasPrevious || left >= previous; + bool rightIsGe = !hasPrevious || right >= previous; + int leftPrimary = ComputePrimary(leftGroup, leftIsGe, previousGroup); + int rightPrimary = ComputePrimary(rightGroup, rightIsGe, previousGroup); - UInt128 distance = hasPrevious ? (value > previous ? value - previous : previous - value) : value; - return new CandidateKey128(primary, false, distance, false, value); - } + int compare = leftPrimary.CompareTo(rightPrimary); + if (compare != 0) + { + return compare; + } - private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY) + bool leftDescending = leftGroup == 3; + bool rightDescending = rightGroup == 3; + + UInt128 leftSecondary; + UInt128 leftTertiary; + if (leftDescending) + { + leftSecondary = left; + leftTertiary = left; + } + else + { + UInt128 leftReference = hasPrevious ? (left > previous ? left - previous : previous - left) : left; + leftSecondary = leftReference; + leftTertiary = left; + } + + UInt128 rightSecondary; + UInt128 rightTertiary; + if (rightDescending) + { + rightSecondary = right; + rightTertiary = right; + } + else + { + UInt128 rightReference = hasPrevious ? (right > previous ? right - previous : previous - right) : right; + rightSecondary = rightReference; + rightTertiary = right; + } + + compare = CompareComponents(leftDescending, leftSecondary, rightDescending, rightSecondary); + if (compare != 0) + { + return compare; + } + + return CompareComponents(leftDescending, leftTertiary, rightDescending, rightTertiary); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] + private static int CompareComponents(bool descendingX, in UInt128 valueX, bool descendingY, in UInt128 valueY) { if (descendingX == descendingY) { @@ -83,22 +109,23 @@ internal static partial class PrimeOrderCalculator return valueX > valueY ? 1 : -1; } + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private static int GetGroupWide(in UInt128 value, in UInt128 prime) { - UInt128 threshold1 = prime >> 3; - if (value <= threshold1) + UInt128 threshold = prime >> 3; + if (value <= threshold) { return 1; } - UInt128 threshold2 = prime >> 2; - if (value <= threshold2) + threshold = prime >> 2; + if (value <= threshold) { return 2; } - UInt128 threshold3 = (UInt128)(((BigInteger)prime * 3) >> 3); - if (value <= threshold3) + threshold = (UInt128)(((BigInteger)prime * 3) >> 3); + if (value <= threshold) { return 3; } @@ -106,8 +133,10 @@ internal static partial class PrimeOrderCalculator return 0; } + [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private static void BuildCandidatesWide(in UInt128 order, FactorEntry128[] factors, int count, List candidates, int limit) { + // TODO: Is this condition ever met on EvenPerfectBitScanner execution paths for this method? Can we remove it? if (count == 0) { return; @@ -118,6 +147,7 @@ internal static partial class PrimeOrderCalculator BuildCandidatesRecursiveWide(order, buffer, 0, UInt128.One, candidates, limit); } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static void BuildCandidatesRecursiveWide( in UInt128 order, in ReadOnlySpan factors, @@ -185,33 +215,14 @@ internal static partial class PrimeOrderCalculator return PartialFactorResult128.Empty; } - Dictionary counts = new(capacity: 8); + PartialFactorResult128 result = PartialFactorResult128.Rent(); + Dictionary counts = result.FactorCounts; + counts.Clear(); uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - UInt128 remaining; + UInt128 remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); - // if (IsGpuPow2Allowed && value <= ulong.MaxValue) - // { - // Dictionary narrowCounts = new(capacity: 8); - // if (TryPopulateSmallPrimeFactorsGpu((ulong)value, limit, narrowCounts, out ulong narrowRemaining)) - // { - // foreach (KeyValuePair entry in narrowCounts) - // { - // counts[(UInt128)entry.Key] = entry.Value; - // } - - // remaining = (UInt128)narrowRemaining; - // } - // else - // { - // remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); - // } - // } - // else - // { - remaining = PopulateSmallPrimeFactorsCpuWide(value, limit, counts); - // } - - List pending = new(); + List pending = result.Pending; + pending.Clear(); if (remaining > UInt128.One) { pending.Add(remaining); @@ -220,11 +231,12 @@ internal static partial class PrimeOrderCalculator if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) { long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - Stack stack = new(); + FixedCapacityStack stack = result.CompositeStack; + stack.Clear(); stack.Push(remaining); pending.Clear(); - long timestamp = 0L; // reused for deadline checks. + long timestamp; // reused for deadline checks. while (stack.Count > 0) { UInt128 composite = stack.Pop(); @@ -276,33 +288,13 @@ internal static partial class PrimeOrderCalculator if (counts.Count == 0 && cofactor == value) { - return new PartialFactorResult128(null, value, false, 0); + result.Cofactor = value; + // Everything else stays at default values initiated in .Rent. + return result; } - ArrayPool? pool = null; - FactorEntry128[]? rented = counts.Count > 0 ? (pool = ThreadStaticPools.FactorEntry128Pool).Rent(counts.Count) : null; - int index = 0; - if (rented is not null) - { - foreach (KeyValuePair entry in counts) - { - rented[index] = new FactorEntry128(entry.Key, entry.Value); - index++; - } - } - - FactorEntry128[]? resultArray = null; - if (pool is not null) - { - Span span = rented.AsSpan(0, index); - span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - resultArray = new FactorEntry128[index]; - span.CopyTo(resultArray); - pool.Return(rented!, clearArray: false); - } - - bool fullyFactored = cofactor == UInt128.One; - return new PartialFactorResult128(resultArray, cofactor, fullyFactored, index); + result.InitializeFromCounts(counts, cofactor, cofactor == UInt128.One); + return result; } private static UInt128 PopulateSmallPrimeFactorsCpuWide(in UInt128 value, uint limit, Dictionary counts) @@ -357,7 +349,7 @@ internal static partial class PrimeOrderCalculator UInt128 c = UInt128.One; UInt128 x = (UInt128)2UL; UInt128 y = x; - long timestamp = 0L; // reused for deadline checks. + long timestamp; // reused for deadline checks. while (true) { @@ -368,8 +360,7 @@ internal static partial class PrimeOrderCalculator } x = AdvancePolynomialWide(x, c, n); - y = AdvancePolynomialWide(y, c, n); - y = AdvancePolynomialWide(y, c, n); + y = AdvancePolynomialTwiceWide(y, c, n); UInt128 diff = x > y ? x - y : y - x; UInt128 d = BinaryGcdWide(diff, n); @@ -396,6 +387,14 @@ internal static partial class PrimeOrderCalculator return (UInt128)value; } + private static UInt128 AdvancePolynomialTwiceWide(in UInt128 x, in UInt128 c, in UInt128 modulus) + { + BigInteger value = (BigInteger)x; + value = (value * value + (BigInteger)c) % (BigInteger)modulus; + value = (value * value + (BigInteger)c) % (BigInteger)modulus; + return (UInt128)value; + } + private static UInt128 BinaryGcdWide(UInt128 a, UInt128 b) { if (a == UInt128.Zero) @@ -479,8 +478,7 @@ internal static partial class PrimeOrderCalculator while (d == UInt128.One) { x = AdvancePolynomialWide(x, c, n); - y = AdvancePolynomialWide(y, c, n); - y = AdvancePolynomialWide(y, c, n); + y = AdvancePolynomialTwiceWide(y, c, n); UInt128 diff = x > y ? x - y : y - x; d = BinaryGcdWide(diff, n); @@ -509,12 +507,13 @@ internal static partial class PrimeOrderCalculator } uint[] smallPrimes = PrimesGenerator.SmallPrimes; + ulong[] smallPrimesPow2 = PrimesGenerator.SmallPrimesPow2; int length = smallPrimes.Length; for (int i = 0; i < length; i++) { uint prime = smallPrimes[i]; UInt128 primeValue = prime; - UInt128 primeSquare = primeValue * primeValue; + UInt128 primeSquare = (UInt128)smallPrimesPow2[i]; if (primeSquare > value) { break; @@ -588,44 +587,104 @@ internal static partial class PrimeOrderCalculator } } - private readonly struct CandidateKey128(int primary, bool secondaryDescending, in UInt128 secondary, bool tertiaryDescending, in UInt128 tertiary) + private sealed class PartialFactorResult128 { - public int Primary { get; } = primary; + [ThreadStatic] + private static PartialFactorResult128? s_poolHead; - public bool SecondaryDescending { get; } = secondaryDescending; + public static readonly PartialFactorResult128 Empty = new() + { + Cofactor = UInt128.One, + FullyFactored = true, + Count = 0 + }; - public UInt128 Secondary { get; } = secondary; + private const int ExponentHardLimit = 256; - public bool TertiaryDescending { get; } = tertiaryDescending; + internal readonly FactorEntry128[] TempFactors = new FactorEntry128[ExponentHardLimit]; + internal readonly List FactorCandidatesList = new(ExponentHardLimit); + internal readonly Dictionary FactorCounts = new(capacity: ExponentHardLimit); + internal readonly List Pending = new(capacity: ExponentHardLimit); + internal readonly FixedCapacityStack CompositeStack = new(ExponentHardLimit); - public UInt128 Tertiary { get; } = tertiary; + private PartialFactorResult128? _next; + + public UInt128 Cofactor; + public int Count; + public readonly FactorEntry128[] Factors = new FactorEntry128[ExponentHardLimit]; + public bool FullyFactored; + public bool HasFactors; + + private PartialFactorResult128() + { + } + + public static PartialFactorResult128 Rent() + { + if (s_poolHead is { } instance) + { + s_poolHead = instance._next; + instance._next = null; + + instance.Cofactor = UInt128.One; + instance.Count = 0; + instance.FullyFactored = false; + instance.HasFactors = false; + + return instance; + } + + return new PartialFactorResult128(); + } + + public void WithAdditionalPrime(UInt128 prime) + { + int index = Count; + + FactorEntry128[] buffer = Factors; + buffer[index] = new FactorEntry128(prime, 1); + index++; + Count = index; + + buffer.AsSpan(0, index).Sort(static (a, b) => a.Value.CompareTo(b.Value)); + Cofactor = UInt128.One; + FullyFactored = true; + HasFactors = true; + } + + public void Dispose() + { + if (this == Empty) + { + return; + } + + _next = s_poolHead; + s_poolHead = this; + } + + internal void InitializeFromCounts(Dictionary counts, UInt128 cofactor, bool fullyFactored) + { + int count = counts.Count; + Count = count; + Cofactor = cofactor; + FullyFactored = fullyFactored; + bool hasFactors = count != 0; + HasFactors = hasFactors; + + if (!hasFactors) + { + return; + } + + int index = 0; + foreach (KeyValuePair entry in counts) + { + Factors[index++] = new FactorEntry128(entry.Key, entry.Value); + } + + // We benefit from static (non-virtual) call to Sort. When JIT finds a method call after a constructor, it optimizes it to a true static call. + new Span(Factors, 0, count).Sort(static (a, b) => a.Value.CompareTo(b.Value)); + } } - - private readonly struct PartialFactorResult128(FactorEntry128[]? factors, UInt128 cofactor, bool fullyFactored, int count) - { - public FactorEntry128[]? Factors { get; } = factors; - - public UInt128 Cofactor { get; } = cofactor; - - public bool FullyFactored { get; } = fullyFactored; - - public int Count { get; } = count; - - public static PartialFactorResult128 Empty => new(null, UInt128.One, true, 0); - - public PartialFactorResult128 WithAdditionalPrime(UInt128 prime) - { - if (Factors is null) - { - FactorEntry128[] local = new FactorEntry128[1]; - local[0] = new FactorEntry128(prime, 1); - return new PartialFactorResult128(local, UInt128.One, true, 1); - } - - FactorEntry128[] extended = new FactorEntry128[Count + 1]; - Array.Copy(Factors, extended, Count); - extended[Count] = new FactorEntry128(prime, 1); - return new PartialFactorResult128(extended, UInt128.One, true, Count + 1); - } - } } diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index 598eab4..3efe013 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -263,12 +263,12 @@ namespace PerfectNumbers.Core public static void Return((ulong, ulong)[] array) => _ulongUlongTupleArrayPool!.Return(array, clearArray: false); [ThreadStatic] - private static FixedCapacityStack>? _primeOrderPendingEntryListPool; + private static FixedCapacityStack>? _primeOrderPendingEntryListPool; [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - internal static List RentPrimeOrderPendingEntryList(int capacityHint) + internal static List RentPrimeOrderPendingEntryList(int capacityHint) { - FixedCapacityStack>? pool = _primeOrderPendingEntryListPool; + FixedCapacityStack>? pool = _primeOrderPendingEntryListPool; if (pool is not null) { if (pool.Pop() is { } list) @@ -286,11 +286,11 @@ namespace PerfectNumbers.Core _primeOrderPendingEntryListPool = new(PerfectNumberConstants.DefaultPoolCapacity); } - return new List(capacityHint); + return new List(capacityHint); } [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] - internal static void ReturnPrimeOrderPendingEntryList(List list) => _primeOrderPendingEntryListPool!.Push(list); + internal static void ReturnPrimeOrderPendingEntryList(List list) => _primeOrderPendingEntryListPool!.Push(list); [ThreadStatic] private static FixedCapacityStack>? _ulongStackPool;