diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index fa72a47..1f95404 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -204,9 +204,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv prime, divisorData, factorCache, - out ulong computedCycle) || computedCycle == 0UL) + out ulong computedCycle, + out bool primeOrderFailed) || computedCycle == 0UL) { - divisorCycle = cycleCache.GetCycleLength(candidate); + divisorCycle = cycleCache.GetCycleLength(candidate, skipPrimeOrderHeuristic: primeOrderFailed); } else { diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index 421154b..ba320f7 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -45,12 +45,38 @@ public sealed class DivisorCycleCache _snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); } - public ulong GetCycleLength(ulong divisor) + public ulong GetCycleLength(ulong divisor, bool skipPrimeOrderHeuristic = false) { - Span result = stackalloc ulong[1]; - ReadOnlySpan singleDivisor = stackalloc ulong[1] { divisor }; - GetCycleLengths(singleDivisor, result); - return result[0]; + if (!skipPrimeOrderHeuristic) + { + Span result = stackalloc ulong[1]; + ReadOnlySpan singleDivisor = stackalloc ulong[1] { divisor }; + GetCycleLengths(singleDivisor, result); + return result[0]; + } + + // The skipPrimeOrderHeuristic path is only used for divisors produced by the CPU by-divisor scanner, + // which always generates values greater than one. Keep this guard disabled here to avoid repeating the + // validation already performed by the general GetCycleLengths path. + // if (divisor <= 1UL) + // { + // throw new InvalidDataException("Divisor must be > 1"); + // } + + ulong[] snapshot = _snapshot; + if (divisor < (ulong)snapshot.Length) + { + ulong cached = snapshot[divisor]; + if (cached == 0UL) + { + throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); + } + + return cached; + } + + MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(divisor); + return MersenneDivisorCycles.CalculateCycleLength(divisor, divisorData, skipPrimeOrderHeuristic: true); } public void GetCycleLengths(ReadOnlySpan divisors, Span cycles) diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index dd72cdf..771fe28 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -94,7 +94,7 @@ public class MersenneDivisorCycles factorCache = cache; } - if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, cache, out ulong cycleLength) && cycleLength != 0UL) + if (TryCalculateCycleLengthForExponent(divisor, exponent, divisorData, cache, out ulong cycleLength, out bool _) && cycleLength != 0UL) { return cycleLength == exponent; } @@ -266,9 +266,11 @@ public class MersenneDivisorCycles ulong exponent, in MontgomeryDivisorData divisorData, Dictionary? factorCache, - out ulong cycleLength) + out ulong cycleLength, + out bool primeOrderFailed) { cycleLength = 0UL; + primeOrderFailed = false; if ((divisor & (divisor - 1UL)) == 0UL) { @@ -294,6 +296,8 @@ public class MersenneDivisorCycles return true; } + primeOrderFailed = true; + ulong phi = divisor - 1UL; // A zero totient would imply divisor == 0, which never occurs on the tested path. Keep the guard documented without // branching in the hot loop. @@ -332,6 +336,7 @@ public class MersenneDivisorCycles { cycleLength = ReduceOrder(divisorData, phi, factorCounts); success = true; + primeOrderFailed = false; } ThreadStaticPools.ReturnFactorCountDictionary(factorCounts); @@ -1024,15 +1029,15 @@ public class MersenneDivisorCycles return pow == 1UL; } - public static ulong CalculateCycleLength(ulong divisor, in MontgomeryDivisorData divisorData) + public static ulong CalculateCycleLength(ulong divisor, in MontgomeryDivisorData divisorData, bool skipPrimeOrderHeuristic = false) { - if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong cycleLength)) + if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong cycleLength, skipPrimeOrderHeuristic)) return cycleLength; return CalculateCycleLengthFallback(divisor); } - internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength) + internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength, bool skipPrimeOrderHeuristic = false) { if ((divisor & (divisor - 1UL)) == 0UL) { @@ -1046,7 +1051,7 @@ public class MersenneDivisorCycles return true; } - if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) + if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) { ulong computedOrder = PrimeOrderCalculator.Calculate( divisor, diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs index 9c60557..9cf94df 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.Cpu.cs @@ -56,17 +56,19 @@ internal static partial class PrimeOrderCalculator ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, divisorData); - if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config)) + PartialFactorResult? orderFactors; + if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, config, out orderFactors)) { return candidateOrder; } if (config.Mode == PrimeOrderMode.Strict) { + orderFactors?.Dispose(); return CalculateByFactorizationCpu(prime, divisorData); } - if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) + if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, orderFactors, out ulong order)) { return order; } @@ -160,80 +162,90 @@ internal static partial class PrimeOrderCalculator return order; } - private static bool TryConfirmOrderCpu(ulong prime, ulong order, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) + private static bool TryConfirmOrderCpu( + ulong prime, + ulong order, + in MontgomeryDivisorData divisorData, + in PrimeOrderSearchConfig config, + out PartialFactorResult? reusableFactorization) { + reusableFactorization = null; + 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) { - if (factorization.Factors is null) + factorization.Dispose(); + return false; + } + + if (!factorization.FullyFactored) + { + if (factorization.Cofactor <= 1UL) { + factorization.Dispose(); return false; } - if (!factorization.FullyFactored) + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); + if (!isPrime) { - 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; + return false; } - 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; - } + PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); + factorization.Dispose(); + factorization = extended; + } - reduced /= primeFactor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } + if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization)) + { + reusableFactorization = factorization; + return false; + } + + factorization.Dispose(); + return true; + } + + private static bool ValidateOrderAgainstFactors( + ulong prime, + ulong order, + in MontgomeryDivisorData divisorData, + PartialFactorResult factorization) + { + 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 (Pow2EqualsOneCpu(reduced, prime, divisorData)) + { + return false; } } + } - return true; - } - finally - { - factorization.Dispose(); - } + return true; } private static bool TryHeuristicFinishCpu( @@ -243,17 +255,18 @@ internal static partial class PrimeOrderCalculator in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, PartialFactorResult phiFactors, + PartialFactorResult? cachedOrderFactors, out ulong result) { result = 0UL; if (order <= 1UL) { + cachedOrderFactors?.Dispose(); return false; } - // TODO: Do we do partial factoring of order multiple times? - // DebugLog("Trying heuristic. Partial factoring order"); - PartialFactorResult orderFactors = PartialFactor(order, config); + // Reuse the partial factorization from TryConfirmOrderCpu when available. + PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactor(order, config); try { if (orderFactors.Factors is null) @@ -410,7 +423,6 @@ internal static partial class PrimeOrderCalculator orderFactors.Dispose(); } } - private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates) { ulong previous = previousOrder ?? 0UL;