From 4ce653d19a301895e012fc9651029dc4d8f099be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Lema=C5=84czyk?= Date: Fri, 3 Oct 2025 23:01:03 +0200 Subject: [PATCH] Accelerate divisor cycle discovery with heuristic order calculator --- PerfectNumbers.Core/MersenneDivisorCycles.cs | 29 +- PerfectNumbers.Core/PrimeOrderCalculator.cs | 857 +++++++++++++++++++ 2 files changed, 883 insertions(+), 3 deletions(-) create mode 100644 PerfectNumbers.Core/PrimeOrderCalculator.cs diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index f8cbd4c..a7a21a3 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -796,6 +796,28 @@ public class MersenneDivisorCycles return 1UL; } + if (divisor <= 3UL) + { + return CalculateCycleLengthFallback(divisor); + } + + if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) + { + PrimeOrderCalculator.PrimeOrderResult orderResult = PrimeOrderCalculator.Calculate( + divisor, + previousOrder: null, + PrimeOrderCalculator.PrimeOrderSearchConfig.HeuristicDefault); + if (orderResult.Status == PrimeOrderCalculator.PrimeOrderStatus.Found && orderResult.Order != 0UL) + { + return orderResult.Order; + } + } + + return CalculateCycleLengthFallback(divisor); + } + + private static ulong CalculateCycleLengthFallback(ulong divisor) + { // Otherwise, find order of 2 mod divisor // TODO: Switch this scalar fallback to the unrolled-hex stepping sequence once // the generator is shared with CPU callers; the benchmark shows the unrolled @@ -803,15 +825,16 @@ public class MersenneDivisorCycles // TODO: Expose a GPU-first branch here so high divisors leverage the ProcessEightBitWindows // kernel measured fastest in CycleLengthGpuVsCpuBenchmarks, returning the result without // storing it in the shared cache. - ulong order = 1UL, pow = 2UL; - while (pow != 1) + ulong order = 1UL; + ulong pow = 2UL; + while (pow != 1UL) { pow <<= 1; if (pow >= divisor) pow -= divisor; order++; - } + } return order; } diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs new file mode 100644 index 0000000..24acb4b --- /dev/null +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -0,0 +1,857 @@ +using System.Buffers; +using System.Buffers.Binary; +using System.Collections.Generic; +using System.Diagnostics; +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Security.Cryptography; +using System.Threading; + +namespace PerfectNumbers.Core; + +internal static class PrimeOrderCalculator +{ + internal enum PrimeOrderStatus + { + Found, + HeuristicUnresolved, + } + + internal readonly struct PrimeOrderResult + { + public PrimeOrderResult(PrimeOrderStatus status, ulong order) + { + Status = status; + Order = order; + } + + public PrimeOrderStatus Status { get; } + + public ulong Order { get; } + } + + internal enum PrimeOrderMode + { + Heuristic, + Strict, + } + + internal readonly struct PrimeOrderSearchConfig + { + public PrimeOrderSearchConfig(uint smallFactorLimit, int pollardRhoMilliseconds, int maxPowChecks, PrimeOrderMode mode) + { + SmallFactorLimit = smallFactorLimit; + PollardRhoMilliseconds = pollardRhoMilliseconds; + MaxPowChecks = maxPowChecks; + Mode = mode; + } + + public uint SmallFactorLimit { get; } + + public int PollardRhoMilliseconds { get; } + + public int MaxPowChecks { get; } + + public PrimeOrderMode Mode { get; } + + public static PrimeOrderSearchConfig HeuristicDefault => new(100_000, 1, 24, PrimeOrderMode.Heuristic); + + public static PrimeOrderSearchConfig StrictDefault => new(1_000_000, 0, 0, PrimeOrderMode.Strict); + } + + public static PrimeOrderResult Calculate(ulong prime, ulong? previousOrder, PrimeOrderSearchConfig config) + { + if (prime <= 3UL) + { + return new PrimeOrderResult(PrimeOrderStatus.Found, prime == 3UL ? 2UL : 1UL); + } + + ulong phi = prime - 1UL; + MontgomeryDivisorData divisorData = MontgomeryDivisorDataCache.Get(prime); + PartialFactorResult phiFactors = PartialFactor(phi, config); + + if (phiFactors.Factors is null) + { + return FinishStrictly(prime, config.Mode); + } + + if (phiFactors.FullyFactored && TrySpecialMax(phi, prime, phiFactors, divisorData)) + { + return new PrimeOrderResult(PrimeOrderStatus.Found, phi); + } + + ulong candidateOrder = InitializeStartingOrder(prime, phi, divisorData); + candidateOrder = ExponentLowering(candidateOrder, prime, phiFactors, divisorData); + + if (TryConfirmOrder(prime, candidateOrder, divisorData, config)) + { + return new PrimeOrderResult(PrimeOrderStatus.Found, candidateOrder); + } + + if (config.Mode == PrimeOrderMode.Strict) + { + return FinishStrictly(prime, PrimeOrderMode.Strict); + } + + if (TryHeuristicFinish(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) + { + return new PrimeOrderResult(PrimeOrderStatus.Found, order); + } + + return new PrimeOrderResult(PrimeOrderStatus.HeuristicUnresolved, 0UL); + } + + private static PrimeOrderResult FinishStrictly(ulong prime, PrimeOrderMode mode) + { + ulong strictOrder = CalculateByDoubling(prime); + return new PrimeOrderResult(mode == PrimeOrderMode.Strict ? PrimeOrderStatus.Found : PrimeOrderStatus.HeuristicUnresolved, strictOrder); + } + + private static bool TrySpecialMax(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 (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + { + return false; + } + } + + return true; + } + + private static ulong InitializeStartingOrder(ulong prime, ulong phi, in MontgomeryDivisorData divisorData) + { + ulong order = phi; + if ((prime & 7UL) == 1UL || (prime & 7UL) == 7UL) + { + ulong half = phi >> 1; + if (half.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + { + order = half; + } + } + + return order; + } + + private static ulong ExponentLowering(ulong order, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData) + { + FactorEntry[]? tempArray = null; + try + { + ReadOnlySpan factorSpan = factors.Factors; + int length = factors.Count; + tempArray = ArrayPool.Shared.Rent(length + 1); + Span buffer = tempArray; + factorSpan.CopyTo(buffer); + + if (!factors.FullyFactored && factors.Cofactor > 1UL && PrimeTester.IsPrimeInternal(factors.Cofactor, CancellationToken.None)) + { + 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 (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + { + order = reduced; + continue; + } + + break; + } + } + + return order; + } + finally + { + if (tempArray is not null) + { + ArrayPool.Shared.Return(tempArray, clearArray: false); + } + } + } + + private static bool TryConfirmOrder(ulong prime, ulong order, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config) + { + if (order == 0UL) + { + return false; + } + + if (order.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) != 1UL) + { + return false; + } + + PartialFactorResult factorization = PartialFactor(order, config); + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL) + { + return false; + } + + if (!PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None)) + { + return false; + } + + factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + } + + ReadOnlySpan span = factorization.Factors; + 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 (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + { + return false; + } + } + } + + return true; + } + + private static bool TryHeuristicFinish( + ulong prime, + ulong order, + ulong? previousOrder, + in MontgomeryDivisorData divisorData, + PrimeOrderSearchConfig config, + PartialFactorResult phiFactors, + out ulong result) + { + result = 0UL; + if (order <= 1UL) + { + return false; + } + + PartialFactorResult orderFactors = PartialFactor(order, config); + if (orderFactors.Factors is null) + { + return false; + } + + if (!orderFactors.FullyFactored) + { + if (orderFactors.Cofactor <= 1UL) + { + return false; + } + + if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None)) + { + return false; + } + + orderFactors = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + } + + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; + List candidates = new(capacity); + FactorEntry[] factorArray = orderFactors.Factors!; + BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); + if (candidates.Count == 0) + { + return false; + } + + SortCandidates(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 (powUsed >= powBudget) + { + break; + } + + ulong candidate = candidates[i]; + powUsed++; + + if (candidate.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) != 1UL) + { + continue; + } + + if (!TryConfirmCandidate(prime, candidate, divisorData, config, ref powUsed, powBudget)) + { + continue; + } + + result = candidate; + return true; + } + + return false; + } + + 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, 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 TryConfirmCandidate(ulong prime, ulong candidate, in MontgomeryDivisorData divisorData, PrimeOrderSearchConfig config, ref int powUsed, int powBudget) + { + PartialFactorResult factorization = PartialFactor(candidate, config); + if (factorization.Factors is null) + { + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL) + { + return false; + } + + if (!PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None)) + { + return false; + } + + factorization = factorization.WithAdditionalPrime(factorization.Cofactor); + } + + 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 (reduced.Pow2MontgomeryModWindowed(divisorData, keepMontgomery: false) == 1UL) + { + return false; + } + } + } + + return true; + } + + private static PartialFactorResult PartialFactor(ulong value, PrimeOrderSearchConfig config) + { + if (value <= 1UL) + { + return PartialFactorResult.Empty; + } + + Dictionary counts = new(capacity: 8); + ulong remaining = value; + uint[] primes = PrimesGenerator.SmallPrimes; + ulong[] squares = PrimesGenerator.SmallPrimesPow2; + int primeCount = primes.Length; + uint limit = config.SmallFactorLimit == 0 ? uint.MaxValue : config.SmallFactorLimit; + + 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; + } + + List pending = new(); + if (remaining > 1UL) + { + pending.Add(remaining); + } + + if (config.PollardRhoMilliseconds > 0 && pending.Count > 0) + { + Stopwatch stopwatch = Stopwatch.StartNew(); + long budgetTicks = TimeSpan.FromMilliseconds(config.PollardRhoMilliseconds).Ticks; + Stack stack = new(); + stack.Push(remaining); + pending.Clear(); + + while (stack.Count > 0) + { + ulong composite = stack.Pop(); + if (composite == 1UL) + { + continue; + } + + if (PrimeTester.IsPrimeInternal(composite, CancellationToken.None)) + { + AddFactor(counts, composite, 1); + continue; + } + + if (stopwatch.ElapsedTicks > budgetTicks) + { + pending.Add(composite); + continue; + } + + if (!TryPollardRho(composite, stopwatch, budgetTicks, out ulong factor)) + { + pending.Add(composite); + continue; + } + + ulong quotient = composite / factor; + stack.Push(factor); + stack.Push(quotient); + } + } + + ulong cofactor = 1UL; + int pendingCount = pending.Count; + for (int i = 0; i < pendingCount; i++) + { + ulong composite = pending[i]; + if (PrimeTester.IsPrimeInternal(composite, CancellationToken.None)) + { + AddFactor(counts, composite, 1); + } + else + { + cofactor = checked(cofactor * composite); + } + } + + if (counts.Count == 0 && cofactor == value) + { + return new PartialFactorResult(null, value, false, 0); + } + + FactorEntry[] factors = ArrayPool.Shared.Rent(counts.Count); + int index = 0; + foreach (KeyValuePair entry in counts) + { + factors[index] = new FactorEntry(entry.Key, entry.Value); + index++; + } + + Span span = factors.AsSpan(0, index); + span.Sort(static (a, b) => a.Value.CompareTo(b.Value)); + + FactorEntry[] resultArray = new FactorEntry[index]; + span.CopyTo(resultArray); + ArrayPool.Shared.Return(factors, clearArray: false); + + bool fullyFactored = cofactor == 1UL; + return new PartialFactorResult(resultArray, cofactor, fullyFactored, index); + } + + private static bool TryPollardRho(ulong n, Stopwatch stopwatch, long budgetTicks, out ulong factor) + { + factor = 0UL; + if ((n & 1UL) == 0UL) + { + factor = 2UL; + return true; + } + + Span buffer = stackalloc byte[8]; + + while (true) + { + if (stopwatch.ElapsedTicks > budgetTicks) + { + return false; + } + + RandomNumberGenerator.Fill(buffer); + ulong c = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 1UL)) + 1UL; + RandomNumberGenerator.Fill(buffer); + ulong x = (BinaryPrimitives.ReadUInt64LittleEndian(buffer) % (n - 2UL)) + 2UL; + ulong y = x; + ulong d = 1UL; + + while (d == 1UL) + { + if (stopwatch.ElapsedTicks > budgetTicks) + { + 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 ulong CalculateByDoubling(ulong prime) + { + ulong order = 1UL; + ulong pow = 2UL; + + while (pow != 1UL) + { + pow <<= 1; + if (pow >= prime) + { + pow -= prime; + } + + order++; + } + + return order; + } + + private readonly struct CandidateKey + { + public CandidateKey(int primary, long secondary, long tertiary) + { + Primary = primary; + Secondary = secondary; + Tertiary = tertiary; + } + + public int Primary { get; } + + public long Secondary { get; } + + public long Tertiary { get; } + } + + private readonly struct FactorEntry + { + public FactorEntry(ulong value, int exponent) + { + Value = value; + Exponent = exponent; + } + + public ulong Value { get; } + + public int Exponent { get; } + } + + private readonly struct PartialFactorResult + { + public PartialFactorResult(FactorEntry[]? factors, ulong cofactor, bool fullyFactored, int count) + { + Factors = factors; + Cofactor = cofactor; + FullyFactored = fullyFactored; + Count = count; + } + + public FactorEntry[]? Factors { get; } + + public ulong Cofactor { get; } + + public bool FullyFactored { get; } + + public int Count { get; } + + public static PartialFactorResult Empty => new(null, 1UL, true, 0); + + public PartialFactorResult WithAdditionalPrime(ulong prime) + { + if (Factors is null) + { + FactorEntry[] local = [new FactorEntry(prime, 1)]; + return new PartialFactorResult(local, 1UL, true, 1); + } + + FactorEntry[] extended = new FactorEntry[Count + 1]; + Array.Copy(Factors, 0, extended, 0, Count); + extended[Count] = new FactorEntry(prime, 1); + Array.Sort(extended, static (a, b) => a.Value.CompareTo(b.Value)); + return new PartialFactorResult(extended, 1UL, true, Count + 1); + } + } +}