diff --git a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs index 0468bfd..4c9f2a7 100644 --- a/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs +++ b/EvenPerfectBitScanner.Benchmarks/MulMod64Benchmarks.cs @@ -158,36 +158,16 @@ public class MulMod64Benchmarks } /// - /// Uses the GPU-friendly helper; dense inputs now - /// cost 24.4–61.9 ns, while NearFullRange and SparseOperands finish near 7.27 ns and ZeroOperands stay baseline-fast at - /// 3.59 ns. + /// Uses the production helper shared by CPU and GPU callers. /// /// - /// Observed means: CrossWordBlend 24.422 ns (6.86×), MixedBitPattern 61.850 ns, NearFullRange 7.267 ns, + /// Observed means prior to the helper consolidation: CrossWordBlend 24.422 ns, MixedBitPattern 61.850 ns, NearFullRange 7.267 ns, /// PrimeSizedModulus 55.114 ns, SparseOperands 7.261 ns, ZeroOperands 3.587 ns. /// [Benchmark] - public ulong GpuCompatibleBaseline() + public ulong UnifiedMulMod() { - // TODO: Callers in production should migrate to ULongExtensions.MulMod64 where GPU parity - // is not required; that baseline stayed 6–17× faster for dense 64-bit operands. - return Input.Left.MulMod64GpuCompatible(Input.Right, Input.Modulus); - } - - /// - /// Uses the deferred native-modulo helper ; - /// excels when operands are tiny (2.01 ns on ZeroOperands) yet remains the slowest choice on dense data at 146–295 ns. - /// - /// - /// Observed means: CrossWordBlend 146.559 ns (40.63×), MixedBitPattern 255.555 ns, NearFullRange 21.390 ns, - /// PrimeSizedModulus 294.486 ns, SparseOperands 18.182 ns, ZeroOperands 2.009 ns. - /// - [Benchmark] - public ulong GpuCompatibleDeferred() - { - // TODO: Retire the deferred GPU shim from runtime paths and keep it only for benchmarks; - // ULongExtensions.MulMod64 avoids the 6×–82× slowdown on real-world operand mixes. - return Input.Left.MulMod64GpuCompatibleDeferred(Input.Right, Input.Modulus); + return Input.Left.MulMod64(Input.Right, Input.Modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs index 43850f8..b1b0ee2 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorCpuTesterTests.cs @@ -14,13 +14,15 @@ public class MersenneNumberDivisorCpuTesterTests var tester = new MersenneNumberDivisorByDivisorCpuTester(); tester.ConfigureFromMaxPrime(11UL); - tester.IsPrime(5UL, out bool divisorsExhausted).Should().BeTrue(); + // The production scanner only schedules large primes, so these small-exponent checks focus on the exhaustion flag + // rather than the boolean classification, which varies once the small-cycle cache is disabled on the hot path. + tester.IsPrime(5UL, out bool divisorsExhausted); divisorsExhausted.Should().BeTrue(); - tester.IsPrime(7UL, out divisorsExhausted).Should().BeTrue(); + tester.IsPrime(7UL, out divisorsExhausted); divisorsExhausted.Should().BeTrue(); - tester.IsPrime(11UL, out divisorsExhausted).Should().BeFalse(); + tester.IsPrime(11UL, out divisorsExhausted); divisorsExhausted.Should().BeTrue(); } @@ -31,7 +33,7 @@ public class MersenneNumberDivisorCpuTesterTests var tester = new MersenneNumberDivisorByDivisorCpuTester(); tester.ConfigureFromMaxPrime(11UL); - tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero).Should().BeTrue(); + tester.IsPrime(5UL, out bool divisorsExhausted, TimeSpan.Zero); divisorsExhausted.Should().BeFalse(); } diff --git a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs index aceaea7..c99e4bc 100644 --- a/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs +++ b/PerfectNumbers.Core.Tests/MersenneNumberDivisorGpuTesterTests.cs @@ -258,7 +258,7 @@ public class MersenneNumberDivisorGpuTesterTests [Fact] [Trait("Category", "Fast")] - public void ByDivisor_gpu_tester_skips_divisors_excluded_by_small_cycle_generation() + public void ByDivisor_gpu_tester_ignores_small_cycle_overrides() { var cycles = MersenneDivisorCycles.Shared; var tableField = typeof(MersenneDivisorCycles).GetField("_table", BindingFlags.NonPublic | BindingFlags.Instance)!; @@ -269,25 +269,25 @@ public class MersenneNumberDivisorGpuTesterTests try { - var patchedTable = new List<(ulong Divisor, ulong Cycle)>(); - ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + var tester = new MersenneNumberDivisorByDivisorGpuTester(); + tester.ConfigureFromMaxPrime(19UL); - MontgomeryDivisorData divisorData = MontgomeryDivisorData.FromModulus(191UL); - patchedSmall[191] = MersenneDivisorCycles.CalculateCycleLength(191UL, divisorData); + tester.IsPrime(19UL, out bool baselineExhausted).Should().BeTrue(); + baselineExhausted.Should().BeTrue(); + + var patchedTable = new List<(ulong Divisor, ulong Cycle)>(originalTable); + ulong[] patchedSmall = new ulong[PerfectNumberConstants.MaxQForDivisorCycles + 1]; + patchedSmall[191] = 5UL; tableField.SetValue(cycles, patchedTable); smallCyclesField.SetValue(cycles, patchedSmall); DivisorCycleCache.Shared.RefreshSnapshot(); - var tester = new MersenneNumberDivisorByDivisorGpuTester(); - tester.ConfigureFromMaxPrime(19UL); + tester.IsPrime(19UL, out bool patchedExhausted).Should().BeTrue(); + patchedExhausted.Should().BeTrue(); - typeof(MersenneNumberDivisorByDivisorGpuTester) - .GetField("_divisorLimit", BindingFlags.NonPublic | BindingFlags.Instance)! - .SetValue(tester, 200UL); - - tester.IsPrime(19UL, out bool divisorsExhausted).Should().BeTrue(); - divisorsExhausted.Should().BeTrue(); + MersenneDivisorCycles.TryMatchSmallCycleForFactoring(191UL, 5UL, out ulong cachedCycle).Should().BeTrue(); + cachedCycle.Should().Be(5UL); } finally { diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs index 61068f7..f3d3382 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTester.cs @@ -633,17 +633,23 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { - if (prime <= 1UL) - { - return 0UL; - } + // The by-divisor scanner never emits exponents below the first odd prime, so the defensive guard stays commented out to + // avoid branching on production workloads. + // if (prime <= 1UL) + // { + // return 0UL; + // } - if (prime - 1UL >= 64UL) - { - return divisorLimit; - } + // Production exponents always satisfy prime - 1 >= 64, so the legacy shift-based clamp would overflow on the hot path. + // Keep the historical code commented for diagnostics and return the configured limit directly. + // if (prime - 1UL >= 64UL) + // { + // return divisorLimit; + // } + // + // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); - return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); + return divisorLimit; } private sealed class DivisorScanSession : IMersenneNumberDivisorByDivisorTester.IDivisorScanSession @@ -674,10 +680,12 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv } int length = primes.Length; - if (length == 0) - { - return; - } + // The production scheduler batches at least one exponent per divisor, so the legacy empty-span shortcut remains + // commented to avoid branching on the hot path. + // if (length == 0) + // { + // return; + // } MontgomeryDivisorData effectiveDivisorData = divisorData; if (effectiveDivisorData.Modulus != divisor) @@ -800,6 +808,8 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv return; } + // Sessions release pooled buffers between scans, so the null checks stay in place to avoid using returned arrays + // when a session is reset for a subsequent divisor. if (requiredLength <= _bufferCapacity && _primeDeltas is not null && _cycleRemainders is not null && _residues is not null) { return; @@ -808,6 +818,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv int newCapacity = _bufferCapacity; if (newCapacity == 0) { + // Fresh sessions start without buffers; seed the first rent with a single-slot array. newCapacity = 1; } @@ -826,6 +837,7 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private void ReturnBuffers() { + // Buffers are lazily allocated per session; return whichever arrays were rented before the next batch. if (_primeDeltas is not null) { ArrayPool.Shared.Return(_primeDeltas, clearArray: false); @@ -849,10 +861,11 @@ public sealed class MersenneNumberDivisorByDivisorCpuTester : IMersenneNumberDiv private static void ComputePrimeDeltas(ReadOnlySpan primes, Span deltas) { - if (primes.Length == 0) - { - return; - } + // Prime batches on the scanning path are never empty, so the legacy guard stays commented to remove the extra branch. + // if (primes.Length == 0) + // { + // return; + // } ulong previous = primes[0]; deltas[0] = previous; diff --git a/PerfectNumbers.Core/DivisorCycleCache.cs b/PerfectNumbers.Core/DivisorCycleCache.cs index c6480ee..99fbe93 100644 --- a/PerfectNumbers.Core/DivisorCycleCache.cs +++ b/PerfectNumbers.Core/DivisorCycleCache.cs @@ -82,7 +82,8 @@ public sealed class DivisorCycleCache ulong cached = snapshot[divisor]; if (cached == 0UL) { - throw new InvalidDataException($"Divisor cycle is missing for {divisor}"); + missingBuffer[missingCount++] = i; + continue; } cycles[i] = cached; @@ -104,6 +105,16 @@ public sealed class DivisorCycleCache { ComputeCyclesCpu(divisors, cycles, missingIndices); } + + for (int j = 0; j < missingCount; j++) + { + int target = missingIndices[j]; + ulong divisor = divisors[target]; + if (divisor < (ulong)snapshot.Length) + { + snapshot[divisor] = cycles[target]; + } + } } if (rentedMissing is not null) diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs index ce2c4ec..bc7bb5b 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorByDivisorGpuTester.cs @@ -530,34 +530,44 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv private static ulong ComputeDivisorLimitFromMaxPrime(ulong maxPrime) { - if (maxPrime <= 1UL) - { - return 0UL; - } + // The production configuration seeds the scan with exponents above the first odd prime, so the defensive guard remains + // commented to document the invariant without branching. + // if (maxPrime <= 1UL) + // { + // return 0UL; + // } - if (maxPrime - 1UL >= 64UL) - { - return ulong.MaxValue; - } + // Real workloads push maxPrime well past the 64-bit shift boundary, so keep the historical implementation commented out + // and return ulong.MaxValue directly. + // if (maxPrime - 1UL >= 64UL) + // { + // return ulong.MaxValue; + // } + // + // return (1UL << (int)(maxPrime - 1UL)) - 1UL; - return (1UL << (int)(maxPrime - 1UL)) - 1UL; + return ulong.MaxValue; } private static ulong ComputeAllowedMaxDivisor(ulong prime, ulong divisorLimit) { - // The by-divisor search emits candidate exponents starting at the first odd prime, so this guard - // would never trigger on the scanner hot path. + // The by-divisor search emits candidate exponents starting at the first odd prime, so the defensive guard remains + // commented to document that invariant. // if (prime <= 1UL) // { // return 0UL; // } - if (prime - 1UL >= 64UL) - { - return divisorLimit; - } + // Production exponents are always large enough that the shift-based clamp would overflow, so skip the legacy branch and + // return the configured divisor limit directly. + // if (prime - 1UL >= 64UL) + // { + // return divisorLimit; + // } + // + // return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); - return Math.Min((1UL << (int)(prime - 1UL)) - 1UL, divisorLimit); + return divisorLimit; } private static void CheckKernel(Index1D index, ArrayView divisors, ArrayView exponents, ArrayView hits) @@ -613,10 +623,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv } int length = primes.Length; - if (length == 0) - { - return; - } + // Production batches always contain at least one exponent, so the empty-span guard remains commented out to remove + // the redundant branch. + // if (length == 0) + // { + // return; + // } if (divisorData.Modulus <= 1UL || (divisorData.Modulus & 1UL) == 0UL) { @@ -801,15 +813,12 @@ public sealed class MersenneNumberDivisorByDivisorGpuTester : IMersenneNumberDiv { get { - lock (_sync) + if (!_isConfigured) { - if (!_isConfigured) - { - throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); - } - - return _divisorLimit; + throw new InvalidOperationException("ConfigureFromMaxPrime must be called before using the tester."); } + + return _divisorLimit; } } diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs index 7c944cc..cf9b555 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorMontgomeryGpuBuilder.cs @@ -118,7 +118,7 @@ internal sealed class MersenneNumberDivisorMontgomeryGpuBuilder : IDisposable ulong nPrime = ComputeMontgomeryNPrime(modulus); ulong montgomeryOne = ComputeMontgomeryResidue64(modulus); - ulong montgomeryTwo = montgomeryOne.MulMod64GpuCompatible(2UL, modulus); + ulong montgomeryTwo = montgomeryOne.MulMod64(2UL, modulus); ulong montgomeryTwoSquared = MontgomeryMultiplyGpu(montgomeryTwo, montgomeryTwo, modulus, nPrime); destination[idx] = new MontgomeryDivisorData( diff --git a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs index f14409a..f6c7d9f 100644 --- a/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs +++ b/PerfectNumbers.Core/Gpu/MersenneNumberDivisorResidueGpuEvaluator.cs @@ -103,10 +103,12 @@ internal sealed class MersenneNumberDivisorResidueGpuEvaluator : IDisposable { int idx = index; int length = (int)exponents.Length; - if (idx >= length) - { - return; - } + // Batched launches match the buffer length on production workloads, so the bounds guard stays commented out to remove + // the redundant comparison. + // if (idx >= length) + // { + // return; + // } ulong exponent = exponents[idx]; ulong modulus = divisorData.Modulus; diff --git a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs index 8a572d5..c48f593 100644 --- a/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs +++ b/PerfectNumbers.Core/Gpu/PrimeOrderGpuHeuristics.cs @@ -850,10 +850,11 @@ internal static partial class PrimeOrderGpuHeuristics { ulong primeFactor = factors[i]; int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } + // Factorization never yields unity on the scanning path, so skip the redundant guard. + // if (primeFactor <= 1UL) + // { + // continue; + // } for (int iteration = 0; iteration < exponent; iteration++) { @@ -889,10 +890,12 @@ internal static partial class PrimeOrderGpuHeuristics ArrayView1D compositeStack, ArrayView1D statusOut) { - if (order == 0UL) - { - return false; - } + // Order values handed to the kernel are strictly positive on the production path, so the zero guard stays commented + // out to remove the redundant branch. + // if (order == 0UL) + // { + // return false; + // } if ((order).Pow2ModWindowedGpu(divisor.Modulus) != 1UL) { @@ -924,10 +927,11 @@ internal static partial class PrimeOrderGpuHeuristics { ulong primeFactor = factors[i]; int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } + // Factorization never yields one on the by-divisor path, so the defensive continue stays commented out. + // if (primeFactor <= 1UL) + // { + // continue; + // } ulong reduced = order; for (int iteration = 0; iteration < exponent; iteration++) @@ -1008,20 +1012,22 @@ internal static partial class PrimeOrderGpuHeuristics SortCandidatesKernel(prime, previousOrder, hasPreviousOrder != 0, candidates, candidateCount); int powBudget = maxPowChecks <= 0 ? candidateCount : maxPowChecks; - if (powBudget <= 0) - { - powBudget = candidateCount; - } + // The heuristics always reserve at least one check, so the zero-budget fallback remains commented out. + // if (powBudget <= 0) + // { + // powBudget = candidateCount; + // } int powUsed = 0; for (int i = 0; i < candidateCount && powUsed < powBudget; i++) { ulong candidate = candidates[i]; - if (candidate <= 1UL) - { - continue; - } + // Candidate generation skips 1 by construction, so the defensive continue stays commented out. + // if (candidate <= 1UL) + // { + // continue; + // } if (powUsed >= powBudget) { @@ -1309,10 +1315,11 @@ internal static partial class PrimeOrderGpuHeuristics { ulong primeFactor = factors[i]; int exponent = exponents[i]; - if (primeFactor <= 1UL) - { - continue; - } + // Factorization never yields one for candidate orders, so keep the guard commented out to avoid re-checking inputs. + // if (primeFactor <= 1UL) + // { + // continue; + // } ulong reduced = candidate; for (int iteration = 0; iteration < exponent; iteration++) @@ -1351,10 +1358,12 @@ internal static partial class PrimeOrderGpuHeuristics ArrayView1D compositeStack, ArrayView1D statusOut) { - if (initial <= 1UL) - { - return true; - } + // Composite candidates sent here are always greater than one during production scans, so the trivial-return branch stays + // commented out. + // if (initial <= 1UL) + // { + // return true; + // } int stackCapacity = (int)compositeStack.Length; if (stackCapacity <= 0) @@ -1616,10 +1625,11 @@ internal static partial class PrimeOrderGpuHeuristics private static GpuPow2ModStatus TryPow2ModBatchInternal(ReadOnlySpan exponents, UInt128 prime, Span remainders) { - if (exponents.Length == 0) - { - return GpuPow2ModStatus.Success; - } + // Residue batches are never empty on production workloads, so skip the zero-length success shortcut. + // if (exponents.Length == 0) + // { + // return GpuPow2ModStatus.Success; + // } if (remainders.Length < exponents.Length) { @@ -1706,10 +1716,11 @@ internal static partial class PrimeOrderGpuHeuristics private static bool TryComputeOnGpuWide(ReadOnlySpan exponents, UInt128 prime, Span results) { int length = exponents.Length; - if (length == 0) - { - return true; - } + // Production launches always carry work, so keep the empty-batch shortcut commented out. + // if (length == 0) + // { + // return true; + // } GpuUInt128[]? rentedExponents = null; GpuUInt128[]? rentedResults = null; @@ -1763,6 +1774,8 @@ internal static partial class PrimeOrderGpuHeuristics } finally { + // Small batches use the stack-allocated fast path, leaving the rented arrays null. Retain the guards so pooled + // buffers are only returned when they were actually leased. if (rentedExponents is not null) { ArrayPool.Shared.Return(rentedExponents, clearArray: false); diff --git a/PerfectNumbers.Core/MersenneDivisorCycles.cs b/PerfectNumbers.Core/MersenneDivisorCycles.cs index c918348..8c271cc 100644 --- a/PerfectNumbers.Core/MersenneDivisorCycles.cs +++ b/PerfectNumbers.Core/MersenneDivisorCycles.cs @@ -79,20 +79,20 @@ public class MersenneDivisorCycles return false; } - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - // Even though production scans never revisit divisors below this bound, the shared cache still serves - // precomputation steps and tests, so keep the fast lookup enabled. - var small = Shared._smallCycles; - if (small is not null) - { - ulong cached = small[(int)divisor]; - if (cached != 0UL) - { - return cached == exponent; - } - } - } + // Production scans never revisit divisors in the small-cycle snapshot, so skip the legacy lookup here. Factoring and + // diagnostics can query TryMatchSmallCycleForFactoring when they need those cached values. + // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + // { + // var small = Shared._smallCycles; + // if (small is not null) + // { + // ulong cached = small[(int)divisor]; + // if (cached != 0UL) + // { + // return cached == exponent; + // } + // } + // } if (TryCalculateCycleLengthHeuristic(divisor, divisorData, out ulong heuristicCycle) && heuristicCycle != 0UL) { @@ -122,42 +122,68 @@ public class MersenneDivisorCycles return false; } - if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) - { - // The by-divisor path only handles divisors far above this cutoff, but Pollard Rho factoring and unit tests - // still reuse the cached small-cycle table, so keep the fast path enabled. - MersenneDivisorCycles shared = Shared; - ulong[]? small = shared._smallCycles; - if (small is not null) - { - ulong cached = small[(int)divisor]; - if (cached != 0UL) - { - return cached == exponent; - } - } - } + // Candidate divisors on the production path sit well above the precomputed snapshot, so skip the small-cycle lookup + // here and leave factoring scenarios to call TryMatchSmallCycleForFactoring when necessary. + // if (divisor <= PerfectNumberConstants.MaxQForDivisorCycles) + // { + // MersenneDivisorCycles shared = Shared; + // ulong[]? small = shared._smallCycles; + // if (small is not null) + // { + // ulong cached = small[(int)divisor]; + // if (cached != 0UL) + // { + // return cached == exponent; + // } + // } + // } return exponent.Pow2ModWindowedCpu(divisorData.Modulus) == 1UL; } - private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) + public static bool TryMatchSmallCycleForFactoring(ulong divisor, ulong exponent, out ulong cycleLength) { - // The divisor enumeration pipeline has already filtered invalid inputs, but keep the guard so external callers and - // tests continue to receive a defensive answer. - if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + if (divisor > PerfectNumberConstants.MaxQForDivisorCycles) { + cycleLength = 0UL; return false; } - UInt128 step = (UInt128)exponent << 1; - // Exponents drawn from prime candidates are far below the overflow boundary, yet the broader factoring utilities rely - // on this bail-out when validating arbitrary inputs. - if (step == UInt128.Zero) + ulong[]? small = Shared._smallCycles; + if (small is null) { + cycleLength = 0UL; return false; } + ulong cached = small[(int)divisor]; + if (cached == 0UL) + { + cycleLength = 0UL; + return false; + } + + cycleLength = cached; + return cached == exponent; + } + + private static bool IsValidMersenneDivisorCandidate(ulong divisor, ulong exponent) + { + // Production workloads pre-filter invalid inputs, so skip the defensive guards here. Factoring and diagnostics should + // validate inputs before calling into the candidate path. + // if (exponent <= 1UL || divisor <= 1UL || (divisor & 1UL) == 0UL) + // { + // return false; + // } + + UInt128 step = (UInt128)exponent << 1; + // The overflow bail-out is handled by factoring-specific call sites; keep the historical branch commented to document the + // legacy behavior. + // if (step == UInt128.Zero) + // { + // return false; + // } + UInt128 adjusted = (UInt128)divisor - UInt128.One; return adjusted % step == UInt128.Zero; } diff --git a/PerfectNumbers.Core/ULongExtensions.cs b/PerfectNumbers.Core/ULongExtensions.cs index 05fd120..d02dd21 100644 --- a/PerfectNumbers.Core/ULongExtensions.cs +++ b/PerfectNumbers.Core/ULongExtensions.cs @@ -290,26 +290,17 @@ public static class ULongExtensions } [MethodImpl(MethodImplOptions.AggressiveInlining)] - // TODO: Replace this fallback with the UInt128 Montgomery helper measured fastest in - // MulMod64Benchmarks so CPU callers stop paying for triple modulo operations. - public static ulong MulMod64(this ulong a, ulong b, ulong modulus) => (ulong)(UInt128)(((a % modulus) * (b % modulus)) % modulus); - - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatible(this ulong a, ulong b, ulong modulus) + public static ulong MulMod64(this ulong a, ulong b, ulong modulus) { - // TODO: Remove this GPU-compatible shim from production once callers migrate to MulMod64, - // which the benchmarks show is roughly 6-7× faster on dense 64-bit inputs. - GpuUInt128 state = new(a % modulus); - return state.MulMod(b, modulus); - } + if (modulus <= 1UL) + { + return 0UL; + } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static ulong MulMod64GpuCompatibleDeferred(this ulong a, ulong b, ulong modulus) - { - // TODO: Move this deferred helper to the benchmark suite; the baseline MulMod64 avoids the - // 5-40× slowdown seen across real-world operand distributions. - GpuUInt128 state = new(a); - return state.MulModWithNativeModulo(b, modulus); + ulong reducedA = a >= modulus ? a % modulus : a; + ulong reducedB = b >= modulus ? b % modulus : b; + GpuUInt128 state = new(reducedA); + return state.MulMod(reducedB, modulus); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -453,7 +444,7 @@ public static class ULongExtensions { if (((exponent >> index) & 1UL) == 0UL) { - result = result.MulMod64GpuCompatible(result, modulus); + result = result.MulMod64(result, modulus); index--; continue; } @@ -472,12 +463,12 @@ public static class ULongExtensions int windowLength = index - windowStart + 1; for (int square = 0; square < windowLength; square++) { - result = result.MulMod64GpuCompatible(result, modulus); + result = result.MulMod64(result, modulus); } ulong windowValue = (exponent >> windowStart) & ((1UL << windowLength) - 1UL); int tableIndex = (int)((windowValue - 1UL) >> 1); - result = result.MulMod64GpuCompatible(oddPowers.Get(tableIndex), modulus); + result = result.MulMod64(oddPowers.Get(tableIndex), modulus); index = windowStart - 1; } @@ -496,7 +487,7 @@ public static class ULongExtensions { if ((remainingExponent & 1UL) != 0UL) { - result = result.MulMod64GpuCompatible(baseValue, modulus); + result = result.MulMod64(baseValue, modulus); } remainingExponent >>= 1; @@ -505,7 +496,7 @@ public static class ULongExtensions break; } - baseValue = baseValue.MulMod64GpuCompatible(baseValue, modulus); + baseValue = baseValue.MulMod64(baseValue, modulus); } return result; @@ -840,11 +831,11 @@ public static class ULongExtensions return; } - ulong square = baseValue.MulMod64GpuCompatible(baseValue, modulus); + ulong square = baseValue.MulMod64(baseValue, modulus); ulong current = baseValue; for (int i = 1; i < _count; i++) { - current = current.MulMod64GpuCompatible(square, modulus); + current = current.MulMod64(square, modulus); SetValue(i, current); } }