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