From b2a6a1c33247f2f666f7c76de8d4fc2fd6f6124a Mon Sep 17 00:00:00 2001 From: mlemanczyk Date: Sat, 24 Jan 2026 01:28:09 +0100 Subject: [PATCH] BitContradiction performance improvements. Exact powers of 2 cases will reuse previous calculation results and transform them for the next q divisor. --- .../BitContradictionSolverWithAMultiplier.cs | 107 ++++++++++++++---- ...erDivisorByDivisorCpuTesterWithTemplate.cs | 64 +++++++++-- 2 files changed, 141 insertions(+), 30 deletions(-) diff --git a/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs b/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs index 58e10cf..9cda88d 100644 --- a/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs +++ b/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs @@ -14,6 +14,14 @@ internal sealed class PrivateWorkSet public ulong[] QMaskWords = []; public ulong[] AOneWin = []; public ulong[] AZeroWin = []; + public int[] PrevOffsets = []; + public int PrevOffsetsLength; + public bool HasPrevOffsets; + public BigInteger LastInvMod2k = BigInteger.Zero; + public bool HasLastInvMod2k; + public BigInteger LastPow2QMod2k = BigInteger.Zero; + public BigInteger LastPow2K = BigInteger.Zero; + public bool HasLastPow2QMod2k; [MethodImpl(MethodImplOptions.AggressiveInlining)] public void EnsureCapacity(int wordCount, int qWordCount) @@ -1288,12 +1296,16 @@ internal static class BitContradictionSolverWithAMultiplier internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, + BigInteger k, + bool isPowerOfTwo, out bool divides, out ContradictionReason reason) #else internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, + BigInteger k, + bool isPowerOfTwo, out bool divides) #endif { @@ -1397,37 +1409,84 @@ internal static class BitContradictionSolverWithAMultiplier BigInteger mask = Pow2Provider.BigIntegerMinusOne(ForcedBits); // Build q (low 1024 bits) from qOneOffsets. q is small anyway, but do it generically. - BigInteger qMod2k = BigInteger.Zero; + BigInteger qMod2k; + bool usedPow2Transform = false; + if (isPowerOfTwo && workSet.HasLastPow2QMod2k && workSet.LastPow2K != BigInteger.Zero && k == (workSet.LastPow2K << 1)) + { + qMod2k = ((workSet.LastPow2QMod2k << 1) - BigInteger.One) & mask; + usedPow2Transform = true; + } + else + { + qMod2k = BigInteger.Zero; + } + int i, t; for (i = 0; i < qOneOffsetsLength; i++) { t = qOneOffsets[i]; - qMod2k |= Pow2Provider.BigInteger(t); + if (!usedPow2Transform) + { + qMod2k |= Pow2Provider.BigInteger(t); + } qMaskWords[t >> 6] |= 1UL << (t & 63); } - // q must be odd for any Mersenne divisor candidate, and we always generate odd q candidates. No need to check it. - // if (qMod2k.IsEven) - // { - // Console.WriteLine("Pruned by qMod2k.IsEven"); - // #if DETAILED_LOG - // reason = ContradictionReason.ParityUnreachable; - // #endif - // return true; - // } - - // Compute inv = q^{-1} mod 2^1024 using Hensel/Newton iteration in Z/2^kZ. - BigInteger aLow = qMod2k; // OK seed for odd q. - - // maxKnownA stands for bits here. We're reusing the variable to limit registry pressure. - int maxKnownA; - for (maxKnownA = 1; maxKnownA < ForcedBits; maxKnownA <<= 1) + if (!usedPow2Transform) { - aLow = (aLow * (BigIntegerNumbers.Two - qMod2k * aLow)) & mask; + qMod2k &= mask; } + // Compute inv = q^{-1} mod 2^1024 using Hensel/Newton iteration in Z/2^kZ. + BigInteger invSeed; + int maxKnownA; + if (workSet.HasLastInvMod2k) + { + invSeed = workSet.LastInvMod2k; + for (int iter = 0; iter < 2; iter++) + { + invSeed = (invSeed * (BigIntegerNumbers.Two - qMod2k * invSeed)) & mask; + } + + if (((qMod2k * invSeed) & mask) != BigInteger.One) + { + invSeed = qMod2k; + for (maxKnownA = 1; maxKnownA < ForcedBits; maxKnownA <<= 1) + { + invSeed = (invSeed * (BigIntegerNumbers.Two - qMod2k * invSeed)) & mask; + } + } + } + else + { + invSeed = qMod2k; + for (maxKnownA = 1; maxKnownA < ForcedBits; maxKnownA <<= 1) + { + invSeed = (invSeed * (BigIntegerNumbers.Two - qMod2k * invSeed)) & mask; + } + } + + BigInteger aLow = (-invSeed) & mask; + + if (isPowerOfTwo) + { + workSet.LastPow2QMod2k = qMod2k; + workSet.LastPow2K = k; + workSet.HasLastPow2QMod2k = true; + } + + workSet.LastInvMod2k = invSeed; + workSet.HasLastInvMod2k = true; + if (workSet.PrevOffsets.Length < qOneOffsetsLength) + { + Array.Resize(ref workSet.PrevOffsets, qOneOffsetsLength); + } + qOneOffsets.CopyTo(workSet.PrevOffsets); + workSet.PrevOffsetsLength = qOneOffsetsLength; + workSet.HasPrevOffsets = true; + // a ≡ -aLow (mod 2^512) - aLow = (-aLow) & mask; + // Apply constraints for bits i in [0..min(511, maxAllowedA)] into knownOne0/knownZero0. // NOTE: now windowSize >= 512, so wordCount is enough to store these bits. @@ -1731,7 +1790,13 @@ internal static class BitContradictionSolverWithAMultiplier internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, - out bool divides) => TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, out divides, out _); + out bool divides) => TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, BigInteger.Zero, false, out divides, out _); +#endif +#if !DETAILED_LOG + internal static bool TryCheckDivisibilityFromOneOffsets( + ReadOnlySpan qOneOffsets, + ulong p, + out bool divides) => TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, BigInteger.Zero, false, out divides); #endif [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs index a5ba4ca..59fd888 100644 --- a/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs +++ b/PerfectNumbers.Core/Cpu/MersenneNumberDivisorByDivisorCpuTesterWithTemplate.cs @@ -1409,7 +1409,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN computedCycle = prime; // } - decided = TryExactBitTreeCheck(computedCycle, divisor, out divides); + decided = TryExactBitTreeCheck(computedCycle, divisor, currentK, out divides); if (!decided) { throw new InvalidOperationException($"BitTree check did not decide for p={prime}, divisor={divisor} (k={currentK})."); @@ -1437,7 +1437,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private bool TryExactBitTreeCheck(ulong exponent, BigInteger divisor, out bool divides) + private bool TryExactBitTreeCheck(ulong exponent, BigInteger divisor, BigInteger currentK, out bool divides) { #if DETAILED_LOG Stopwatch sw = Stopwatch.StartNew(); @@ -4617,7 +4617,7 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN sw.Reset(); sw.Start(); - bool decided = TryExactBitContradictionCheck(prime, divisor, out divides); + bool decided = TryExactBitContradictionCheck(prime, divisor, currentK, out divides); sw.Stop(); if (!decided) { @@ -4711,9 +4711,17 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN [ThreadStatic] private static int[]? _qBits; + [ThreadStatic] + private static int[]? _pow2Offsets; + [ThreadStatic] + private static int _pow2OffsetCount; + [ThreadStatic] + private static BigInteger _pow2K; + [ThreadStatic] + private static bool _hasPow2K; [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static bool TryExactBitContradictionCheck(ulong prime, BigInteger divisor, out bool divides) + private static bool TryExactBitContradictionCheck(ulong prime, BigInteger divisor, BigInteger currentK, out bool divides) { divides = false; // const ulong DebugTopDownDivisor = 1209708008767UL; @@ -4774,14 +4782,52 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN int[] oneOffsets = _qBits!; // Span oneOffsets = stackalloc int[qBitLength]; int oneCount = 0; - for (int i = 0; i < qBitLength; i++) + bool isPowerOfTwo = currentK.IsPowerOfTwo; + bool usedPow2Offsets = false; + int[]? pow2Offsets = _pow2Offsets; + int pow2OffsetCount = _pow2OffsetCount; + BigInteger pow2K = _pow2K; + bool hasPow2K = _hasPow2K; + if (isPowerOfTwo && hasPow2K && pow2Offsets != null && pow2OffsetCount > 0 && currentK == (pow2K << 1)) { - if (((divisor >> i) & BigInteger.One) == BigInteger.One) + int needed = pow2OffsetCount + 1; + if (needed <= oneOffsets.Length) { - oneOffsets[oneCount++] = i; + oneOffsets[0] = 0; + for (int i = 0; i < pow2OffsetCount; i++) + { + oneOffsets[i + 1] = pow2Offsets[i] + 1; + } + oneCount = needed; + usedPow2Offsets = true; } } + if (!usedPow2Offsets) + { + for (int i = 0; i < qBitLength; i++) + { + if (((divisor >> i) & BigInteger.One) == BigInteger.One) + { + oneOffsets[oneCount++] = i; + } + } + } + + if (isPowerOfTwo) + { + int[]? pow2OffsetsTarget = _pow2Offsets; + if (pow2OffsetsTarget == null || pow2OffsetsTarget.Length < oneCount) + { + pow2OffsetsTarget = new int[oneCount]; + _pow2Offsets = pow2OffsetsTarget; + } + Array.Copy(oneOffsets, pow2OffsetsTarget, oneCount); + _pow2OffsetCount = oneCount; + _pow2K = currentK; + _hasPow2K = true; + } + if (oneCount == 0) { divides = false; @@ -4791,9 +4837,9 @@ public struct MersenneNumberDivisorByDivisorCpuTesterWithTemplate() : IMersenneN ReadOnlySpan oneOffsetsSlice = oneOffsets[..oneCount]; // BitContradictionSolverWithQMultiplier.SetDebugEnabled(true); #if DETAILED_LOG - bool decided = BitContradictionSolverWithAMultiplier.TryCheckDivisibilityFromOneOffsets(oneOffsetsSlice, prime, out divides, out var reason); + bool decided = BitContradictionSolverWithAMultiplier.TryCheckDivisibilityFromOneOffsets(oneOffsetsSlice, prime, currentK, isPowerOfTwo, out divides, out var reason); #else - bool decided = BitContradictionSolverWithAMultiplier.TryCheckDivisibilityFromOneOffsets(oneOffsetsSlice, prime, out divides); + bool decided = BitContradictionSolverWithAMultiplier.TryCheckDivisibilityFromOneOffsets(oneOffsetsSlice, prime, currentK, isPowerOfTwo, out divides); #endif // BitContradictionSolverWithQMultiplier.SetDebugEnabled(false); if (decided && divides)