Clarify skip heuristic guard and tighten order factor init

This commit is contained in:
Marek Lemańczyk committed 2025-10-17 15:05:11 +02:00
1 parent 43cf4a978f
commit 8428e6c5bd
4 files changed
+77 -33

No files matched your search

@@ -204,9 +204,10 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
prime, prime,
divisorData, divisorData,
factorCache, 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 else
{ {
+27 -1
View File
@@ -45,7 +45,9 @@ public sealed class DivisorCycleCache
_snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot(); _snapshot = MersenneDivisorCycles.Shared.ExportSmallCyclesSnapshot();
} }
public ulong GetCycleLength(ulong divisor) public ulong GetCycleLength(ulong divisor, bool skipPrimeOrderHeuristic = false)
{
if (!skipPrimeOrderHeuristic)
{ {
Span<ulong> result = stackalloc ulong[1]; Span<ulong> result = stackalloc ulong[1];
ReadOnlySpan<ulong> singleDivisor = stackalloc ulong[1] { divisor }; ReadOnlySpan<ulong> singleDivisor = stackalloc ulong[1] { divisor };
@@ -53,6 +55,30 @@ public sealed class DivisorCycleCache
return result[0]; 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<ulong> divisors, Span<ulong> cycles) public void GetCycleLengths(ReadOnlySpan<ulong> divisors, Span<ulong> cycles)
{ {
ulong[] snapshot = _snapshot; ulong[] snapshot = _snapshot;
+11 -6
View File
@@ -94,7 +94,7 @@ public class MersenneDivisorCycles
factorCache = cache; 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; return cycleLength == exponent;
} }
@@ -266,9 +266,11 @@ public class MersenneDivisorCycles
ulong exponent, ulong exponent,
in MontgomeryDivisorData divisorData, in MontgomeryDivisorData divisorData,
Dictionary<ulong, FactorCacheEntry>? factorCache, Dictionary<ulong, FactorCacheEntry>? factorCache,
out ulong cycleLength) out ulong cycleLength,
out bool primeOrderFailed)
{ {
cycleLength = 0UL; cycleLength = 0UL;
primeOrderFailed = false;
if ((divisor & (divisor - 1UL)) == 0UL) if ((divisor & (divisor - 1UL)) == 0UL)
{ {
@@ -294,6 +296,8 @@ public class MersenneDivisorCycles
return true; return true;
} }
primeOrderFailed = true;
ulong phi = divisor - 1UL; ulong phi = divisor - 1UL;
// A zero totient would imply divisor == 0, which never occurs on the tested path. Keep the guard documented without // A zero totient would imply divisor == 0, which never occurs on the tested path. Keep the guard documented without
// branching in the hot loop. // branching in the hot loop.
@@ -332,6 +336,7 @@ public class MersenneDivisorCycles
{ {
cycleLength = ReduceOrder(divisorData, phi, factorCounts); cycleLength = ReduceOrder(divisorData, phi, factorCounts);
success = true; success = true;
primeOrderFailed = false;
} }
ThreadStaticPools.ReturnFactorCountDictionary(factorCounts); ThreadStaticPools.ReturnFactorCountDictionary(factorCounts);
@@ -1024,15 +1029,15 @@ public class MersenneDivisorCycles
return pow == 1UL; 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 cycleLength;
return CalculateCycleLengthFallback(divisor); 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) if ((divisor & (divisor - 1UL)) == 0UL)
{ {
@@ -1046,7 +1051,7 @@ public class MersenneDivisorCycles
return true; return true;
} }
if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) if (!skipPrimeOrderHeuristic && PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
{ {
ulong computedOrder = PrimeOrderCalculator.Calculate( ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor, divisor,
+36 -24
View File
@@ -56,17 +56,19 @@ internal static partial class PrimeOrderCalculator
ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData); ulong candidateOrder = InitializeStartingOrderCpu(prime, phi, divisorData);
candidateOrder = ExponentLoweringCpu(candidateOrder, prime, phiFactors, 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; return candidateOrder;
} }
if (config.Mode == PrimeOrderMode.Strict) if (config.Mode == PrimeOrderMode.Strict)
{ {
orderFactors?.Dispose();
return CalculateByFactorizationCpu(prime, divisorData); 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; return order;
} }
@@ -160,28 +162,29 @@ internal static partial class PrimeOrderCalculator
return order; 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) if (order == 0UL)
{ {
return false; 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)) if (!Pow2EqualsOneCpu(order, prime, divisorData))
{ {
return false; return false;
} }
// DebugLog("Partial factoring order");
// TODO: Do we do partial factoring of order multiple times?
PartialFactorResult factorization = PartialFactor(order, config); PartialFactorResult factorization = PartialFactor(order, config);
try
{
if (factorization.Factors is null) if (factorization.Factors is null)
{ {
factorization.Dispose();
return false; return false;
} }
@@ -189,25 +192,39 @@ internal static partial class PrimeOrderCalculator
{ {
if (factorization.Cofactor <= 1UL) if (factorization.Cofactor <= 1UL)
{ {
factorization.Dispose();
return false; return false;
} }
// DebugLog("Cofactor > 1. Testing primality of cofactor");
bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor); bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(factorization.Cofactor);
// bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None);
if (!isPrime) if (!isPrime)
{ {
factorization.Dispose();
return false; return false;
} }
// DebugLog("Adding cofactor as prime factor");
PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor);
factorization.Dispose(); factorization.Dispose();
factorization = extended; factorization = extended;
} }
ReadOnlySpan<FactorEntry> span = factorization.Factors; if (!ValidateOrderAgainstFactors(prime, order, divisorData, factorization))
// DebugLog("Verifying prime-power reductions"); {
reusableFactorization = factorization;
return false;
}
factorization.Dispose();
return true;
}
private static bool ValidateOrderAgainstFactors(
ulong prime,
ulong order,
in MontgomeryDivisorData divisorData,
PartialFactorResult factorization)
{
ReadOnlySpan<FactorEntry> span = factorization.Factors!;
int length = factorization.Count; int length = factorization.Count;
for (int i = 0; i < length; i++) for (int i = 0; i < length; i++)
{ {
@@ -230,11 +247,6 @@ internal static partial class PrimeOrderCalculator
return true; return true;
} }
finally
{
factorization.Dispose();
}
}
private static bool TryHeuristicFinishCpu( private static bool TryHeuristicFinishCpu(
ulong prime, ulong prime,
@@ -243,17 +255,18 @@ internal static partial class PrimeOrderCalculator
in MontgomeryDivisorData divisorData, in MontgomeryDivisorData divisorData,
in PrimeOrderSearchConfig config, in PrimeOrderSearchConfig config,
PartialFactorResult phiFactors, PartialFactorResult phiFactors,
PartialFactorResult? cachedOrderFactors,
out ulong result) out ulong result)
{ {
result = 0UL; result = 0UL;
if (order <= 1UL) if (order <= 1UL)
{ {
cachedOrderFactors?.Dispose();
return false; return false;
} }
// TODO: Do we do partial factoring of order multiple times? // Reuse the partial factorization from TryConfirmOrderCpu when available.
// DebugLog("Trying heuristic. Partial factoring order"); PartialFactorResult orderFactors = cachedOrderFactors ?? PartialFactor(order, config);
PartialFactorResult orderFactors = PartialFactor(order, config);
try try
{ {
if (orderFactors.Factors is null) if (orderFactors.Factors is null)
@@ -410,7 +423,6 @@ internal static partial class PrimeOrderCalculator
orderFactors.Dispose(); orderFactors.Dispose();
} }
} }
private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates) private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates)
{ {
ulong previous = previousOrder ?? 0UL; ulong previous = previousOrder ?? 0UL;