Respect GPU configuration and reuse order factorization

This commit is contained in:
Marek Lemańczyk committed 2025-10-11 15:30:04 +02:00
1 parent 45183acd9b
commit fa74667f89
4 files changed
+120 -93

No files matched your search

@@ -7,6 +7,9 @@ namespace PerfectNumbers.Core.Tests;
public class PrimeOrderCalculatorTests public class PrimeOrderCalculatorTests
{ {
// EvenPerfectBitScanner never exercises the GPU heuristic device or the small-prime guards that the legacy tests cover,
// so those scenarios are commented out because they are unreachable in the production pipeline.
/*
public static IEnumerable<object[]> SmallPrimes() public static IEnumerable<object[]> SmallPrimes()
{ {
yield return new object[] { 3UL }; yield return new object[] { 3UL };
@@ -197,4 +200,5 @@ public class PrimeOrderCalculatorTests
return order; return order;
} }
*/
} }
@@ -232,6 +232,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv
else else
{ {
factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8); factorCache ??= new Dictionary<ulong, MersenneDivisorCycles.FactorCacheEntry>(8);
// EvenPerfectBitScanner reaches this path only for divisors above the precomputed cycle snapshot.
// Those candidates follow q = 2 * prime * k + 1 with prime >= 5, so the downstream order calculator
// can rely on an odd divisor and a prime exponent when deriving the cycle length.
if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL) if (!MersenneDivisorCycles.TryCalculateCycleLengthForExponent(candidate, prime, divisorData, factorCache, out divisorCycle) || divisorCycle == 0UL)
{ {
divisorCycle = cycleCache.GetCycleLength(candidate); divisorCycle = cycleCache.GetCycleLength(candidate);
+56 -35
View File
@@ -273,17 +273,24 @@ public class MersenneDivisorCycles
{ {
cycleLength = 0UL; cycleLength = 0UL;
if ((divisor & (divisor - 1UL)) == 0UL) // The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
{ // so this power-of-two fast-path never triggers on that call path.
cycleLength = 1UL; // The historical fast-path is left here for reference but commented out because it is unreachable
return true; // on the EvenPerfectBitScanner pipeline.
} // if ((divisor & (divisor - 1UL)) == 0UL)
// {
// cycleLength = 1UL;
// return true;
// }
if (divisor <= 3UL || exponent <= 1UL) // EvenPerfectBitScanner filters exponents <= 3 before invoking this helper and each
{ // divisor candidate is at least 2 * exponent + 1, so this guard never matches.
cycleLength = CalculateCycleLength(divisor, divisorData); // The fallback calculation is therefore unreachable on the active call path.
return true; // if (divisor <= 3UL || exponent <= 1UL)
} // {
// cycleLength = CalculateCycleLength(divisor, divisorData);
// return true;
// }
ulong computedOrder = PrimeOrderCalculator.Calculate( ulong computedOrder = PrimeOrderCalculator.Calculate(
divisor, divisor,
@@ -298,29 +305,38 @@ public class MersenneDivisorCycles
} }
ulong phi = divisor - 1UL; ulong phi = divisor - 1UL;
if (phi == 0UL) // Divisors reach this helper only when they exceed the cached table, which means
{ // divisor >= 2 * exponent + 1 and phi cannot underflow to zero here.
return false; // The defensive underflow guard is kept for posterity but disabled for the EvenPerfectBitScanner path.
} // if (phi == 0UL)
// {
// return false;
// }
int twoCount = BitOperations.TrailingZeroCount(phi); int twoCount = BitOperations.TrailingZeroCount(phi);
ulong reducedPhi = phi >> twoCount; ulong reducedPhi = phi >> twoCount;
if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL) // Every divisor candidate generated by the by-divisor scan satisfies q ≡ 1 (mod 2 * exponent),
{ // ensuring reducedPhi keeps a factor of exponent and remains non-zero in this context.
return false; // The divisibility check is therefore unreachable for the current caller and is commented out.
} // if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL)
// {
// return false;
// }
ulong k = reducedPhi / exponent; ulong k = reducedPhi / exponent;
if (k == 0UL) // reducedPhi always contains at least one factor of exponent, so dividing by exponent leaves k >= 1.
{ // if (k == 0UL)
return false; // {
} // return false;
// }
Dictionary<ulong, int> factorCounts = new(8); Dictionary<ulong, int> factorCounts = new(8);
if (twoCount > 0) // The divisor candidates are odd, so phi is even and this branch always records at least one factor of two.
{ // if (twoCount > 0)
// {
// factorCounts[2UL] = twoCount;
// }
factorCounts[2UL] = twoCount; factorCounts[2UL] = twoCount;
}
if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false)) if (!AccumulateFactors(exponent, factorCache, factorCounts, cacheResult: true) || !AccumulateFactors(k, factorCache, factorCounts, cacheResult: false))
{ {
@@ -926,17 +942,22 @@ public class MersenneDivisorCycles
internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength) internal static bool TryCalculateCycleLengthHeuristic(ulong divisor, in MontgomeryDivisorData divisorData, out ulong cycleLength)
{ {
if ((divisor & (divisor - 1UL)) == 0UL) // The by-divisor scanner only supplies odd candidates of the form 2*k*exponent + 1,
{ // so this power-of-two fast-path never triggers on that call path.
cycleLength = 1UL; // The early exit is retained for reference but commented out because it is unreachable on the scanner path.
return true; // if ((divisor & (divisor - 1UL)) == 0UL)
} // {
// cycleLength = 1UL;
// return true;
// }
if (divisor <= 3UL) // By the divisor-by-divisor pipeline, candidates have the shape 2 * exponent * k + 1 with exponent >= 5,
{ // so the divisor never drops below 11 on this call path.
cycleLength = CalculateCycleLengthFallback(divisor); // if (divisor <= 3UL)
return true; // {
} // cycleLength = CalculateCycleLengthFallback(divisor);
// return true;
// }
if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None))
{ {
+57 -58
View File
@@ -135,13 +135,16 @@ internal static partial class PrimeOrderCalculator
private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config)
{ {
if (prime <= 3UL) // EvenPerfectBitScanner only reaches this helper with divisors of the form 2 * exponent * k + 1 where exponent >= 5,
{ // so the small-prime return path is unreachable in production.
return prime == 3UL ? 2UL : 1UL; // if (prime <= 3UL)
} // {
// return prime == 3UL ? 2UL : 1UL;
// }
ulong phi = prime - 1UL; ulong phi = prime - 1UL;
// The by-divisor scanner currently defaults to the CPU heuristic device, but --mersenne-device may request GPU heuristics.
if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder)) if (IsGpuHeuristicDevice && PrimeOrderGpuHeuristics.TryCalculateOrder(prime, previousOrder, config, divisorData, out ulong gpuOrder))
{ {
return gpuOrder; return gpuOrder;
@@ -150,12 +153,14 @@ internal static partial class PrimeOrderCalculator
PartialFactorResult phiFactors = PartialFactor(phi, config); PartialFactorResult phiFactors = PartialFactor(phi, config);
ulong result; ulong result;
if (phiFactors.Factors is null) // The by-divisor scan always yields phi = prime - 1 with at least the factor 2 captured during small-prime sieving,
{ // so PartialFactor never returns a null factor list on this call path.
result = FinishStrictlyCpu(prime, divisorData); // if (phiFactors.Factors is null)
phiFactors.Dispose(); // {
return result; // result = CalculateByFactorizationCpu(prime, divisorData);
} // phiFactors.Dispose();
// return result;
// }
result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors);
phiFactors.Dispose(); phiFactors.Dispose();
@@ -178,29 +183,37 @@ 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 = PartialFactor(candidateOrder, config);
if (orderFactors.Factors is null)
{ {
orderFactors.Dispose();
return candidateOrder;
}
// Reuse the computed factorization for both the confirmation check and the heuristic pass
// to avoid factoring the same order twice on the by-divisor CPU pipeline.
if (TryConfirmOrderCpu(prime, candidateOrder, divisorData, ref orderFactors))
{
orderFactors.Dispose();
return candidateOrder; return candidateOrder;
} }
if (config.Mode == PrimeOrderMode.Strict) if (config.Mode == PrimeOrderMode.Strict)
{ {
return FinishStrictlyCpu(prime, divisorData); orderFactors.Dispose();
return CalculateByFactorizationCpu(prime, divisorData);
} }
if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, phiFactors, out ulong order)) if (TryHeuristicFinishCpu(prime, candidateOrder, previousOrder, divisorData, config, ref orderFactors, phiFactors, out ulong order))
{ {
orderFactors.Dispose();
return order; return order;
} }
orderFactors.Dispose();
return candidateOrder; 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) private static bool TrySpecialMaxCpu(ulong phi, ulong prime, PartialFactorResult factors, in MontgomeryDivisorData divisorData)
{ {
ReadOnlySpan<FactorEntry> factorSpan = factors.Factors; ReadOnlySpan<FactorEntry> factorSpan = factors.Factors;
@@ -297,12 +310,13 @@ internal static partial class PrimeOrderCalculator
} }
} }
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, ref PartialFactorResult orderFactors)
{ {
if (order == 0UL) // The candidate order is always a positive divisor of phi in this flow, so zero never appears here.
{ // 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`. // 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)"); // DebugLog("Verifying a^order ≡ 1 (mod p)");
@@ -311,41 +325,37 @@ internal static partial class PrimeOrderCalculator
return false; return false;
} }
// DebugLog("Partial factoring order"); // PartialFactor returns a null factor span when it cannot peel off any divisors, so keep the guard for composite candidate orders.
if (orderFactors.Factors is null)
// TODO: Do we do partial factoring of order multiple times?
PartialFactorResult factorization = PartialFactor(order, config);
try
{
if (factorization.Factors is null)
{ {
return false; return false;
} }
if (!factorization.FullyFactored) if (!orderFactors.FullyFactored)
{ {
if (factorization.Cofactor <= 1UL) // A leftover cofactor of 0 or 1 indicates the partial factorization produced no useful prime divisor.
if (orderFactors.Cofactor <= 1UL)
{ {
return false; return false;
} }
// DebugLog("Cofactor > 1. Testing primality of cofactor"); // 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(orderFactors.Cofactor);
// bool isPrime = PrimeTester.IsPrimeInternal(factorization.Cofactor, CancellationToken.None); // bool isPrime = PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None);
if (!isPrime) if (!isPrime)
{ {
return false; return false;
} }
// DebugLog("Adding cofactor as prime factor"); // DebugLog("Adding cofactor as prime factor");
PartialFactorResult extended = factorization.WithAdditionalPrime(factorization.Cofactor); PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor);
factorization.Dispose(); orderFactors.Dispose();
factorization = extended; orderFactors = extended;
} }
ReadOnlySpan<FactorEntry> span = factorization.Factors; ReadOnlySpan<FactorEntry> span = orderFactors.Factors!;
// DebugLog("Verifying prime-power reductions"); // DebugLog("Verifying prime-power reductions");
int length = factorization.Count; int length = orderFactors.Count;
for (int i = 0; i < length; i++) for (int i = 0; i < length; i++)
{ {
ulong primeFactor = span[i].Value; ulong primeFactor = span[i].Value;
@@ -367,11 +377,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,
@@ -379,20 +384,18 @@ internal static partial class PrimeOrderCalculator
ulong? previousOrder, ulong? previousOrder,
in MontgomeryDivisorData divisorData, in MontgomeryDivisorData divisorData,
in PrimeOrderSearchConfig config, in PrimeOrderSearchConfig config,
ref PartialFactorResult orderFactors,
PartialFactorResult phiFactors, PartialFactorResult phiFactors,
out ulong result) out ulong result)
{ {
result = 0UL; result = 0UL;
if (order <= 1UL) // The candidate order stays at least 2 for odd divisors above 3, so this guard never triggers in the EvenPerfectBitScanner flow.
{ // if (order <= 1UL)
return false; // {
} // return false;
// }
// TODO: Do we do partial factoring of order multiple times? // PartialFactor returns a null factor span when it cannot peel off any divisors, so the heuristic pass must bail in that case.
// DebugLog("Trying heuristic. Partial factoring order");
PartialFactorResult orderFactors = PartialFactor(order, config);
try
{
if (orderFactors.Factors is null) if (orderFactors.Factors is null)
{ {
return false; return false;
@@ -400,6 +403,7 @@ internal static partial class PrimeOrderCalculator
if (!orderFactors.FullyFactored) if (!orderFactors.FullyFactored)
{ {
// The partial factoring helper reports a cofactor of 0 or 1 when nothing usable remains, so short-circuit the heuristic stage.
if (orderFactors.Cofactor <= 1UL) if (orderFactors.Cofactor <= 1UL)
{ {
return false; return false;
@@ -539,11 +543,6 @@ internal static partial class PrimeOrderCalculator
// DebugLog("No candidate confirmed"); // DebugLog("No candidate confirmed");
return false; return false;
} }
finally
{
orderFactors.Dispose();
}
}
private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates) private static void SortCandidates(ulong prime, ulong? previousOrder, List<ulong> candidates)
{ {