From ca0f498d0d9f44249880f8d52903cbf19b3685fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Fri, 17 Oct 2025 02:34:41 +0200 Subject: [PATCH] Recycle small prime stepper arrays --- .../PrimeOrderCalculator.Cpu.cs | 1525 +++++++++++++++ .../PrimeOrderCalculator.Gpu.cs | 48 + ...Wide.cs => PrimeOrderCalculator.Shared.cs} | 108 +- PerfectNumbers.Core/PrimeOrderCalculator.cs | 1681 +---------------- PerfectNumbers.Core/ThreadStaticPools.cs | 24 + 5 files changed, 1744 insertions(+), 1642 deletions(-) create mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs create mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs rename PerfectNumbers.Core/{PrimeOrderCalculator.Wide.cs => PrimeOrderCalculator.Shared.cs} (92%) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs new file mode 100644 index 0000000..b537f3a --- /dev/null +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -0,0 +1,1525 @@ +using System; +using System.Buffers; +using System.Buffers.Binary; +using System.Collections.Generic; +using System.Diagnostics; +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; +using System.Threading; +using PerfectNumbers.Core.Gpu; + +namespace PerfectNumbers.Core; + +internal static partial class PrimeOrderCalculator +{ + [ThreadStatic] + private static CycleRemainderStepper[]? s_smallPrimeCycleSteppers; + + [ThreadStatic] + private static ulong[]? s_smallPrimeStepperPreviousValues; + + [ThreadStatic] + private static ulong[]? s_smallPrimeStepperRemainders; + + [ThreadStatic] + private static bool[]? s_smallPrimeStepperHasState; + + private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) + { + // TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31. + if (prime <= 3UL) + { + return prime == 3UL ? 2UL : 1UL; + } + + ulong phi = prime - 1UL; + + if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) + { + return gpuOrder; + } + + PartialFactorResult phiFactors = PartialFactor(phi, config); + + ulong result; + if (phiFactors.Factors is null) + { + result = CalculateByFactorizationCpu(prime, divisorData); + phiFactors.Dispose(); + return result; + } + + result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); + phiFactors.Dispose(); + return result; + } + + private static ulong RunHeuristicPipelineCpu( + ulong prime, + ulong? previousOrder, + in PrimeOrderSearchConfig config, + in MontgomeryDivisorData divisorData, + ulong phi, + PartialFactorResult phiFactors) + { + if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) + { + return phi; + } + + ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); + candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); + + if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config)) + { + return candidateOrder; + } + + if (config.Mode == PrimeOrderMode.Strict) + { + return CalculateByFactorizationCpu(prime, divisorData); + } + + if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) + { + return order; + } + + return candidateOrder; + } + + + private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + { + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + for (int i = 0; i < length; i++) + { + ulong factor = factorSpan[i].Value; + ulong reduced = phi / factor; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + return false; + } + } + + return true; + } + + private static ulong InitializeStartingOrderCpu(ulong prime, ulong phi, in MontgomeryDivisorData divisorData) + { + ulong order = phi; + if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) + { + ulong half = phi >> 1; + if (Pow2EqualsOneCpu(half, prime, divisorData)) + { + order = half; + } + } + + return order; + } + + private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + { + ArrayPool pool = ThreadStaticPools.FactorEntryPool; + + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + FactorEntry[] tempArray = pool.Rent(length + 1); + + Span buffer = tempArray; + factorSpan.CopyTo(buffer); + + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None); + if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime) + { + buffer[length] = new FactorEntry(factors.Cofactor, 1); + length++; + } + + buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value)); + + for (int i = 0; i < length; i++) + { + ulong primeFactor = buffer[i].Value; + int exponent = buffer[i].Exponent; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != 0UL) + { + break; + } + + ulong reduced = order / primeFactor; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + order = reduced; + continue; + } + + break; + } + } + + pool.Return(tempArray, clearArray: false); + return order; + } + + private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) + { + if (order == 0UL) + { + return false; + } + + // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`. + // DebugLog("Verifying a^order ≡ 1 (mod p)"); + if (!Pow2EqualsOneCpu(order, prime, divisorData)) + { + return false; + } + + // DebugLog("Partial factoring order"); + + // TODO: Do we do partial factoring of order multiple times? + PartialFactorResult factorization = PartialFactor(order, config); + try + { + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL) + { + return false; + } + + // DebugLog("Cofactor > 1. Testing primality of cofactor"); + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); + if (!isPrime) + { + return false; + } + + // DebugLog("Adding cofactor as prime factor"); + PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.Dispose(); + factorization = extended; + } + + ReadOnlySpan span = factorization.Factors; + // DebugLog("Verifying prime-power reductions"); + int length = factorization.Count; + for (int i = 0; i < length; i++) + { + ulong primeFactor = span[i].Value; + ulong reduced = order; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + return false; + } + } + } + + return true; + } + finally + { + factorization.Dispose(); + } + } + + private static bool TryHeuristicFinishCpu( + ulong prime, + ulong order, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderSearchConfig config, + PartialFactorResult phiFactors, + out ulong result) + { + result = 0UL; + if (order <= 1UL) + { + return false; + } + + // TODO: Do we do partial factoring of order multiple times? + // DebugLog("Trying heuristic. Partial factoring order"); + PartialFactorResult orderFactors = PartialFactor(order, config); + try + { + if (orderFactors.Factors is null) + { + return false; + } + + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= 1UL) + { + return false; + } + + var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor); + // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None)) + if (!isPrime) + { + return false; + } + + PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + orderFactors.Dispose(); + orderFactors = extended; + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; + List candidates = ThreadStaticPools.RentUlongList(capacity); + FactorEntry[] factorArray = orderFactors.Factors!; + // DebugLog("Building candidates list"); + BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); + if (candidates.Count == 0) + { + ThreadStaticPools.ReturnUlongList(candidates); + return false; + } + + // DebugLog("Sorting candidates"); + SortCandidates(prime, previousOrder, candidates); + + int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; + int powUsed = 0; + int candidateCount = candidates.Count; + bool allowGpuBatch = true; + Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + + // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); + int index = 0; + const int MaxGpuBatchSize = 256; + const int StackGpuBatchSize = 64; + Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize]; + ArrayPool pool = ThreadStaticPools.UlongPool; + while (index < candidateCount && powUsed < powBudget) + { + int remaining = candidateCount - index; + int budgetRemaining = powBudget - powUsed; + int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize)); + if (batchSize <= 0) + { + break; + } + + ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); + ulong[]? gpuPool = null; + Span pooledGpuRemainders = default; + bool gpuSuccess = false; + bool gpuStackRemainders = false; + GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; + + if (allowGpuBatch && IsGpuPow2Allowed) + { + if (batchSize <= StackGpuBatchSize) + { + Span localRemainders = stackGpuRemainders.Slice(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) + { + gpuSuccess = true; + gpuStackRemainders = true; + } + } + else + { + gpuPool = pool.Rent(batchSize); + Span pooledRemainders = gpuPool.AsSpan(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) + { + pooledGpuRemainders = pooledRemainders; + gpuSuccess = true; + } + else + { + pool.Return(gpuPool, clearArray: false); + gpuPool = null; + } + } + + if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) + { + allowGpuBatch = false; + } + } + + for (int i = 0; i < batchSize && powUsed < powBudget; i++) + { + ulong candidate = batch[i]; + powUsed++; + + bool equalsOne; + if (gpuSuccess) + { + ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; + equalsOne = remainderValue == 1UL; + } + else + { + equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData); + } + if (!equalsOne) + { + continue; + } + + if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + continue; + } + + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + ThreadStaticPools.ReturnUlongList(candidates); + result = candidate; + return true; + } + + if (gpuPool is not null) + { + pool.Return(gpuPool, clearArray: false); + } + + index += batchSize; + } + + ThreadStaticPools.ReturnUlongList(candidates); + // DebugLog("No candidate confirmed"); + return false; + } + finally + { + orderFactors.Dispose(); + } + } + + private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) + { + ulong previous = previousOrder ?? 0UL; + int previousGroup = previousOrder.HasValue ? GetGroup(previousOrder.Value, prime) : 1; + bool hasPrevious = previousOrder.HasValue; + + candidates.Sort((x, y) => + { + CandidateKey keyX = BuildKey(x, prime, previous, previousGroup, hasPrevious); + CandidateKey keyY = BuildKey(y, prime, previous, previousGroup, hasPrevious); + int compare = keyX.Primary.CompareTo(keyY.Primary); + if (compare != 0) + { + return compare; + } + + compare = keyX.Secondary.CompareTo(keyY.Secondary); + if (compare != 0) + { + return compare; + } + + return keyX.Tertiary.CompareTo(keyY.Tertiary); + }); + } + + private static CandidateKey BuildKey(ulong value, ulong prime, ulong previous, int previousGroup, bool hasPrevious) + { + int group = GetGroup(value, prime); + if (group == 0) + { + return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); + } + + bool isGe = !hasPrevious || value >= previous; + int primary = ComputePrimary(group, isGe, previousGroup); + long secondary; + long tertiary; + + if (group == 3) + { + secondary = -(long)value; + tertiary = -(long)value; + } + else + { + ulong reference = hasPrevious ? previous : 0UL; + ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; + secondary = (long)distance; + tertiary = (long)value; + } + + return new CandidateKey(primary, secondary, tertiary); + } + + private static int ComputePrimary(int group, bool isGe, int previousGroup) + { + int groupOffset = group switch + { + 1 => 0, + 2 => 2, + 3 => 4, + _ => 6, + }; + + if (group == previousGroup) + { + if (group == 3) + { + return groupOffset + (isGe ? 0 : 3); + } + + return groupOffset + (isGe ? 0 : 1); + } + + return groupOffset + (isGe ? 0 : 1); + } + + private static int GetGroup(ulong value, ulong prime) + { + ulong threshold1 = prime >> 3; + if (value <= threshold1) + { + return 1; + } + + ulong threshold2 = prime >> 2; + if (value <= threshold2) + { + return 2; + } + + ulong threshold3 = (prime * 3UL) >> 3; + if (value <= threshold3) + { + return 3; + } + + return 0; + } + + private static void BuildCandidates(ulong order, FactorEntry[] factors, int count, List candidates, int limit) + { + if (count == 0) + { + return; + } + + Span buffer = factors.AsSpan(0, count); + buffer.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + BuildCandidatesRecursive(order, buffer, 0, 1UL, candidates, limit); + } + + private static void BuildCandidatesRecursive(ulong order, in ReadOnlySpan factors, int index, ulong divisorProduct, List candidates, int limit) + { + if (candidates.Count >= limit) + { + return; + } + + if (index >= factors.Length) + { + if (divisorProduct == 1UL || divisorProduct == order) + { + return; + } + + ulong candidate = order / divisorProduct; + if (candidate > 1UL && candidate < order) + { + candidates.Add(candidate); + } + + return; + } + + FactorEntry factor = factors[index]; + ulong primeFactor = factor.Value; + ulong contribution = 1UL; + for (int exponent = 0; exponent <= factor.Exponent; exponent++) + { + ulong nextDivisor = divisorProduct * contribution; + if (nextDivisor > order) + { + break; + } + + BuildCandidatesRecursive(order, factors, index + 1, nextDivisor, candidates, limit); + if (candidates.Count >= limit) + { + return; + } + + if (exponent == factor.Exponent) + { + break; + } + + if (contribution > order / primeFactor) + { + break; + } + + contribution *= primeFactor; + } + } + + private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) + { + PartialFactorResult factorization = PartialFactor(candidate, config); + try + { + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL) + { + return false; + } + + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); + if (!isPrime) + { + return false; + } + + PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.Dispose(); + factorization = extended; + } + + ReadOnlySpan span = factorization.Factors; + int length = factorization.Count; + for (int i = 0; i < length; i++) + { + ulong primeFactor = span[i].Value; + ulong reduced = candidate; + for (int iteration = 0; iteration < span[i].Exponent; iteration++) + { + if ((reduced % primeFactor) != 0UL) + { + break; + } + + reduced /= primeFactor; + if (powUsed >= powBudget && powBudget > 0) + { + return false; + } + + powUsed++; + if (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + return false; + } + } + } + + return true; + } + finally + { + factorization.Dispose(); + } + } + + private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) + { + if (IsGpuPow2Allowed) + { + GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); + if (status == GpuPow2ModStatus.Success) + { + return remainder == 1UL; + } + } + + return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; + } + + [ThreadStatic] + private static Dictionary? _primalityCache; + + private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } + + const int FactorSlotCount = GpuSmallPrimeFactorSlots; + + // stackalloc is faster than pooling + // ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount); + // int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); + // Span primeSlots = primeSlotsArray.AsSpan(0, FactorSlotCount); + // Span exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount); + Span primeSlots = stackalloc ulong[FactorSlotCount]; + Span exponentSlots = stackalloc int[FactorSlotCount]; + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); + + int factorCount = 0; + Dictionary? counts = null; + bool useDictionary = false; + + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + ulong remaining = value; + + bool gpuFactored = false; + if (IsGpuHeuristicDevice) + { + gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); + } + + List pending = ThreadStaticPools.RentUlongList(2); + Stack? compositeStack = null; + PartialFactorResult result; + var cache = _primalityCache ??= []; + + if (!gpuFactored) + { + // We don't need to worry about leftovers, because we always use indexes within the calculated counts + // primeSlots.Clear(); + // exponentSlots.Clear(); + if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) + { + useDictionary = true; + counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(factorCount, 8)); + remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); + } + } + + if (remaining > 1UL) + { + pending.Add(remaining); + } + + if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) + { + long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); + compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count << 1, 4)); + compositeStack.Push(remaining); + pending.Clear(); + + bool pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; + if (!pollardRhoDeadlineReached) + { + while (compositeStack.Count > 0) + { + ulong composite = compositeStack.Pop(); + // composite will never be smaller than 1 on the execution path + // if (composite <= 1UL) + // { + // continue; + // } + + if (!cache.TryGetValue(composite, out bool isPrime)) + { + isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + cache.Add(composite, isPrime); + } + + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + continue; + } + + if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) + { + pending.Add(composite); + continue; + } + + ulong quotient = composite / factor; + compositeStack.Push(factor); + compositeStack.Push(quotient); + } + + pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; + } + + if (pollardRhoDeadlineReached) + { + while (compositeStack.Count > 0) + { + pending.Add(compositeStack.Pop()); + } + } + } + + ulong cofactor = 1UL; + int pendingCount = pending.Count; + for (int i = 0; i < pendingCount; i++) + { + ulong composite = pending[i]; + if (!cache.TryGetValue(composite, out bool isPrime)) + { + isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); + cache.Add(composite, isPrime); + } + + if (isPrime) + { + AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); + } + // composite is never smaller on the execution path + // else if (composite > 1UL) + else + { + cofactor = checked(cofactor * composite); + } + } + + ArrayPool pool = ThreadStaticPools.FactorEntryPool; + if (useDictionary) + { + if ((counts is null || counts.Count == 0) && cofactor == value) + { + result = PartialFactorResult.Rent(null, value, false, 0); + goto ReturnResult; + } + + if (counts is null || counts.Count == 0) + { + result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + goto ReturnResult; + } + + FactorEntry[] factors = pool.Rent(counts.Count); + int index = 0; + foreach (KeyValuePair entry in counts) + { + factors[index] = new FactorEntry(entry.Key, entry.Value); + index++; + } + + Array.Sort(factors, static (a, b) => a.Value.CompareTo(b.Value)); + + bool fullyFactored = cofactor == 1UL; + result = PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); + goto ReturnResult; + } + + // This will never happen in production code. We'll always get at least 1 factor + // if (factorCount == 0) + // { + // if (cofactor == value) + // { + // result = PartialFactorResult.Rent(null, value, false, 0); + // goto ReturnResult; + // } + + // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); + // goto ReturnResult; + // } + + FactorEntry[] array = pool.Rent(factorCount); + int arrayIndex = 0; + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeSlots[i]; + int exponentValue = 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[arrayIndex] = new FactorEntry(primeValue, exponentValue); + arrayIndex++; + } + + Span arraySpan = array.AsSpan(0, arrayIndex); + arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + bool fullyFactoredArray = cofactor == 1UL; + result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); + + ReturnResult: + // ThreadStaticPools.UlongPool.Return(primeSlotsArray); + // ThreadStaticPools.IntPool.Return(exponentSlotsArray); + + if (counts is not null) + { + ThreadStaticPools.ReturnUlongIntDictionary(counts); + } + + if (compositeStack is not null) + { + compositeStack.Clear(); + ThreadStaticPools.ReturnUlongStack(compositeStack); + } + + pending.Clear(); + ThreadStaticPools.ReturnUlongList(pending); + return result; + } + + private static CycleRemainderStepper[] EnsureSmallPrimeCycleSteppers(uint[] primes) + { + int primeCount = primes.Length; + CycleRemainderStepper[]? steppers = s_smallPrimeCycleSteppers; + ArrayPool stepperPool = ThreadStaticPools.CycleRemainderStepperPool; + ArrayPool ulongPool = ThreadStaticPools.UlongPool; + ArrayPool boolPool = ThreadStaticPools.BoolPool; + if (steppers is null || steppers.Length < primeCount) + { + CycleRemainderStepper[]? previousSteppers = steppers; + steppers = stepperPool.Rent(primeCount); + int stepperLength = steppers.Length; + for (int i = 0; i < primeCount; i++) + { + steppers[i] = new CycleRemainderStepper(primes[i]); + } + for (int i = primeCount; i < stepperLength; i++) + { + steppers[i] = default; + } + + if (previousSteppers is not null) + { + stepperPool.Return(previousSteppers, clearArray: false); + } + + if (s_smallPrimeStepperPreviousValues is not null) + { + ulongPool.Return(s_smallPrimeStepperPreviousValues, clearArray: false); + } + + if (s_smallPrimeStepperRemainders is not null) + { + ulongPool.Return(s_smallPrimeStepperRemainders, clearArray: false); + } + + if (s_smallPrimeStepperHasState is not null) + { + boolPool.Return(s_smallPrimeStepperHasState, clearArray: false); + } + + ulong[] previousValues = ulongPool.Rent(stepperLength); + Array.Clear(previousValues, 0, stepperLength); + ulong[] remainders = ulongPool.Rent(stepperLength); + Array.Clear(remainders, 0, stepperLength); + bool[] hasState = boolPool.Rent(stepperLength); + Array.Clear(hasState, 0, stepperLength); + + s_smallPrimeCycleSteppers = steppers; + s_smallPrimeStepperPreviousValues = previousValues; + s_smallPrimeStepperRemainders = remainders; + s_smallPrimeStepperHasState = hasState; + return steppers; + } + + int stepperCapacity = steppers.Length; + if (s_smallPrimeStepperPreviousValues is null || s_smallPrimeStepperPreviousValues.Length != stepperCapacity) + { + ulong[] newPrevious = ulongPool.Rent(stepperCapacity); + if (s_smallPrimeStepperPreviousValues is { } existingPrevious) + { + int copyLength = Math.Min(existingPrevious.Length, stepperCapacity); + System.Array.Copy(existingPrevious, newPrevious, copyLength); + if (copyLength < stepperCapacity) + { + Array.Clear(newPrevious, copyLength, stepperCapacity - copyLength); + } + ulongPool.Return(existingPrevious, clearArray: false); + } + else + { + Array.Clear(newPrevious, 0, stepperCapacity); + } + + s_smallPrimeStepperPreviousValues = newPrevious; + } + + if (s_smallPrimeStepperRemainders is null || s_smallPrimeStepperRemainders.Length != stepperCapacity) + { + ulong[] newRemainders = ulongPool.Rent(stepperCapacity); + if (s_smallPrimeStepperRemainders is { } existingRemainders) + { + int copyLength = Math.Min(existingRemainders.Length, stepperCapacity); + System.Array.Copy(existingRemainders, newRemainders, copyLength); + if (copyLength < stepperCapacity) + { + Array.Clear(newRemainders, copyLength, stepperCapacity - copyLength); + } + ulongPool.Return(existingRemainders, clearArray: false); + } + else + { + Array.Clear(newRemainders, 0, stepperCapacity); + } + + s_smallPrimeStepperRemainders = newRemainders; + } + + if (s_smallPrimeStepperHasState is null || s_smallPrimeStepperHasState.Length != stepperCapacity) + { + bool[] newHasState = boolPool.Rent(stepperCapacity); + if (s_smallPrimeStepperHasState is { } existingHasState) + { + int copyLength = Math.Min(existingHasState.Length, stepperCapacity); + System.Array.Copy(existingHasState, newHasState, copyLength); + if (copyLength < stepperCapacity) + { + Array.Clear(newHasState, copyLength, stepperCapacity - copyLength); + } + boolPool.Return(existingHasState, clearArray: false); + } + else + { + Array.Clear(newHasState, 0, stepperCapacity); + } + + s_smallPrimeStepperHasState = newHasState; + } + + return steppers; + } + + + private static bool TryPopulateSmallPrimeFactorsCpu( + ulong value, + uint limit, + Span primeTargets, + Span exponentTargets, + out int factorCount, + out ulong remaining) + { + factorCount = 0; + remaining = value; + + // value will never <= 1 in production code + // if (value <= 1UL) + // { + // return true; + // } + + int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); + // capacity will never equal 0 in production code + // if (capacity == 0) + // { + // return false; + // } + + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + ulong remainingLocal = value; + uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; + CycleRemainderStepper[] cycleSteppers = EnsureSmallPrimeCycleSteppers(primes); + ulong[] stepperPreviousValues = s_smallPrimeStepperPreviousValues!; + ulong[] stepperRemainders = s_smallPrimeStepperRemainders!; + bool[] stepperHasState = s_smallPrimeStepperHasState!; + ulong valueForSteppers = value; + + for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > effectiveLimit) + { + break; + } + + ulong primeSquare = squares[i]; + // primeSquare will never == 0 in production code + // if (primeSquare != 0 && primeSquare > remainingLocal) + if (primeSquare > remainingLocal) + { + break; + } + + // primeCandidate will never equal 0 in production code + // if (primeCandidate == 0UL) + // { + // continue; + // } + + ref CycleRemainderStepper stepper = ref cycleSteppers[i]; + ulong remainder; + + if (!stepperHasState[i]) + { + stepper.Reset(); + remainder = stepper.Initialize(valueForSteppers); + stepperHasState[i] = true; + } + else + { + ulong previousValue = stepperPreviousValues[i]; + if (valueForSteppers > previousValue) + { + remainder = stepper.ComputeNext(valueForSteppers); + } + else if (valueForSteppers < previousValue) + { + stepper.Reset(); + remainder = stepper.Initialize(valueForSteppers); + } + else + { + remainder = stepperRemainders[i]; + } + } + + stepperPreviousValues[i] = valueForSteppers; + stepperRemainders[i] = remainder; + + if (remainder != 0UL) + { + continue; + } + + ulong primeValue = primeCandidate; + int exponent = ExtractSmallPrimeExponent(ref remainingLocal, primeValue); + if (exponent == 0) + { + continue; + } + + if (factorCount >= capacity) + { + throw new InvalidOperationException($"Capacity is smaller than factor count"); + } + + primeTargets[factorCount] = primeValue; + exponentTargets[factorCount] = exponent; + factorCount++; + } + + remaining = remainingLocal; + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + // Unroll the initial divisions because small exponents dominate in practice. + private static int ExtractSmallPrimeExponent(ref ulong value, ulong primeValue) + { + ulong dividend = value; + ulong quotient = dividend / primeValue; + ulong remainder = dividend - (quotient * primeValue); + if (remainder != 0UL) + { + return 0; + } + + dividend = quotient; + int exponent = 1; + + quotient = dividend / primeValue; + remainder = dividend - (quotient * primeValue); + if (remainder != 0UL) + { + value = dividend; + return exponent; + } + + dividend = quotient; + exponent++; + + quotient = dividend / primeValue; + remainder = dividend - (quotient * primeValue); + if (remainder != 0UL) + { + value = dividend; + return exponent; + } + + dividend = quotient; + exponent++; + + quotient = dividend / primeValue; + remainder = dividend - (quotient * primeValue); + if (remainder != 0UL) + { + value = dividend; + return exponent; + } + + dividend = quotient; + exponent++; + + while (true) + { + quotient = dividend / primeValue; + remainder = dividend - (quotient * primeValue); + if (remainder != 0UL) + { + break; + } + + dividend = quotient; + exponent++; + } + + value = dividend; + return exponent; + } + + private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary counts) + { + ulong remaining = value; + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + + for (int i = 0; i < primeCount; i++) + { + uint primeCandidate = primes[i]; + if (primeCandidate > limit) + { + break; + } + + if (squares[i] > remaining) + { + break; + } + + ulong primeValue = primeCandidate; + if ((remaining % primeValue) != 0UL) + { + continue; + } + + int exponent = 0; + do + { + remaining /= primeValue; + exponent++; + } + while ((remaining % primeValue) == 0UL); + + counts[primeValue] = exponent; + } + + return remaining; + } + + private static long CreateDeadlineTimestamp(int milliseconds) + { + if (milliseconds <= 0) + { + return long.MaxValue; + } + + long stopwatchTicks = ConvertMillisecondsToStopwatchTicks(milliseconds); + long startTimestamp = Stopwatch.GetTimestamp(); + long deadline = startTimestamp + stopwatchTicks; + if (deadline < startTimestamp) + { + return long.MaxValue; + } + + return deadline; + } + + private static long ConvertMillisecondsToStopwatchTicks(int milliseconds) + { + long frequency = Stopwatch.Frequency; + long baseTicks = (frequency / 1000L) * milliseconds; + long remainder = (frequency % 1000L) * milliseconds; + return baseTicks + (remainder / 1000L); + } + + private static bool TryPollardRho(ulong n, long deadlineTimestamp, out ulong factor) + { + factor = 0UL; + if ((n & 1UL) == 0UL) + { + factor = 2UL; + return true; + } + + long timestamp = 0L; // reused for deadline checks. + while (true) + { + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + return false; + } + + ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + timestamp = Stopwatch.GetTimestamp(); + if (timestamp > deadlineTimestamp) + { + return false; + } + + x = AdvancePolynomial(x, c, n); + y = AdvancePolynomial(y, c, n); + y = AdvancePolynomial(y, c, n); + ulong diff = x > y ? x - y : y - x; + d = BinaryGcd(diff, n); + } + + if (d != n) + { + factor = d; + return true; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) + { + UInt128 value = (UInt128)x * x + c; + return (ulong)(value % modulus); + } + + private static ulong BinaryGcd(ulong a, ulong b) + { + if (a == 0UL) + { + return b; + } + + if (b == 0UL) + { + return a; + } + + int shift = BitOperations.TrailingZeroCount(a | b); + a >>= BitOperations.TrailingZeroCount(a); + + while (true) + { + b >>= BitOperations.TrailingZeroCount(b); + if (a > b) + { + (a, b) = (b, a); + } + + b -= a; + if (b == 0UL) + { + return a << shift; + } + } + } + + private static void AddFactor(Dictionary counts, ulong prime, int exponent) + { + if (counts.TryGetValue(prime, out int existing)) + { + counts[prime] = existing + exponent; + } + else + { + counts[prime] = exponent; + } + } + + private static bool TryAppendFactor(Span primes, Span exponents, ref int count, ulong prime, int exponent) + { + if (exponent <= 0) + { + return true; + } + + for (int i = 0; i < count; i++) + { + if (primes[i] == prime) + { + exponents[i] += exponent; + return true; + } + } + + if (count >= primes.Length) + { + return false; + } + + primes[count] = prime; + exponents[count] = exponent; + count++; + return true; + } + + private static void AddFactorToCollector( + ref bool useDictionary, + ref Dictionary? counts, + Span primes, + Span exponents, + ref int count, + ulong prime, + int exponent) + { + if (exponent <= 0) + { + return; + } + + if (useDictionary) + { + AddFactor(counts!, prime, exponent); + return; + } + + if (TryAppendFactor(primes, exponents, ref count, prime, exponent)) + { + return; + } + + counts ??= ThreadStaticPools.RentUlongIntDictionary(Math.Max(count, 8)); + CopyFactorsToDictionary(primes, exponents, count, counts); + count = 0; + useDictionary = true; + AddFactor(counts, prime, exponent); + } + + private static void CopyFactorsToDictionary( + in ReadOnlySpan primes, + in ReadOnlySpan exponents, + int count, + Dictionary target) + { + for (int i = 0; i < count; i++) + { + ulong prime = primes[i]; + int exponent = exponents[i]; + if (prime == 0UL || exponent == 0) + { + continue; + } + + target[prime] = exponent; + } + } + + private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData) + { + ulong phi = prime - 1UL; + Dictionary counts = new(capacity: 8); + FactorCompletely(phi, counts); + if (counts.Count == 0) + { + return phi; + } + + List> entries = [.. counts]; + entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); + + ulong order = phi; + int entryCount = entries.Count; + for (int i = 0; i < entryCount; i++) + { + ulong primeFactor = entries[i].Key; + int exponent = entries[i].Value; + for (int iteration = 0; iteration < exponent; iteration++) + { + if ((order % primeFactor) != 0UL) + { + break; + } + + ulong candidate = order / primeFactor; + if (Pow2EqualsOneCpu(candidate, prime, divisorData)) + { + order = candidate; + continue; + } + + break; + } + } + + return order; + } + + private static void FactorCompletely(ulong value, Dictionary counts) + { + if (value <= 1UL) + { + return; + } + + if (Open.Numeric.Primes.Prime.Numbers.IsPrime(value)) + { + AddFactor(counts, value, 1); + return; + } + + ulong factor = PollardRhoStrict(value); + FactorCompletely(factor, counts); + FactorCompletely(value / factor, counts); + } + + private static ulong PollardRhoStrict(ulong n) + { + if ((n & 1UL) == 0UL) + { + return 2UL; + } + + while (true) + { + ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; + ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + x = AdvancePolynomial(x, c, n); + y = AdvancePolynomial(y, c, n); + y = AdvancePolynomial(y, c, n); + ulong diff = x > y ? x - y : y - x; + d = BinaryGcd(diff, n); + } + + if (d != n) + { + return d; + } + } + } + + private readonly struct CandidateKey(int primary, long secondary, long tertiary) + { + public readonly int Primary = primary; + public readonly long Secondary = secondary; + public readonly long Tertiary = tertiary; + } +} +} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs new file mode 100644 index 0000000..42b977a --- /dev/null +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Gpu.cs @@ -0,0 +1,48 @@ +using System; +using System.Buffers; +using System.Collections.Generic; +using PerfectNumbers.Core.Gpu; + +namespace PerfectNumbers.Core; + +internal static partial class PrimeOrderCalculator +{ + private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; + + private const int GpuSmallPrimeFactorSlots = 64; + + private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) + { + ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); + int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); + Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); + primeBuffer.Clear(); + exponentBuffer.Clear(); + + if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) + { + remaining = value; + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); + return false; + } + + remaining = gpuRemaining; + for (int i = 0; i < factorCount; i++) + { + ulong primeValue = primeBuffer[i]; + int exponent = exponentBuffer[i]; + if (primeValue == 0UL || exponent == 0) + { + continue; + } + + counts[primeValue] = exponent; + } + + ThreadStaticPools.UlongPool.Return(primeBufferArray); + ThreadStaticPools.IntPool.Return(exponentBufferArray); + return true; + } +} diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs similarity index 92% rename from PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs rename to PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs index 4850853..5459159 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Wide.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Shared.cs @@ -1,54 +1,90 @@ +using System; using System.Buffers; using System.Collections.Generic; using System.Diagnostics; using System.Numerics; +using System.Runtime.CompilerServices; using System.Threading; -using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - public static UInt128 Calculate( - in UInt128 prime, - in UInt128? previousOrder, - in PrimeOrderSearchConfig config, - in PrimeOrderHeuristicDevice device) - { - var scope = UsePow2Mode(device); - MontgomeryDivisorData divisorData; - UInt128 result; - if (prime <= ulong.MaxValue) - { - ulong? previous = null; - if (previousOrder.HasValue) - { - UInt128 previousValue = previousOrder.Value; - if (previousValue <= ulong.MaxValue) - { - previous = (ulong)previousValue; - } - else - { - previous = ulong.MaxValue; - } - } + [ThreadStatic] + private static bool s_pow2ModeInitialized; - ulong prime64 = (ulong)prime; - divisorData = MontgomeryDivisorData.FromModulus(prime64); - ulong order64 = Calculate(prime64, previous, divisorData, config, device); - result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; - } - else + [ThreadStatic] + private static bool s_allowGpuPow2; + + [ThreadStatic] + private static PrimeOrderHeuristicDevice s_deviceMode; + + [ThreadStatic] + private static bool s_debugLoggingEnabled; + + private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; + + private readonly struct Pow2ModeScope { - divisorData = default; - result = CalculateWideInternal(prime, previousOrder, divisorData, config); + private readonly bool _previousInitialized; + private readonly bool _previousAllow; + private readonly PrimeOrderHeuristicDevice _previousDevice; + + public Pow2ModeScope(PrimeOrderHeuristicDevice device) + { + _previousInitialized = s_pow2ModeInitialized; + _previousAllow = s_allowGpuPow2; + _previousDevice = s_deviceMode; + s_pow2ModeInitialized = true; + s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; + s_deviceMode = device; + } + + public void Dispose() + { + s_allowGpuPow2 = _previousAllow; + s_pow2ModeInitialized = _previousInitialized; + s_deviceMode = _previousDevice; + } } - scope.Dispose(); - return result; - } + private readonly struct DebugLoggingScope + { + private readonly bool _previous; + public DebugLoggingScope(bool enabled) + { + _previous = s_debugLoggingEnabled; + s_debugLoggingEnabled = enabled; + } + + public void Dispose() + { + s_debugLoggingEnabled = _previous; + } + } + + private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device); + + private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled); + + [Conditional("DEBUG")] + private static void DebugLog(string message) + { + if (s_debugLoggingEnabled) + { + Console.WriteLine(message); + } + } + + [Conditional("DEBUG")] + private static void DebugLog(Func messageFactory) + { + if (s_debugLoggingEnabled) + { + Console.WriteLine(messageFactory()); + } + } private static UInt128 CalculateWideInternal(in UInt128 prime, in UInt128? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { if (prime <= UInt128.One) diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index a823bdf..fcd17a9 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -1,1614 +1,83 @@ -using System.Buffers; -using System.Collections.Generic; -using System.Buffers.Binary; -using System.Diagnostics; using System.Numerics; -using System.Runtime.CompilerServices; -using System.Runtime.InteropServices; -using System.Threading; using PerfectNumbers.Core.Gpu; namespace PerfectNumbers.Core; internal static partial class PrimeOrderCalculator { - internal enum PrimeOrderMode - { - Heuristic, - Strict, - } - - internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) - { - public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 10_240, maxPowChecks: 24, PrimeOrderMode.Heuristic); - public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); - - public readonly uint SmallFactorLimit = smallFactorLimit; - public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; - public readonly int MaxPowChecks = maxPowChecks; - public readonly PrimeOrderMode Mode = mode; - } - - internal enum PrimeOrderHeuristicDevice - { - Cpu, - Gpu, - } - - [ThreadStatic] - private static bool s_pow2ModeInitialized; - - [ThreadStatic] - private static bool s_allowGpuPow2; - - private static bool IsGpuPow2Allowed => s_pow2ModeInitialized && s_allowGpuPow2; - - [ThreadStatic] - private static PrimeOrderHeuristicDevice s_deviceMode; - - [ThreadStatic] - private static bool s_debugLoggingEnabled; - [ThreadStatic] - private static CycleRemainderStepper[]? s_smallPrimeCycleSteppers; - - [ThreadStatic] - private static ulong[]? s_smallPrimeStepperPreviousValues; - - [ThreadStatic] - private static ulong[]? s_smallPrimeStepperRemainders; - - [ThreadStatic] - private static bool[]? s_smallPrimeStepperHasState; - - - - - - - private static bool IsGpuHeuristicDevice => s_pow2ModeInitialized && s_deviceMode == PrimeOrderHeuristicDevice.Gpu; - - private readonly struct Pow2ModeScope - { - private readonly bool _previousInitialized; - private readonly bool _previousAllow; - private readonly PrimeOrderHeuristicDevice _previousDevice; - - public Pow2ModeScope(PrimeOrderHeuristicDevice device) - { - _previousInitialized = s_pow2ModeInitialized; - _previousAllow = s_allowGpuPow2; - _previousDevice = s_deviceMode; - s_pow2ModeInitialized = true; - s_allowGpuPow2 = device == PrimeOrderHeuristicDevice.Gpu; - s_deviceMode = device; - } - - public void Dispose() - { - s_allowGpuPow2 = _previousAllow; - s_pow2ModeInitialized = _previousInitialized; - s_deviceMode = _previousDevice; - } - } - - private readonly struct DebugLoggingScope - { - private readonly bool _previous; - - public DebugLoggingScope(bool enabled) - { - _previous = s_debugLoggingEnabled; - s_debugLoggingEnabled = enabled; - } - - public void Dispose() - { - s_debugLoggingEnabled = _previous; - } - } - - private static Pow2ModeScope UsePow2Mode(PrimeOrderHeuristicDevice device) => new(device); - - private static DebugLoggingScope UseDebugLogging(bool enabled) => new(enabled); - - [Conditional("DEBUG")] - private static void DebugLog(string message) - { - if (s_debugLoggingEnabled) - { - Console.WriteLine(message); - } - } - - [Conditional("DEBUG")] - private static void DebugLog(Func messageFactory) - { - if (s_debugLoggingEnabled) - { - Console.WriteLine(messageFactory()); - } - } - - public static ulong Calculate( - ulong prime, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderSearchConfig config, - PrimeOrderHeuristicDevice device) - { - var scope = UsePow2Mode(device); - ulong order = CalculateInternal(prime, previousOrder, divisorData, config); - scope.Dispose(); - return order; - } - - private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) - { - // TODO: Is this condition ever met on EvenPerfectBitScanner's execution path? If not, we can add a clarification comment and comment out the entire block. We want to support p candidates at least greater or equal to 31. - if (prime <= 3UL) - { - return prime == 3UL ? 2UL : 1UL; - } - - ulong phi = prime - 1UL; - - if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) - { - return gpuOrder; - } - - PartialFactorResult phiFactors = PartialFactor(phi, config); - - ulong result; - if (phiFactors.Factors is null) - { - result = FinishStrictlyCpu(prime, divisorData); - phiFactors.Dispose(); - return result; - } - - result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); - phiFactors.Dispose(); - return result; - } - - private static ulong RunHeuristicPipelineCpu( - ulong prime, - ulong? previousOrder, - in PrimeOrderSearchConfig config, - in MontgomeryDivisorData divisorData, - ulong phi, - PartialFactorResult phiFactors) - { - if (phiFactors.FullyFactored && TrySpecialMaxCpu(phi, prime, phiFactors, divisorData)) - { - return phi; - } - - ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); - candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); - - if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config)) - { - return candidateOrder; - } - - if (config.Mode == PrimeOrderMode.Strict) - { - return FinishStrictlyCpu(prime, divisorData); - } - - if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) - { - return order; - } - - return candidateOrder; - } - - private static ulong FinishStrictlyCpu(ulong prime, in MontgomeryDivisorData divisorData) - { - return CalculateByFactorizationCpu(prime, divisorData); - } - - private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - for (int i = 0; i < length; i++) - { - ulong factor = factorSpan[i].Value; - ulong reduced = phi / factor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } - } - - return true; - } - - private static ulong InitializeStartingOrderCpu(ulong prime, ulong phi, in MontgomeryDivisorData divisorData) - { - ulong order = phi; - if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) - { - ulong half = phi >> 1; - if (Pow2EqualsOneCpu(half, prime, divisorData)) - { - order = half; - } - } - - return order; - } - - private static ulong ExponentLoweringCpu(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) - { - ArrayPool pool = ThreadStaticPools.FactorEntryPool; - - ReadOnlySpan factorSpan = factors.Factors; - int length = factors.Count; - FactorEntry[] tempArray = pool.Rent(length + 1); - - Span buffer = tempArray; - factorSpan.CopyTo(buffer); - - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factors.Cofactor); - // bool isPrime = PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None); - if (!factors.FullyFactored && factors.Cofactor > 1UL && isPrime) - { - buffer[length] = new FactorEntry(factors.Cofactor, 1); - length++; - } - - buffer.Slice(0, length).Sort(static (a, b) => a.Value.CompareTo(b.Value)); - - for (int i = 0; i < length; i++) - { - ulong primeFactor = buffer[i].Value; - int exponent = buffer[i].Exponent; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != 0UL) - { - break; - } - - ulong reduced = order / primeFactor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - order = reduced; - continue; - } - - break; - } - } - - pool.Return(tempArray, clearArray: false); - return order; - } - - private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) - { - if (order == 0UL) - { - return false; - } - - // Calculating `a^order ≡ 1 (mod p)` is a prerequisite for `order` being the actual order of 2 modulo `p`. - // DebugLog("Verifying a^order ≡ 1 (mod p)"); - if (!Pow2EqualsOneCpu(order, prime, divisorData)) - { - return false; - } - - // DebugLog("Partial factoring order"); - - // TODO: Do we do partial factoring of order multiple times? - PartialFactorResult factorization = PartialFactor(order, config); - try - { - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL) - { - return false; - } - - // DebugLog("Cofactor > 1. Testing primality of cofactor"); - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); - // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); - if (!isPrime) - { - return false; - } - - // DebugLog("Adding cofactor as prime factor"); - PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); - factorization.Dispose(); - factorization = extended; - } - - ReadOnlySpan span = factorization.Factors; - // DebugLog("Verifying prime-power reductions"); - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - ulong primeFactor = span[i].Value; - ulong reduced = order; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } - } - } - - return true; - } - finally - { - factorization.Dispose(); - } - } - - private static bool TryHeuristicFinishCpu( - ulong prime, - ulong order, - ulong? previousOrder, - in MontgomeryDivisorData divisorData, - in PrimeOrderSearchConfig config, - PartialFactorResult phiFactors, - out ulong result) - { - result = 0UL; - if (order <= 1UL) - { - return false; - } - - // TODO: Do we do partial factoring of order multiple times? - // DebugLog("Trying heuristic. Partial factoring order"); - PartialFactorResult orderFactors = PartialFactor(order, config); - try - { - if (orderFactors.Factors is null) - { - return false; - } - - if (!orderFactors.FullyFactored) - { - if (orderFactors.Cofactor <= 1UL) - { - return false; - } - - var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor); - // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None)) - if (!isPrime) - { - return false; - } - - PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); - orderFactors.Dispose(); - orderFactors = extended; - } - - int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; - List candidates = ThreadStaticPools.RentUlongList(capacity); - FactorEntry[] factorArray = orderFactors.Factors!; - // DebugLog("Building candidates list"); - BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); - if (candidates.Count == 0) - { - ThreadStaticPools.ReturnUlongList(candidates); - return false; - } - - // DebugLog("Sorting candidates"); - SortCandidates(prime, previousOrder, candidates); - - int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; - int powUsed = 0; - int candidateCount = candidates.Count; - bool allowGpuBatch = true; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); - - // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); - int index = 0; - const int MaxGpuBatchSize = 256; - const int StackGpuBatchSize = 64; - Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize]; - ArrayPool pool = ThreadStaticPools.UlongPool; - while (index < candidateCount && powUsed < powBudget) - { - int remaining = candidateCount - index; - int budgetRemaining = powBudget - powUsed; - int batchSize = Math.Min(remaining, Math.Min(budgetRemaining, MaxGpuBatchSize)); - if (batchSize <= 0) - { - break; - } - - ReadOnlySpan batch = candidateSpan.Slice(index, batchSize); - ulong[]? gpuPool = null; - Span pooledGpuRemainders = default; - bool gpuSuccess = false; - bool gpuStackRemainders = false; - GpuPow2ModStatus status = GpuPow2ModStatus.Unavailable; - - if (allowGpuBatch && IsGpuPow2Allowed) - { - if (batchSize <= StackGpuBatchSize) - { - Span localRemainders = stackGpuRemainders.Slice(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - gpuSuccess = true; - gpuStackRemainders = true; - } - } - else - { - gpuPool = pool.Rent(batchSize); - Span pooledRemainders = gpuPool.AsSpan(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - pooledGpuRemainders = pooledRemainders; - gpuSuccess = true; - } - else - { - pool.Return(gpuPool, clearArray: false); - gpuPool = null; - } - } - - if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) - { - allowGpuBatch = false; - } - } - - for (int i = 0; i < batchSize && powUsed < powBudget; i++) - { - ulong candidate = batch[i]; - powUsed++; - - bool equalsOne; - if (gpuSuccess) - { - ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; - equalsOne = remainderValue == 1UL; - } - else - { - equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData); - } - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - ThreadStaticPools.ReturnUlongList(candidates); - result = candidate; - return true; - } - - if (gpuPool is not null) - { - pool.Return(gpuPool, clearArray: false); - } - - index += batchSize; - } - - ThreadStaticPools.ReturnUlongList(candidates); - // DebugLog("No candidate confirmed"); - return false; - } - finally - { - orderFactors.Dispose(); - } - } - - private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) - { - ulong previous = previousOrder ?? 0UL; - int previousGroup = previousOrder.HasValue ? GetGroup(previousOrder.Value, prime) : 1; - bool hasPrevious = previousOrder.HasValue; - - candidates.Sort((x, y) => - { - CandidateKey keyX = BuildKey(x, prime, previous, previousGroup, hasPrevious); - CandidateKey keyY = BuildKey(y, prime, previous, previousGroup, hasPrevious); - int compare = keyX.Primary.CompareTo(keyY.Primary); - if (compare != 0) - { - return compare; - } - - compare = keyX.Secondary.CompareTo(keyY.Secondary); - if (compare != 0) - { - return compare; - } - - return keyX.Tertiary.CompareTo(keyY.Tertiary); - }); - } - - private static CandidateKey BuildKey(ulong value, ulong prime, ulong previous, int previousGroup, bool hasPrevious) - { - int group = GetGroup(value, prime); - if (group == 0) - { - return new CandidateKey(int.MaxValue, long.MaxValue, long.MaxValue); - } - - bool isGe = !hasPrevious || value >= previous; - int primary = ComputePrimary(group, isGe, previousGroup); - long secondary; - long tertiary; - - if (group == 3) - { - secondary = -(long)value; - tertiary = -(long)value; - } - else - { - ulong reference = hasPrevious ? previous : 0UL; - ulong distance = hasPrevious ? (value > reference ? value - reference : reference - value) : value; - secondary = (long)distance; - tertiary = (long)value; - } - - return new CandidateKey(primary, secondary, tertiary); - } - - private static int ComputePrimary(int group, bool isGe, int previousGroup) - { - int groupOffset = group switch - { - 1 => 0, - 2 => 2, - 3 => 4, - _ => 6, - }; - - if (group == previousGroup) - { - if (group == 3) - { - return groupOffset + (isGe ? 0 : 3); - } - - return groupOffset + (isGe ? 0 : 1); - } - - return groupOffset + (isGe ? 0 : 1); - } - - private static int GetGroup(ulong value, ulong prime) - { - ulong threshold1 = prime >> 3; - if (value <= threshold1) - { - return 1; - } - - ulong threshold2 = prime >> 2; - if (value <= threshold2) - { - return 2; - } - - ulong threshold3 = (prime * 3UL) >> 3; - if (value <= threshold3) - { - return 3; - } - - return 0; - } - - private static void BuildCandidates(ulong order, FactorEntry[] factors, int count, List candidates, int limit) - { - if (count == 0) - { - return; - } - - Span buffer = factors.AsSpan(0, count); - buffer.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - BuildCandidatesRecursive(order, buffer, 0, 1UL, candidates, limit); - } - - private static void BuildCandidatesRecursive(ulong order, in ReadOnlySpan factors, int index, ulong divisorProduct, List candidates, int limit) - { - if (candidates.Count >= limit) - { - return; - } - - if (index >= factors.Length) - { - if (divisorProduct == 1UL || divisorProduct == order) - { - return; - } - - ulong candidate = order / divisorProduct; - if (candidate > 1UL && candidate < order) - { - candidates.Add(candidate); - } - - return; - } - - FactorEntry factor = factors[index]; - ulong primeFactor = factor.Value; - ulong contribution = 1UL; - for (int exponent = 0; exponent <= factor.Exponent; exponent++) - { - ulong nextDivisor = divisorProduct * contribution; - if (nextDivisor > order) - { - break; - } - - BuildCandidatesRecursive(order, factors, index + 1, nextDivisor, candidates, limit); - if (candidates.Count >= limit) - { - return; - } - - if (exponent == factor.Exponent) - { - break; - } - - if (contribution > order / primeFactor) - { - break; - } - - contribution *= primeFactor; - } - } - - private static bool TryConfirmCandidateCpu(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, ref int powUsed, int powBudget) - { - PartialFactorResult factorization = PartialFactor(candidate, config); - try - { - if (factorization.Factors is null) - { - return false; - } - - if (!factorization.FullyFactored) - { - if (factorization.Cofactor <= 1UL) - { - return false; - } - - bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); - // bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); - if (!isPrime) - { - return false; - } - - PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); - factorization.Dispose(); - factorization = extended; - } - - ReadOnlySpan span = factorization.Factors; - int length = factorization.Count; - for (int i = 0; i < length; i++) - { - ulong primeFactor = span[i].Value; - ulong reduced = candidate; - for (int iteration = 0; iteration < span[i].Exponent; iteration++) - { - if ((reduced % primeFactor) != 0UL) - { - break; - } - - reduced /= primeFactor; - if (powUsed >= powBudget && powBudget > 0) - { - return false; - } - - powUsed++; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } - } - } - - return true; - } - finally - { - factorization.Dispose(); - } - } - - private static bool Pow2EqualsOneCpu(ulong exponent, ulong prime, in MontgomeryDivisorData divisorData) - { - if (IsGpuPow2Allowed) - { - GpuPow2ModStatus status = PrimeOrderGpuHeuristics.TryPow2Mod(exponent, prime, out ulong remainder, divisorData); - if (status == GpuPow2ModStatus.Success) - { - return remainder == 1UL; - } - } - - return exponent.Pow2MontgomeryModWindowedCpu(divisorData, keepMontgomery: false) == 1UL; - } - - [ThreadStatic] - private static Dictionary? _primalityCache; - - private static PartialFactorResult PartialFactor(ulong value, in PrimeOrderSearchConfig config) - { - if (value <= 1UL) - { - return PartialFactorResult.Empty; - } - - const int FactorSlotCount = GpuSmallPrimeFactorSlots; - - // stackalloc is faster than pooling - // ulong[] primeSlotsArray = ThreadStaticPools.UlongPool.Rent(FactorSlotCount); - // int[] exponentSlotsArray = ThreadStaticPools.IntPool.Rent(FactorSlotCount); - // Span primeSlots = primeSlotsArray.AsSpan(0, FactorSlotCount); - // Span exponentSlots = exponentSlotsArray.AsSpan(0, FactorSlotCount); - Span primeSlots = stackalloc ulong[FactorSlotCount]; - Span exponentSlots = stackalloc int[FactorSlotCount]; - // We don't need to worry about leftovers, because we always use indexes within the calculated counts - // primeSlots.Clear(); - // exponentSlots.Clear(); - - int factorCount = 0; - Dictionary? counts = null; - bool useDictionary = false; - - uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; - ulong remaining = value; - - bool gpuFactored = false; - if (IsGpuHeuristicDevice) - { - gpuFactored = PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeSlots, exponentSlots, out factorCount, out remaining, out _); - } - - List pending = ThreadStaticPools.RentUlongList(2); - Stack? compositeStack = null; - PartialFactorResult result; - var cache = _primalityCache ??= []; - - if (!gpuFactored) - { - // We don't need to worry about leftovers, because we always use indexes within the calculated counts - // primeSlots.Clear(); - // exponentSlots.Clear(); - if (!TryPopulateSmallPrimeFactorsCpu(value, limit, primeSlots, exponentSlots, out factorCount, out remaining)) - { - useDictionary = true; - counts = ThreadStaticPools.RentUlongIntDictionary(Math.Max(factorCount, 8)); - remaining = PopulateSmallPrimeFactorsCpu(value, limit, counts); - } - } - - if (remaining > 1UL) - { - pending.Add(remaining); - } - - if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) - { - long deadlineTimestamp = CreateDeadlineTimestamp(config.PollardRhoMilliseconds); - compositeStack = ThreadStaticPools.RentUlongStack(Math.Max(pending.Count << 1, 4)); - compositeStack.Push(remaining); - pending.Clear(); - - bool pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; - if (!pollardRhoDeadlineReached) - { - while (compositeStack.Count > 0) - { - ulong composite = compositeStack.Pop(); - // composite will never be smaller than 1 on the execution path - // if (composite <= 1UL) - // { - // continue; - // } - - if (!cache.TryGetValue(composite, out bool isPrime)) - { - isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - cache.Add(composite, isPrime); - } - - if (isPrime) - { - AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); - continue; - } - - if (!TryPollardRho(composite, deadlineTimestamp, out ulong factor)) - { - pending.Add(composite); - continue; - } - - ulong quotient = composite / factor; - compositeStack.Push(factor); - compositeStack.Push(quotient); - } - - pollardRhoDeadlineReached = Stopwatch.GetTimestamp() > deadlineTimestamp; - } - - if (pollardRhoDeadlineReached) - { - while (compositeStack.Count > 0) - { - pending.Add(compositeStack.Pop()); - } - } - } - - ulong cofactor = 1UL; - int pendingCount = pending.Count; - for (int i = 0; i < pendingCount; i++) - { - ulong composite = pending[i]; - if (!cache.TryGetValue(composite, out bool isPrime)) - { - isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(composite); - cache.Add(composite, isPrime); - } - - if (isPrime) - { - AddFactorToCollector(ref useDictionary, ref counts, primeSlots, exponentSlots, ref factorCount, composite, 1); - } - // composite is never smaller on the execution path - // else if (composite > 1UL) - else - { - cofactor = checked(cofactor * composite); - } - } - - ArrayPool pool = ThreadStaticPools.FactorEntryPool; - if (useDictionary) - { - if ((counts is null || counts.Count == 0) && cofactor == value) - { - result = PartialFactorResult.Rent(null, value, false, 0); - goto ReturnResult; - } - - if (counts is null || counts.Count == 0) - { - result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - goto ReturnResult; - } - - FactorEntry[] factors = pool.Rent(counts.Count); - int index = 0; - foreach (KeyValuePair entry in counts) - { - factors[index] = new FactorEntry(entry.Key, entry.Value); - index++; - } - - Array.Sort(factors, static (a, b) => a.Value.CompareTo(b.Value)); - - bool fullyFactored = cofactor == 1UL; - result = PartialFactorResult.Rent(factors, cofactor, fullyFactored, index); - goto ReturnResult; - } - - // This will never happen in production code. We'll always get at least 1 factor - // if (factorCount == 0) - // { - // if (cofactor == value) - // { - // result = PartialFactorResult.Rent(null, value, false, 0); - // goto ReturnResult; - // } - - // result = PartialFactorResult.Rent(null, cofactor, cofactor == 1UL, 0); - // goto ReturnResult; - // } - - FactorEntry[] array = pool.Rent(factorCount); - int arrayIndex = 0; - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeSlots[i]; - int exponentValue = 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[arrayIndex] = new FactorEntry(primeValue, exponentValue); - arrayIndex++; - } - - Span arraySpan = array.AsSpan(0, arrayIndex); - arraySpan.Sort(static (a, b) => a.Value.CompareTo(b.Value)); - bool fullyFactoredArray = cofactor == 1UL; - result = PartialFactorResult.Rent(array, cofactor, fullyFactoredArray, arrayIndex); - - ReturnResult: - // ThreadStaticPools.UlongPool.Return(primeSlotsArray); - // ThreadStaticPools.IntPool.Return(exponentSlotsArray); - - if (counts is not null) - { - ThreadStaticPools.ReturnUlongIntDictionary(counts); - } - - if (compositeStack is not null) - { - compositeStack.Clear(); - ThreadStaticPools.ReturnUlongStack(compositeStack); - } - - pending.Clear(); - ThreadStaticPools.ReturnUlongList(pending); - return result; - } - - private static CycleRemainderStepper[] EnsureSmallPrimeCycleSteppers(uint[] primes) - { - int primeCount = primes.Length; - CycleRemainderStepper[]? steppers = s_smallPrimeCycleSteppers; - if (steppers is null || steppers.Length < primeCount) - { - int length = primeCount; - steppers = new CycleRemainderStepper[length]; - for (int i = 0; i < length; i++) - { - steppers[i] = new CycleRemainderStepper(primes[i]); - } - - s_smallPrimeCycleSteppers = steppers; - s_smallPrimeStepperPreviousValues = new ulong[length]; - s_smallPrimeStepperRemainders = new ulong[length]; - s_smallPrimeStepperHasState = new bool[length]; - return steppers; - } - - if (s_smallPrimeStepperPreviousValues is null || s_smallPrimeStepperPreviousValues.Length != steppers.Length) - { - var previousValues = new ulong[steppers.Length]; - if (s_smallPrimeStepperPreviousValues is { } existingPrevious) - { - System.Array.Copy(existingPrevious, previousValues, Math.Min(existingPrevious.Length, previousValues.Length)); - } - - s_smallPrimeStepperPreviousValues = previousValues; - } - - if (s_smallPrimeStepperRemainders is null || s_smallPrimeStepperRemainders.Length != steppers.Length) - { - var remainders = new ulong[steppers.Length]; - if (s_smallPrimeStepperRemainders is { } existingRemainders) - { - System.Array.Copy(existingRemainders, remainders, Math.Min(existingRemainders.Length, remainders.Length)); - } - - s_smallPrimeStepperRemainders = remainders; - } - - if (s_smallPrimeStepperHasState is null || s_smallPrimeStepperHasState.Length != steppers.Length) - { - var hasState = new bool[steppers.Length]; - if (s_smallPrimeStepperHasState is { } existingHasState) - { - System.Array.Copy(existingHasState, hasState, Math.Min(existingHasState.Length, hasState.Length)); - } - - s_smallPrimeStepperHasState = hasState; - } - - return steppers; - } - - private const int GpuSmallPrimeFactorSlots = 64; - - private static bool TryPopulateSmallPrimeFactorsCpu( - ulong value, - uint limit, - Span primeTargets, - Span exponentTargets, - out int factorCount, - out ulong remaining) - { - factorCount = 0; - remaining = value; - - // value will never <= 1 in production code - // if (value <= 1UL) - // { - // return true; - // } - - int capacity = Math.Min(primeTargets.Length, exponentTargets.Length); - // capacity will never equal 0 in production code - // if (capacity == 0) - // { - // return false; - // } - - uint[] primes = PrimesGenerator.SmallPrimes; - ulong[] squares = PrimesGenerator.SmallPrimesPow2; - int primeCount = primes.Length; - ulong remainingLocal = value; - uint effectiveLimit = limit == 0 ? uint.MaxValue : limit; - CycleRemainderStepper[] cycleSteppers = EnsureSmallPrimeCycleSteppers(primes); - ulong[] stepperPreviousValues = s_smallPrimeStepperPreviousValues!; - ulong[] stepperRemainders = s_smallPrimeStepperRemainders!; - bool[] stepperHasState = s_smallPrimeStepperHasState!; - ulong valueForSteppers = value; - - for (int i = 0; i < primeCount && remainingLocal > 1UL; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > effectiveLimit) - { - break; - } - - ulong primeSquare = squares[i]; - // primeSquare will never == 0 in production code - // if (primeSquare != 0 && primeSquare > remainingLocal) - if (primeSquare > remainingLocal) - { - break; - } - - // primeCandidate will never equal 0 in production code - // if (primeCandidate == 0UL) - // { - // continue; - // } - - ref CycleRemainderStepper stepper = ref cycleSteppers[i]; - ulong remainder; - - if (!stepperHasState[i]) - { - stepper.Reset(); - remainder = stepper.Initialize(valueForSteppers); - stepperHasState[i] = true; - } - else - { - ulong previousValue = stepperPreviousValues[i]; - if (valueForSteppers > previousValue) - { - remainder = stepper.ComputeNext(valueForSteppers); - } - else if (valueForSteppers < previousValue) - { - stepper.Reset(); - remainder = stepper.Initialize(valueForSteppers); - } - else - { - remainder = stepperRemainders[i]; - } - } - - stepperPreviousValues[i] = valueForSteppers; - stepperRemainders[i] = remainder; - - if (remainder != 0UL) - { - continue; - } - - ulong primeValue = primeCandidate; - int exponent = ExtractSmallPrimeExponent(ref remainingLocal, primeValue); - if (exponent == 0) - { - continue; - } - - if (factorCount >= capacity) - { - throw new InvalidOperationException($"Capacity is smaller than factor count"); - } - - primeTargets[factorCount] = primeValue; - exponentTargets[factorCount] = exponent; - factorCount++; - } - - remaining = remainingLocal; - return true; - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - // Unroll the initial divisions because small exponents dominate in practice. - private static int ExtractSmallPrimeExponent(ref ulong value, ulong primeValue) - { - ulong dividend = value; - ulong quotient = dividend / primeValue; - ulong remainder = dividend - (quotient * primeValue); - if (remainder != 0UL) - { - return 0; - } - - dividend = quotient; - int exponent = 1; - - quotient = dividend / primeValue; - remainder = dividend - (quotient * primeValue); - if (remainder != 0UL) - { - value = dividend; - return exponent; - } - - dividend = quotient; - exponent++; - - quotient = dividend / primeValue; - remainder = dividend - (quotient * primeValue); - if (remainder != 0UL) - { - value = dividend; - return exponent; - } - - dividend = quotient; - exponent++; - - quotient = dividend / primeValue; - remainder = dividend - (quotient * primeValue); - if (remainder != 0UL) - { - value = dividend; - return exponent; - } - - dividend = quotient; - exponent++; - - while (true) - { - quotient = dividend / primeValue; - remainder = dividend - (quotient * primeValue); - if (remainder != 0UL) - { - break; - } - - dividend = quotient; - exponent++; - } - - value = dividend; - return exponent; - } - - private static ulong PopulateSmallPrimeFactorsCpu(ulong value, uint limit, Dictionary counts) - { - ulong remaining = value; - uint[] primes = PrimesGenerator.SmallPrimes; - ulong[] squares = PrimesGenerator.SmallPrimesPow2; - int primeCount = primes.Length; - - for (int i = 0; i < primeCount; i++) - { - uint primeCandidate = primes[i]; - if (primeCandidate > limit) - { - break; - } - - if (squares[i] > remaining) - { - break; - } - - ulong primeValue = primeCandidate; - if ((remaining % primeValue) != 0UL) - { - continue; - } - - int exponent = 0; - do - { - remaining /= primeValue; - exponent++; - } - while ((remaining % primeValue) == 0UL); - - counts[primeValue] = exponent; - } - - return remaining; - } - - private static bool TryPopulateSmallPrimeFactorsGpu(ulong value, uint limit, Dictionary counts, out ulong remaining) - { - ulong[] primeBufferArray = ThreadStaticPools.UlongPool.Rent(GpuSmallPrimeFactorSlots); - int[] exponentBufferArray = ThreadStaticPools.IntPool.Rent(GpuSmallPrimeFactorSlots); - Span primeBuffer = primeBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - Span exponentBuffer = exponentBufferArray.AsSpan(0, GpuSmallPrimeFactorSlots); - primeBuffer.Clear(); - exponentBuffer.Clear(); - - if (!PrimeOrderGpuHeuristics.TryPartialFactor(value, limit, primeBuffer, exponentBuffer, out int factorCount, out ulong gpuRemaining, out _)) - { - remaining = value; - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); - return false; - } - - remaining = gpuRemaining; - for (int i = 0; i < factorCount; i++) - { - ulong primeValue = primeBuffer[i]; - int exponent = exponentBuffer[i]; - if (primeValue == 0UL || exponent == 0) - { - continue; - } - - counts[primeValue] = exponent; - } - - ThreadStaticPools.UlongPool.Return(primeBufferArray); - ThreadStaticPools.IntPool.Return(exponentBufferArray); - return true; - } - private static long CreateDeadlineTimestamp(int milliseconds) - { - if (milliseconds <= 0) - { - return long.MaxValue; - } - - long stopwatchTicks = ConvertMillisecondsToStopwatchTicks(milliseconds); - long startTimestamp = Stopwatch.GetTimestamp(); - long deadline = startTimestamp + stopwatchTicks; - if (deadline < startTimestamp) - { - return long.MaxValue; - } - - return deadline; - } - - private static long ConvertMillisecondsToStopwatchTicks(int milliseconds) - { - long frequency = Stopwatch.Frequency; - long baseTicks = (frequency / 1000L) * milliseconds; - long remainder = (frequency % 1000L) * milliseconds; - return baseTicks + (remainder / 1000L); - } - - private static bool TryPollardRho(ulong n, long deadlineTimestamp, out ulong factor) - { - factor = 0UL; - if ((n & 1UL) == 0UL) - { - factor = 2UL; - return true; - } - - long timestamp = 0L; // reused for deadline checks. - while (true) - { - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - return false; - } - - ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; - ulong y = x; - ulong d = 1UL; - - while (d == 1UL) - { - timestamp = Stopwatch.GetTimestamp(); - if (timestamp > deadlineTimestamp) - { - return false; - } - - x = AdvancePolynomial(x, c, n); - y = AdvancePolynomial(y, c, n); - y = AdvancePolynomial(y, c, n); - ulong diff = x > y ? x - y : y - x; - d = BinaryGcd(diff, n); - } - - if (d != n) - { - factor = d; - return true; - } - } - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static ulong AdvancePolynomial(ulong x, ulong c, ulong modulus) - { - UInt128 value = (UInt128)x * x + c; - return (ulong)(value % modulus); - } - - private static ulong BinaryGcd(ulong a, ulong b) - { - if (a == 0UL) - { - return b; - } - - if (b == 0UL) - { - return a; - } - - int shift = BitOperations.TrailingZeroCount(a | b); - a >>= BitOperations.TrailingZeroCount(a); - - while (true) - { - b >>= BitOperations.TrailingZeroCount(b); - if (a > b) - { - (a, b) = (b, a); - } - - b -= a; - if (b == 0UL) - { - return a << shift; - } - } - } - - private static void AddFactor(Dictionary counts, ulong prime, int exponent) - { - if (counts.TryGetValue(prime, out int existing)) - { - counts[prime] = existing + exponent; - } - else - { - counts[prime] = exponent; - } - } - - private static bool TryAppendFactor(Span primes, Span exponents, ref int count, ulong prime, int exponent) - { - if (exponent <= 0) - { - return true; - } - - for (int i = 0; i < count; i++) - { - if (primes[i] == prime) - { - exponents[i] += exponent; - return true; - } - } - - if (count >= primes.Length) - { - return false; - } - - primes[count] = prime; - exponents[count] = exponent; - count++; - return true; - } - - private static void AddFactorToCollector( - ref bool useDictionary, - ref Dictionary? counts, - Span primes, - Span exponents, - ref int count, - ulong prime, - int exponent) - { - if (exponent <= 0) - { - return; - } - - if (useDictionary) - { - AddFactor(counts!, prime, exponent); - return; - } - - if (TryAppendFactor(primes, exponents, ref count, prime, exponent)) - { - return; - } - - counts ??= ThreadStaticPools.RentUlongIntDictionary(Math.Max(count, 8)); - CopyFactorsToDictionary(primes, exponents, count, counts); - count = 0; - useDictionary = true; - AddFactor(counts, prime, exponent); - } - - private static void CopyFactorsToDictionary( - in ReadOnlySpan primes, - in ReadOnlySpan exponents, - int count, - Dictionary target) - { - for (int i = 0; i < count; i++) - { - ulong prime = primes[i]; - int exponent = exponents[i]; - if (prime == 0UL || exponent == 0) - { - continue; - } - - target[prime] = exponent; - } - } - - private static ulong CalculateByFactorizationCpu(ulong prime, in MontgomeryDivisorData divisorData) - { - ulong phi = prime - 1UL; - Dictionary counts = new(capacity: 8); - FactorCompletely(phi, counts); - if (counts.Count == 0) - { - return phi; - } - - List> entries = [.. counts]; - entries.Sort(static (a, b) => a.Key.CompareTo(b.Key)); - - ulong order = phi; - int entryCount = entries.Count; - for (int i = 0; i < entryCount; i++) - { - ulong primeFactor = entries[i].Key; - int exponent = entries[i].Value; - for (int iteration = 0; iteration < exponent; iteration++) - { - if ((order % primeFactor) != 0UL) - { - break; - } - - ulong candidate = order / primeFactor; - if (Pow2EqualsOneCpu(candidate, prime, divisorData)) - { - order = candidate; - continue; - } - - break; - } - } - - return order; - } - - private static void FactorCompletely(ulong value, Dictionary counts) - { - if (value <= 1UL) - { - return; - } - - if (Open.Numeric.Primes.Prime.Numbers.IsPrime(value)) - { - AddFactor(counts, value, 1); - return; - } - - ulong factor = PollardRhoStrict(value); - FactorCompletely(factor, counts); - FactorCompletely(value / factor, counts); - } - - private static ulong PollardRhoStrict(ulong n) - { - if ((n & 1UL) == 0UL) - { - return 2UL; - } - - while (true) - { - ulong c = (DeterministicRandom.NextUInt64() % (n - 1UL)) + 1UL; - ulong x = (DeterministicRandom.NextUInt64() % (n - 2UL)) + 2UL; - ulong y = x; - ulong d = 1UL; - - while (d == 1UL) - { - x = AdvancePolynomial(x, c, n); - y = AdvancePolynomial(y, c, n); - y = AdvancePolynomial(y, c, n); - ulong diff = x > y ? x - y : y - x; - d = BinaryGcd(diff, n); - } - - if (d != n) - { - return d; - } - } - } - - private readonly struct CandidateKey(int primary, long secondary, long tertiary) - { - public readonly int Primary = primary; - public readonly long Secondary = secondary; - public readonly long Tertiary = tertiary; - } + internal enum PrimeOrderMode + { + Heuristic, + Strict, + } + + internal readonly struct PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) + { + public static readonly PrimeOrderSearchConfig HeuristicDefault = new(smallFactorLimit: 4_000_000, pollardRhoMilliseconds: 24 * 10_240, maxPowChecks: 24, PrimeOrderMode.Heuristic); + public static readonly PrimeOrderSearchConfig StrictDefault = new(smallFactorLimit: 1_000_000, pollardRhoMilliseconds: 0, maxPowChecks: 0, PrimeOrderMode.Strict); + + public readonly uint SmallFactorLimit = smallFactorLimit; + public readonly int PollardRhoMilliseconds = pollardRhoMilliseconds; + public readonly int MaxPowChecks = maxPowChecks; + public readonly PrimeOrderMode Mode = mode; + } + + internal enum PrimeOrderHeuristicDevice + { + Cpu, + Gpu, + } + + public static ulong Calculate( + ulong prime, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + in PrimeOrderSearchConfig config, + PrimeOrderHeuristicDevice device) + { + var scope = UsePow2Mode(device); + ulong order = CalculateInternal(prime, previousOrder, divisorData, config); + scope.Dispose(); + return order; + } + + public static UInt128 Calculate( + in UInt128 prime, + in UInt128? previousOrder, + in PrimeOrderSearchConfig config, + in PrimeOrderHeuristicDevice device) + { + var scope = UsePow2Mode(device); + MontgomeryDivisorData divisorData; + UInt128 result; + if (prime <= ulong.MaxValue) + { + ulong? previous = null; + if (previousOrder.HasValue) + { + UInt128 previousValue = previousOrder.Value; + if (previousValue <= ulong.MaxValue) + { + previous = (ulong)previousValue; + } + else + { + previous = ulong.MaxValue; + } + } + + ulong prime64 = (ulong)prime; + divisorData = MontgomeryDivisorData.FromModulus(prime64); + ulong order64 = Calculate(prime64, previous, divisorData, config, device); + result = order64 == 0UL ? UInt128.Zero : (UInt128)order64; + } + else + { + divisorData = default; + result = CalculateWideInternal(prime, previousOrder, divisorData, config); + } + + scope.Dispose(); + return result; + } } diff --git a/PerfectNumbers.Core/ThreadStaticPools.cs b/PerfectNumbers.Core/ThreadStaticPools.cs index ed00820..8a58c2e 100644 --- a/PerfectNumbers.Core/ThreadStaticPools.cs +++ b/PerfectNumbers.Core/ThreadStaticPools.cs @@ -11,6 +11,18 @@ namespace PerfectNumbers.Core [ThreadStatic] private static ArrayPool? _factorEntryPool; + [ThreadStatic] + private static ArrayPool? _cycleRemainderStepperPool; + + public static ArrayPool CycleRemainderStepperPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _cycleRemainderStepperPool ??= ArrayPool.Create(); + } + } + public static ArrayPool FactorEntryPool { [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -70,6 +82,18 @@ namespace PerfectNumbers.Core [ThreadStatic] private static ArrayPool? _ulongPool; + [ThreadStatic] + private static ArrayPool? _boolPool; + + public static ArrayPool BoolPool + { + [MethodImpl(MethodImplOptions.AggressiveInlining)] + get + { + return _boolPool ??= ArrayPool.Create(); + } + } + public static ArrayPool UlongPool { [MethodImpl(MethodImplOptions.AggressiveInlining)]