diff --git a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs index 843131c..ac6c7f0 100644 --- a/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs +++ b/PerfectNumbers.Core.Tests/PrimeOrderCalculatorTests.cs @@ -7,6 +7,9 @@ namespace PerfectNumbers.Core.Tests; 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 SmallPrimes() { yield return new object[] { 3UL }; @@ -197,4 +200,5 @@ public class PrimeOrderCalculatorTests return order; } +*/ } diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 49baab4..f5862af 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -232,6 +232,9 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv else { factorCache ??= new Dictionary(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) { divisorCycle = cycleCache.GetCycleLength(candidate); diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index 410978a..5efc86a 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -273,17 +273,24 @@ public class MersenneDivisorCycles { cycleLength = 0UL; - if ((divisor & (divisor - 1UL)) == 0UL) - { - cycleLength = 1UL; - return true; - } + // 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. + // The historical fast-path is left here for reference but commented out because it is unreachable + // on the EvenPerfectBitScanner pipeline. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // cycleLength = 1UL; + // return true; + // } - if (divisor <= 3UL || exponent <= 1UL) - { - cycleLength = CalculateCycleLength(divisor, divisorData); - return true; - } + // EvenPerfectBitScanner filters exponents <= 3 before invoking this helper and each + // divisor candidate is at least 2 * exponent + 1, so this guard never matches. + // The fallback calculation is therefore unreachable on the active call path. + // if (divisor <= 3UL || exponent <= 1UL) + // { + // cycleLength = CalculateCycleLength(divisor, divisorData); + // return true; + // } ulong computedOrder = PrimeOrderCalculator.Calculate( divisor, @@ -298,29 +305,38 @@ public class MersenneDivisorCycles } ulong phi = divisor - 1UL; - if (phi == 0UL) - { - return false; - } + // Divisors reach this helper only when they exceed the cached table, which means + // divisor >= 2 * exponent + 1 and phi cannot underflow to zero here. + // The defensive underflow guard is kept for posterity but disabled for the EvenPerfectBitScanner path. + // if (phi == 0UL) + // { + // return false; + // } int twoCount = BitOperations.TrailingZeroCount(phi); ulong reducedPhi = phi >> twoCount; - if (reducedPhi == 0UL || (reducedPhi % exponent) != 0UL) - { - return false; - } + // 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. + // 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; - if (k == 0UL) - { - return false; - } + // reducedPhi always contains at least one factor of exponent, so dividing by exponent leaves k >= 1. + // if (k == 0UL) + // { + // return false; + // } Dictionary factorCounts = new(8); - if (twoCount > 0) - { - factorCounts[2UL] = twoCount; - } + // 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; 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) { - if ((divisor & (divisor - 1UL)) == 0UL) - { - cycleLength = 1UL; - return true; - } + // 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. + // The early exit is retained for reference but commented out because it is unreachable on the scanner path. + // if ((divisor & (divisor - 1UL)) == 0UL) + // { + // cycleLength = 1UL; + // return true; + // } - if (divisor <= 3UL) - { - cycleLength = CalculateCycleLengthFallback(divisor); - return true; - } + // 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. + // if (divisor <= 3UL) + // { + // cycleLength = CalculateCycleLengthFallback(divisor); + // return true; + // } if (PrimeTester.IsPrimeInternal(divisor, CancellationToken.None)) { diff --git a/PerfectNumbers.Core/PrimeOrderCalculator.cs b/PerfectNumbers.Core/PrimeOrderCalculator.cs index edc4f35..c89155e 100644 --- a/PerfectNumbers.Core/PrimeOrderCalculator.cs +++ b/PerfectNumbers.Core/PrimeOrderCalculator.cs @@ -135,13 +135,16 @@ internal static partial class PrimeOrderCalculator private static ulong CalculateInternal(ulong prime, ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config) { - if (prime <= 3UL) - { - return prime == 3UL ? 2UL : 1UL; - } + // 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. + // if (prime <= 3UL) + // { + // return prime == 3UL ? 2UL : 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)) { return gpuOrder; @@ -150,12 +153,14 @@ internal static partial class PrimeOrderCalculator PartialFactorResult phiFactors = PartialFactor(phi, config); ulong result; - if (phiFactors.Factors is null) - { - result = FinishStrictlyCpu(prime, divisorData); - phiFactors.Dispose(); - return result; - } + // 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. + // if (phiFactors.Factors is null) + // { + // result = CalculateByFactorizationCpu(prime, divisorData); + // phiFactors.Dispose(); + // return result; + // } result = RunHeuristicPipelineCpu(prime, previousOrder, config, divisorData, phi, phiFactors); phiFactors.Dispose(); @@ -178,29 +183,37 @@ 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 = 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; } 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; } + orderFactors.Dispose(); 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; @@ -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) - { - return false; - } + // The candidate order is always a positive divisor of phi in this flow, so zero never appears here. + // 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)"); @@ -311,66 +325,57 @@ internal static partial class PrimeOrderCalculator return false; } - // DebugLog("Partial factoring order"); - - // TODO: Do we do partial factoring of order multiple times? - PartialFactorResult factorization = PartialFactor(order, config); - try + // 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) { - if (factorization.Factors is null) + return false; + } + + if (!orderFactors.FullyFactored) + { + // A leftover cofactor of 0 or 1 indicates the partial factorization produced no useful prime divisor. + if (orderFactors.Cofactor <= 1UL) { return false; } - if (!factorization.FullyFactored) + // DebugLog("Cofactor > 1. Testing primality of cofactor"); + bool isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor); + // bool isPrime = PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None); + 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; - } - - reduced /= primeFactor; - if (Pow2EqualsOneCpu(reduced, prime, divisorData)) - { - return false; - } - } - } - - return true; + // DebugLog("Adding cofactor as prime factor"); + PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + orderFactors.Dispose(); + orderFactors = extended; } - finally + + ReadOnlySpan span = orderFactors.Factors!; + // DebugLog("Verifying prime-power reductions"); + int length = orderFactors.Count; + for (int i = 0; i < length; i++) { - factorization.Dispose(); + 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; } private static bool TryHeuristicFinishCpu( @@ -379,170 +384,164 @@ internal static partial class PrimeOrderCalculator ulong? previousOrder, in MontgomeryDivisorData divisorData, in PrimeOrderSearchConfig config, + ref PartialFactorResult orderFactors, PartialFactorResult phiFactors, out ulong result) { 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; + // } + + // PartialFactor returns a null factor span when it cannot peel off any divisors, so the heuristic pass must bail in that case. + if (orderFactors.Factors is null) { 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.FullyFactored) { - if (orderFactors.Factors is null) + // 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) { 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 = new(capacity); - FactorEntry[] factorArray = orderFactors.Factors!; - // DebugLog("Building candidates list"); - BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); - if (candidates.Count == 0) + var isPrime = Open.Numeric.Primes.Prime.Numbers.IsPrime(orderFactors.Cofactor); + // if (!PrimeTester.IsPrimeInternal(orderFactors.Cofactor, CancellationToken.None)) + if (!isPrime) { return false; } - // DebugLog("Sorting candidates"); - SortCandidates(prime, previousOrder, candidates); + PartialFactorResult extended = orderFactors.WithAdditionalPrime(orderFactors.Cofactor); + orderFactors.Dispose(); + orderFactors = extended; + } - int powBudget = config.MaxPowChecks <= 0 ? candidates.Count : config.MaxPowChecks; - int powUsed = 0; - int candidateCount = candidates.Count; - bool allowGpuBatch = true; - Span candidateSpan = CollectionsMarshal.AsSpan(candidates); + int capacity = config.MaxPowChecks <= 0 ? 64 : config.MaxPowChecks * 4; + List candidates = new(capacity); + FactorEntry[] factorArray = orderFactors.Factors!; + // DebugLog("Building candidates list"); + BuildCandidates(order, factorArray, orderFactors.Count, candidates, capacity); + if (candidates.Count == 0) + { + return false; + } - // DebugLog(() => $"Checking candidates ({candidateCount} candidates, {powBudget} pow budget)"); - int index = 0; - const int MaxGpuBatchSize = 1024; - const int StackGpuBatchSize = 1024; - Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize]; - ArrayPool uLongPool = ULongPool; - while (index < candidateCount && powUsed < powBudget) + // 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 = 1024; + const int StackGpuBatchSize = 1024; + Span stackGpuRemainders = stackalloc ulong[StackGpuBatchSize]; + ArrayPool uLongPool = 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) { - 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) { - 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) { - Span localRemainders = stackGpuRemainders.Slice(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, localRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - gpuSuccess = true; - gpuStackRemainders = true; - } - } - else - { - gpuPool = uLongPool.Rent(batchSize); - Span pooledRemainders = gpuPool.AsSpan(0, batchSize); - status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); - if (status == GpuPow2ModStatus.Success) - { - pooledGpuRemainders = pooledRemainders; - gpuSuccess = true; - } - else - { - uLongPool.Return(gpuPool, clearArray: false); - gpuPool = null; - } - } - - if (!gpuSuccess && (status == GpuPow2ModStatus.Overflow || status == GpuPow2ModStatus.Unavailable)) - { - allowGpuBatch = false; + gpuSuccess = true; + gpuStackRemainders = true; } } - - for (int i = 0; i < batchSize && powUsed < powBudget; i++) + else { - ulong candidate = batch[i]; - powUsed++; - - bool equalsOne; - if (gpuSuccess) + gpuPool = uLongPool.Rent(batchSize); + Span pooledRemainders = gpuPool.AsSpan(0, batchSize); + status = PrimeOrderGpuHeuristics.TryPow2ModBatch(batch, prime, pooledRemainders, divisorData); + if (status == GpuPow2ModStatus.Success) { - ulong remainderValue = gpuStackRemainders ? stackGpuRemainders[i] : pooledGpuRemainders[i]; - equalsOne = remainderValue == 1UL; + pooledGpuRemainders = pooledRemainders; + gpuSuccess = true; } else - { - equalsOne = Pow2EqualsOneCpu(candidate, prime, divisorData); - } - if (!equalsOne) - { - continue; - } - - if (!TryConfirmCandidateCpu(prime, candidate, divisorData, config, ref powUsed, powBudget)) - { - continue; - } - - if (gpuPool is not null) { uLongPool.Return(gpuPool, clearArray: false); + gpuPool = null; } + } - result = candidate; - return true; + 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) { - uLongPool.Return(gpuPool, clearArray: false); + uLongPool.Return(gpuPool, clearArray: false); } - index += batchSize; + result = candidate; + return true; } - // DebugLog("No candidate confirmed"); - return false; - } - finally - { - orderFactors.Dispose(); + if (gpuPool is not null) + { + uLongPool.Return(gpuPool, clearArray: false); + } + + index += batchSize; } + + // DebugLog("No candidate confirmed"); + return false; } private static void SortCandidates(ulong prime, ulong? previousOrder, List candidates)