using System.Runtime.CompilerServices; using System.Numerics; using PerfectNumbers.Core; using ILGPU; namespace PerfectNumbers.Core.ByDivisor; internal sealed class PrivateWorkSet { // Active hot-path scratch buffer. public ulong[] KnownOne1 = []; public int[] CarryDropActiveAIndices = []; #if false private const int QPrefixReverseCacheSize = 1024; private const int ForcedALowBits = PerfectNumberConstants.BitContradictionForcedALowBits; // Currently disabled buffers retained for quick restore. public ulong[] KnownOne0 = []; public ulong[] KnownZero0 = []; public ulong[] KnownZero1 = []; public ulong[] QMaskWords = []; public ulong[] AOneWin = []; public ulong[] AZeroWin = []; public int[] PrevOffsets = []; public int PrevOffsetsLength; public bool HasPrevOffsets; public BigInteger[] InvByKMod = []; public bool[] HasInvByKMod = []; public bool HasLastInvMod2k; public BigInteger LastInvMod2k = BigInteger.Zero; public BigInteger LastPow2K = BigInteger.Zero; public BigInteger LastPow2QMod2k = BigInteger.Zero; public bool HasLastPow2QMod2k; public byte[] QMod2kBytes = []; public ulong[] QPrefixReverseKeys0 = []; public ulong[] QPrefixReverseKeys1 = []; public ulong[] QPrefixReverseValues0 = []; public ulong[] QPrefixReverseValues1 = []; public ulong[] QPrefixReverseValues2 = []; public ulong[] QPrefixReverseValues3 = []; public ulong[] QPrefixReverseValues4 = []; public ulong[] QPrefixReverseValues5 = []; public ulong[] QPrefixReverseValues6 = []; public ulong[] QPrefixReverseValues7 = []; public ulong[] QPrefixReverseValues8 = []; public ulong[] QPrefixReverseValues9 = []; public ulong[] QPrefixReverseValues10 = []; public ulong[] QPrefixReverseValues11 = []; public ulong[] QPrefixReverseValues12 = []; public ulong[] QPrefixReverseValues13 = []; public ulong[] QPrefixReverseValues14 = []; public ulong[] QPrefixReverseValues15 = []; public bool[] QPrefixReverseValid = []; public ulong[] LastQPrefix = []; public ulong[] LastQReverse = []; public bool HasLastQPrefix1024; #endif [MethodImpl(MethodImplOptions.AggressiveInlining)] public void EnsureCapacity(int wordCount, int qWordCount, int carryDropActiveCount) { if (KnownOne1.Length < wordCount) { EnsureUlong(ref KnownOne1, wordCount); } if (CarryDropActiveAIndices.Length < carryDropActiveCount) { EnsureInt(ref CarryDropActiveAIndices, carryDropActiveCount); } #if false if (KnownOne0.Length < wordCount) { EnsureUlong(ref KnownOne0, wordCount); EnsureUlong(ref KnownZero0, wordCount); EnsureUlong(ref KnownZero1, wordCount); } if (QMaskWords.Length < qWordCount) { EnsureUlong(ref QMaskWords, qWordCount); EnsureUlong(ref AOneWin, qWordCount); EnsureUlong(ref AZeroWin, qWordCount); } if (InvByKMod.Length != ForcedALowBits) { InvByKMod = new BigInteger[ForcedALowBits]; HasInvByKMod = new bool[ForcedALowBits]; } if (QMod2kBytes.Length != (ForcedALowBits >> 3)) { QMod2kBytes = new byte[ForcedALowBits >> 3]; } if (QPrefixReverseKeys0.Length != QPrefixReverseCacheSize) { QPrefixReverseKeys0 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseKeys1 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues0 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues1 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues2 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues3 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues4 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues5 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues6 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues7 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues8 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues9 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues10 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues11 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues12 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues13 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues14 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValues15 = new ulong[QPrefixReverseCacheSize]; QPrefixReverseValid = new bool[QPrefixReverseCacheSize]; } if (LastQPrefix.Length != 16) { LastQPrefix = new ulong[16]; LastQReverse = new ulong[16]; } #endif } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void EnsureUlong(ref ulong[] arr, int needed) => Array.Resize(ref arr, needed); [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void EnsureInt(ref int[] arr, int needed) => Array.Resize(ref arr, needed); } /// /// Helper utilities for the BitContradiction divisor search. These methods implement the /// interval-based carry propagation described in the BitContradictions plan document. /// The solver itself will build on top of these primitives. /// [DeviceDependentTemplate(typeof(ComputationDevice))] internal static class BitContradictionSolverWithAMultiplier { #if DETAILED_LOG private const ulong DebugTargetP = 138000001UL; private static readonly BigInteger DebugTargetK = new(3); private static readonly ulong DebugTargetQ = checked((ulong)((DebugTargetK * DebugTargetP) << 1)); // private const ulong DebugTargetQ = 1209708008767UL; [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool IsTargetCase(ulong p, System.Numerics.BigInteger k) => p == DebugTargetP && k == DebugTargetK; [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void Dbg(ulong p, System.Numerics.BigInteger k, string msg) { if (IsTargetCase(p, k)) Console.WriteLine(msg); } #endif #if WITH_PRUNE_STATS private static volatile int _pruneStatsConsoleCounter; private static ulong _prunedWithGetAKnownStateRowAware; private static ulong _prunedWithTryRunHighBitAndBorrowPrefiltersCombined; private static ulong _prunedWithTargetJ; private static ulong _prunedWithInterMinInterMax; private static ulong _prunedWithTryProcessBottomUpBlockRowAware; private static ulong _prunedWithTryAdvanceStableUnknown; private static ulong _prunedWithBestCaseCarryOutOneStep; private static ulong _prunedWithCouldStillReachTopCarryZero_UpperBound; private static ulong _prunedWithTryFindSingleUnknownInQ_Prefix; private static ulong _prunedWithTryPropagateCarry1; private static ulong _prunedWithTryPropagateCarry2; private static ulong _prunedWithTopCarry; private static ulong _prunedWithPowmod2FinalShift; private static ulong _prunedWithTruncatedLength1; private static ulong _prunedWithTruncatedLength2; private static ulong _prunedWithTryHighBitCarryRangeToColumnJ; private static ulong _prunedWithCarryLowNonZeroFinal; private static ulong _prunedWithEvenQAfterBuild; private static ulong _prunedWithPowmod2FinalResidue; private static ulong _prunedWithPowmod2Direct; private static ulong _prunedWithPowmod2Tracked; private static ulong _prunedWithPowmod2TrackedOptimized; private static ulong _prunedWithPowmod2ZeroBeforeLastBlock; private static ulong _prunedWithTryDynamicModularRangePrune; private static ulong _prunedWithTopPrefixRangeConflict; private static ulong _prunedWithQTimesAInterval; private static ulong _prunedWithTopStateConflict; private static ulong _prunedWithTopBothZeroConflict; private static ulong _prunedWithCarryParityPrecheck; private static ulong _prunedWithKnownBitsModResidue; private static ulong _prunedWithCarryDropHeuristic; private static ulong _prunedWithCarryDominance67; private static ulong _prunedWithCarryDominance8910; private static ulong _prunedWithCarryDominance15161718; private static ulong _prunedWithCarryDominanceReached10AndStayedAbove5; private static ulong _prunedWithCarryDominance1011121314; private static ulong _prunedWithHighDensityWithShortLowStreaksCarry; #endif private const int PingPongBlockColumns = 4096; private const int TailCarryBatchColumns = 16384; private const int TailCarryBatchRepeats = 1; private const int TopDownBorrowColumns = 4096; private const int TopDownBorrowMinUnknown = 12; private const int PrefixRunCacheOffsets = 256; private const int PrefixRunCacheCapacity = 256; private const int ModLocalPrefilterMaxABits = 20; private const int TopDownCacheCapacity = 16; private const int HighBitCarryPrefilterColumns = 16; private const int MiniPingPongPrefilterColumns = 60; private const int DynamicModPeriodicColumns = 1024; private const int DynamicModPeriodicMask = DynamicModPeriodicColumns - 1; private const int TopPrefixProbeBits = 24; private const int ForcedALowBits = PerfectNumberConstants.BitContradictionForcedALowBits; private const int QPrefixReverseCacheSize = 1024; private const int QPrefixReverseCacheMask = QPrefixReverseCacheSize - 1; private static readonly int[] ModLocalPrefilterModuli = { 65537, 65521 }; private static readonly int[] KnownBitsModResidueModuli = { 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 }; private const int CarryDropHeuristicBlockBits = 64; private const int Powmod2AdaptiveTargetBlocks = 512; private const int Powmod2AdaptiveMinRatio = 16; private const int CarryDropHeuristicMaxBits = 1024; private const bool EnableCarryDropHeuristic = true; private static readonly ulong[] BitMask64 = CreateBitMask64(); private static readonly byte[] ReverseByteLut = CreateReverseByteLut(); #if WITH_PRUNE_STATS [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void IncrementPruneCounter(ref ulong counter) { Atomic.Add(ref counter, 1UL); _pruneStatsConsoleCounter++; if (_pruneStatsConsoleCounter < PerfectNumberConstants.ConsoleInterval) { return; } PrintPruneStats(); _pruneStatsConsoleCounter = 0; } private static void PrintPruneStats() { Console.WriteLine("Prune stats:"); Console.WriteLine( $"TryPropagateCarry1={_prunedWithTryPropagateCarry1}, " + $"Powmod2ZeroBeforeLastBlock={_prunedWithPowmod2ZeroBeforeLastBlock}, " + $"Powmod2Tracked={_prunedWithPowmod2Tracked}, " + $"Powmod2TrackedOptimized={_prunedWithPowmod2TrackedOptimized}, " + $"CarryDropHeuristic={_prunedWithCarryDropHeuristic}, " + $"CarryDominance67={_prunedWithCarryDominance67}, " + $"CarryDominance8910AndMore={_prunedWithCarryDominance8910}, " + $"CarryDominance15161718={_prunedWithCarryDominance15161718}, " + $"CarryDominanceReached10AndStayedAbove5={_prunedWithCarryDominanceReached10AndStayedAbove5}, " + $"CarryDominance1011121314={_prunedWithCarryDominance1011121314}, " + $"HighDensityWithShortLowStreaksCarry={_prunedWithHighDensityWithShortLowStreaksCarry}" ); } #endif [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ValidateHeuristicPruneAgainstPowmod2OrThrow(BigInteger q, int p, string filterName) { if (BigInteger.ModPow(BigIntegerNumbers.Two, p, q) == BigInteger.One) { throw new InvalidOperationException($"Heuristic prune '{filterName}' rejected a valid divisor for p={p}, q={q}."); } } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static long BestCaseCarryOutOneStep(long carryIn, int forced, int possible, int requiredBit /*1*/) { // Choose sum in [forced..possible] s.t. (carryIn + sum) % 2 == requiredBit, // minimize carryOut = (carryIn + sum - requiredBit) / 2. int needSumParity = (requiredBit ^ ((int)carryIn & 1)) & 1; int sum = forced; if ((sum & 1) != needSumParity) sum++; if (sum > possible) return -1; // parity infeasible in this column return (carryIn + sum - requiredBit) >> 1; } // --- ADD near other helpers in BitContradictionSolverWithAMultiplier.cs --- [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool CouldStillReachTopCarryZero_UpperBound( CarryRange carryIn, int currentColumn, int highColumn, int maxPossibleOnesPerColumn) // this is "possible" = forced+remaining for *this* column; { // Quick conservative check: // If even the *maximum* per-column ones cannot prevent carry from staying >=1 by the time we hit highColumn, // we can prune. We use a crude but safe bound: // // carry_{t+1} >= floor(carry_t / 2) // and carry_{t+1} <= floor((carry_t + maxPossibleOnesPerColumn) / 2) // // For a SAFE prune we need a LOWER bound on the best case (smallest achievable carry). // The smallest carry you can get after one step is at least floor(carryMin / 2). // After L steps: at least floor(carryMin / 2^L). // // If that lower bound is still >=1 at the top, impossible. int stepsLeft = highColumn - currentColumn + 1; if (stepsLeft <= 0) return true; long minCarry = carryIn.Min; if (minCarry <= 0) return true; // floor(minCarry / 2^stepsLeft) if (stepsLeft >= 63) return false; // 2^63 exceeds long shifts; if minCarry>0 it will go to 0, so be conservative long minPossibleAtTop = minCarry >> stepsLeft; return minPossibleAtTop == 0; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong[] CreateBitMask64() { ulong[] masks = new ulong[64]; ulong value = 1UL; for (int i = 0; i < 64; i++) { masks[i] = value; value <<= 1; } return masks; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static byte[] CreateReverseByteLut() { byte[] lut = new byte[256]; for (int i = 0; i < lut.Length; i++) { int v = i; int r = 0; for (int b = 0; b < 8; b++) { r = (r << 1) | (v & 1); v >>= 1; } lut[i] = (byte)r; } return lut; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool IntervalHasParity(long lo, long hi, int parity) // parity: 0=even,1=odd { if (lo > hi) return false; return ((lo & 1L) == parity) || (lo + 1 <= hi); } #if false [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong ReverseBits64(ulong value) { return ((ulong)ReverseByteLut[(byte)value] << 56) | ((ulong)ReverseByteLut[(byte)(value >> 8)] << 48) | ((ulong)ReverseByteLut[(byte)(value >> 16)] << 40) | ((ulong)ReverseByteLut[(byte)(value >> 24)] << 32) | ((ulong)ReverseByteLut[(byte)(value >> 32)] << 24) | ((ulong)ReverseByteLut[(byte)(value >> 40)] << 16) | ((ulong)ReverseByteLut[(byte)(value >> 48)] << 8) | ReverseByteLut[(byte)(value >> 56)]; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ReverseBits1024(in ReadOnlySpan prefix, Span reversed) { reversed[0] = ReverseBits64(prefix[15]); reversed[1] = ReverseBits64(prefix[14]); reversed[2] = ReverseBits64(prefix[13]); reversed[3] = ReverseBits64(prefix[12]); reversed[4] = ReverseBits64(prefix[11]); reversed[5] = ReverseBits64(prefix[10]); reversed[6] = ReverseBits64(prefix[9]); reversed[7] = ReverseBits64(prefix[8]); reversed[8] = ReverseBits64(prefix[7]); reversed[9] = ReverseBits64(prefix[6]); reversed[10] = ReverseBits64(prefix[5]); reversed[11] = ReverseBits64(prefix[4]); reversed[12] = ReverseBits64(prefix[3]); reversed[13] = ReverseBits64(prefix[2]); reversed[14] = ReverseBits64(prefix[1]); reversed[15] = ReverseBits64(prefix[0]); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void GetReversePrefix1024(PrivateWorkSet workSet, in ReadOnlySpan prefix, Span reversed) { ulong p0 = prefix[0]; ulong p1 = prefix[1]; ulong[] keys0 = workSet.QPrefixReverseKeys0; ulong[] keys1 = workSet.QPrefixReverseKeys1; ulong[] values0 = workSet.QPrefixReverseValues0; ulong[] values1 = workSet.QPrefixReverseValues1; ulong[] values2 = workSet.QPrefixReverseValues2; ulong[] values3 = workSet.QPrefixReverseValues3; ulong[] values4 = workSet.QPrefixReverseValues4; ulong[] values5 = workSet.QPrefixReverseValues5; ulong[] values6 = workSet.QPrefixReverseValues6; ulong[] values7 = workSet.QPrefixReverseValues7; ulong[] values8 = workSet.QPrefixReverseValues8; ulong[] values9 = workSet.QPrefixReverseValues9; ulong[] values10 = workSet.QPrefixReverseValues10; ulong[] values11 = workSet.QPrefixReverseValues11; ulong[] values12 = workSet.QPrefixReverseValues12; ulong[] values13 = workSet.QPrefixReverseValues13; ulong[] values14 = workSet.QPrefixReverseValues14; ulong[] values15 = workSet.QPrefixReverseValues15; bool[] valid = workSet.QPrefixReverseValid; int slot = (int)p0 & QPrefixReverseCacheMask; if (valid[slot] && keys0[slot] == p0 && keys1[slot] == p1) { reversed[0] = values0[slot]; reversed[1] = values1[slot]; reversed[2] = values2[slot]; reversed[3] = values3[slot]; reversed[4] = values4[slot]; reversed[5] = values5[slot]; reversed[6] = values6[slot]; reversed[7] = values7[slot]; reversed[8] = values8[slot]; reversed[9] = values9[slot]; reversed[10] = values10[slot]; reversed[11] = values11[slot]; reversed[12] = values12[slot]; reversed[13] = values13[slot]; reversed[14] = values14[slot]; reversed[15] = values15[slot]; return; } ReverseBits1024(prefix, reversed); keys0[slot] = p0; keys1[slot] = p1; values0[slot] = reversed[0]; values1[slot] = reversed[1]; values2[slot] = reversed[2]; values3[slot] = reversed[3]; values4[slot] = reversed[4]; values5[slot] = reversed[5]; values6[slot] = reversed[6]; values7[slot] = reversed[7]; values8[slot] = reversed[8]; values9[slot] = reversed[9]; values10[slot] = reversed[10]; values11[slot] = reversed[11]; values12[slot] = reversed[12]; values13[slot] = reversed[13]; values14[slot] = reversed[14]; values15[slot] = reversed[15]; valid[slot] = true; } #endif #if DETAILED_LOG [ThreadStatic] private static long _statStableAdvanceCalls; [ThreadStatic] private static long _statStableAdvanceCols; [ThreadStatic] private static long _statStableAdvanceColsBatched; [ThreadStatic] private static long _statScalarCols; [ThreadStatic] private static long _statDeltaRequests; [ThreadStatic] private static long _statDeltaHits; [ThreadStatic] private static long _statDeltaBuilds; [ThreadStatic] private static long _statStableAdvanceBigShift; // how many times columnCount>=64 in TryAdvanceStableUnknown [ThreadStatic] private static long _statBatchEntry; [ThreadStatic] private static long _statBatchCols; [ThreadStatic] private static long _statCleanChecks; [ThreadStatic] private static long _statCleanFalse; [ThreadStatic] private static int _statFirstDirtyColumn; [ThreadStatic] private static int _statFirstDirtyWord; [ThreadStatic] private static ulong _statFirstDirtyQ; [ThreadStatic] private static ulong _statFirstDirtyKnown; public readonly struct TopDownPruneFailure(long column, long carryMin, long carryMax, long unknown) { public readonly long Column = column; public readonly long CarryMin = carryMin; public readonly long CarryMax = carryMax; public readonly long Unknown = unknown; } private readonly struct TopDownCacheValue(bool success, TopDownPruneFailure? failure) { public readonly bool Success = success; public readonly TopDownPruneFailure? Failure = failure; } [ThreadStatic] private static Dictionary? _topDownCache; [ThreadStatic] private static TopDownPruneFailure? _lastTopDownFailure; internal static TopDownPruneFailure? LastTopDownFailure => _lastTopDownFailure; private static void PrintStats() { Console.WriteLine( $"[STATS] stableCalls={_statStableAdvanceCalls} stableCols={_statStableAdvanceColsBatched} " + $"scalarCols={_statScalarCols} bigShift={_statStableAdvanceBigShift} " + $"deltaReq={_statDeltaRequests} deltaHits={_statDeltaHits} deltaBuilds={_statDeltaBuilds}"); Console.WriteLine( $"[CLEAN] checks={_statCleanChecks} false={_statCleanFalse} " + $"firstDirtyCol={_statFirstDirtyColumn} word={_statFirstDirtyWord} " + $"qWord=0x{_statFirstDirtyQ:X16} knownWord=0x{_statFirstDirtyKnown:X16}"); } #else [ThreadStatic] private static Dictionary? _topDownCache; #endif [ThreadStatic] private static long _dynamicModChecks; [ThreadStatic] private static long _dynamicModFailures; [ThreadStatic] private static Dictionary? _qModCache; [ThreadStatic] private static List? _topDownCacheOrder; [ThreadStatic] private static Dictionary? _prefixRunCache; [ThreadStatic] private static List? _prefixRunCacheOrder; [ThreadStatic] private static Dictionary? _pow2DiffModCache; [ThreadStatic] private static PrivateWorkSet _privateWorkSet; private readonly struct QModCacheKey(long hash, int mod) { public readonly long Hash = hash; public readonly int Mod = mod; public bool Equals(QModCacheKey other) => Hash == other.Hash && Mod == other.Mod; public override bool Equals(object? obj) => obj is QModCacheKey other && Equals(other); public override int GetHashCode() => HashCode.Combine(Hash, Mod); } private readonly struct TopDownCacheKey(long hash, int length, long maxAllowedA, long pLong) { public readonly long Hash = hash; public readonly int Length = length; public readonly long MaxAllowedA = maxAllowedA; public readonly long PLong = pLong; public bool Equals(TopDownCacheKey other) => Hash == other.Hash && Length == other.Length && MaxAllowedA == other.MaxAllowedA && PLong == other.PLong; public override bool Equals(object? obj) => obj is TopDownCacheKey other && Equals(other); public override int GetHashCode() => HashCode.Combine(Hash, Length, MaxAllowedA, PLong); } private readonly struct PrefixRunCacheEntry( int prefixCount, int runCount, int lastOffset, int lastRunStart, int lastRunLength, int[] runStarts, int[] runLengths) { public readonly int PrefixCount = prefixCount; public readonly int RunCount = runCount; public readonly int LastOffset = lastOffset; public readonly int LastRunStart = lastRunStart; public readonly int LastRunLength = lastRunLength; public readonly int[] RunStarts = runStarts; public readonly int[] RunLengths = runLengths; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static byte AllowedCarryInMask4FromCarryOut(byte carryOutMask4, byte onesMask4) { byte m = 0; for (int co = 0; co < 4; co++) { if (((carryOutMask4 >> co) & 1) == 0) continue; int r = ((co << 1) + 1) & 3; // (2*carryOut+1) mod 4 for (int o = 0; o < 4; o++) { if (((onesMask4 >> o) & 1) == 0) continue; int ci = (r - o) & 3; // parzystość: carryIn + ones musi być nieparzyste if (((ci + o) & 1) == 0) continue; m |= (byte)(1 << ci); } } return m; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static byte NextCarryMask4(byte carryInMask4, byte onesMask4, int requiredResultBit) { int b = requiredResultBit & 1; byte outMask = 0; for (int ci = 0; ci < 4; ci++) { if (((carryInMask4 >> ci) & 1) == 0) continue; for (int o = 0; o < 4; o++) { if (((onesMask4 >> o) & 1) == 0) continue; // warunek bitu: (carryIn + ones) ≡ b (mod 2) if (((ci + o) & 1) != b) continue; // s = carryIn + ones - b jest parzyste; bazujemy na resztach mod 4, // ale UWAGA: (carryIn + 4) zmienia (s/2) o +2 => trzeba domknąć o +2. int s = (ci + o - b) & 3; // s mod 4, zawsze 0 albo 2 int baseCo = ((s >> 1) & 3); // 0 albo 1 outMask |= (byte)(1 << baseCo); outMask |= (byte)(1 << (baseCo ^ 2)); // +2 (mod4) dla bezpieczeństwa } } return outMask; } private static Vector SimdPopCount(Vector value) { // A classic parallel popcount algorithm (SWAR) var m1 = new Vector(0x5555555555555555); // 0101... var m2 = new Vector(0x3333333333333333); // 0011... var m4 = new Vector(0x0F0F0F0F0F0F0F0F); // 00001111... var h01 = new Vector(0x0101010101010101); value -= (value >> 1) & m1; value = (value & m2) + ((value >> 2) & m2); value = (value + (value >> 4)) & m4; return (value * h01) >> 56; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ComputeForcedRemainingAndClean( in ulong[] qMask, in ulong[] oneWin, in ulong[] zeroWin, int wordCount, out int forced, out int remaining, out bool clean) { // Use scalar implementation for small inputs or when SIMD is not available. int i; if (!Vector.IsHardwareAccelerated || wordCount < Vector.Count) { forced = 0; remaining = 0; clean = true; for (i = 0; i < wordCount; i++) { ulong q = qMask[i]; ulong one = oneWin[i]; ulong known = one | zeroWin[i]; one &= q; ulong unknownBits = q & ~known; forced += PopCount(one); remaining += PopCount(unknownBits); if ((q & known) != 0) { clean = false; } } return; } // Vectorized implementation var vForced = Vector.Zero; var vRemaining = Vector.Zero; var vContradictions = Vector.Zero; int limit = wordCount - (wordCount % Vector.Count); for (i = 0; i < limit; i += Vector.Count) { var vq = new Vector(qMask, i); var vOne = new Vector(oneWin, i); var vZero = new Vector(zeroWin, i); var vKnown = vOne | vZero; var vForcedBits = vOne & vq; var vUnknownBits = vq & ~vKnown; vForced += SimdPopCount(vForcedBits); vRemaining += SimdPopCount(vUnknownBits); vContradictions |= (vq & vKnown); } forced = (int)Vector.Sum(vForced); remaining = (int)Vector.Sum(vRemaining); clean = (vContradictions == Vector.Zero); // Process the remainder for (; i < wordCount; i++) { ulong q = qMask[i]; ulong one = oneWin[i]; ulong known = one | zeroWin[i]; one &= q; ulong unknownBits = q & ~known; forced += PopCount(one); remaining += PopCount(unknownBits); if (clean && (q & known) != 0) { clean = false; } } } private static readonly byte[] _popCount16Cache = CreatePopCount16Cache(); private static byte[] CreatePopCount16Cache() { const int CacheSize = 65536; byte[] cache = new byte[CacheSize]; for (int i = 0; i < CacheSize; i++) { cache[i] = (byte)BitOperations.PopCount((uint)i); } return cache; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int PopCount(ulong value) { byte[] cache = _popCount16Cache; int result = cache[(ushort)value]; value >>= 16; result += cache[(ushort)value]; value >>= 16; result += cache[(ushort)value]; value >>= 16; result += cache[(ushort)value]; return result; } private static long ComputeQOffsetsHash(in ReadOnlySpan qOneOffsets) { const long fnvOffset = unchecked(1469598103934665603); const long fnvPrime = 1099511628211; long hash = fnvOffset; int qOneOffsetsLength = qOneOffsets.Length; for (int i = 0; i < qOneOffsetsLength; i++) { hash ^= qOneOffsets[i]; hash *= fnvPrime; } return hash; } private static int GetQModCached(ReadOnlySpan qOneOffsets, long qHash, int mod) { Dictionary? cache = _qModCache; if (cache is null) { cache = []; _qModCache = cache; } var key = new QModCacheKey(qHash, mod); if (cache.TryGetValue(key, out int cached)) { return cached; } int value = ComputeQMod(qOneOffsets, mod); cache[key] = value; return value; } private static long ComputeQOffsetsPrefixHash(in ReadOnlySpan qOneOffsets, int count) { const long fnvOffset = unchecked(1469598103934665603); const long fnvPrime = 1099511628211; long hash = fnvOffset; for (int i = 0; i < count; i++) { hash ^= qOneOffsets[i]; hash *= fnvPrime; } hash ^= count; hash *= fnvPrime; return hash; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool IsCleanInQ(in ulong[] qMask, in ulong[] oneWin, in ulong[] zeroWin, int wordCount) { for (int i = 0; i < wordCount; i++) { if ((qMask[i] & (oneWin[i] | zeroWin[i])) != 0) return false; } return true; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static byte MaskFromIntervalMod4(long lo, long hi) { if (lo > hi) return 0; // Jeżeli długość >= 4, to mamy wszystkie klasy reszt. if (hi - lo >= 3) return 0b1111; byte m = 0; for (long x = lo; x <= hi; x++) m |= (byte)(1 << (int)(x & 3)); return m; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static long AlignUpToAllowedMod4(long value, byte mask) { // przesuwa do najbliższej wartości >=value której (value mod 4) jest w masce for (int d = 0; d < 4; d++) { long v = value + d; if (((mask >> (int)(v & 3)) & 1) != 0) return v; } return long.MaxValue; // brak } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static long AlignDownToAllowedMod4(long value, byte mask) { for (int d = 0; d < 4; d++) { long v = value - d; if (((mask >> (int)(v & 3)) & 1) != 0) return v; } return long.MinValue; } private static bool TryGetPrefixRunCache(long key, out PrefixRunCacheEntry entry) { var cache = _prefixRunCache; if (cache is null) { cache = []; _prefixRunCache = cache; } return cache.TryGetValue(key, out entry); } private static void StorePrefixRunCache(long key, in PrefixRunCacheEntry entry) { _prefixRunCacheOrder ??= new List(PrefixRunCacheCapacity + 1); if (_prefixRunCache!.ContainsKey(key)) { _prefixRunCache[key] = entry; return; } if (_prefixRunCacheOrder.Count >= PrefixRunCacheCapacity) { long evict = _prefixRunCacheOrder[0]; _prefixRunCacheOrder.RemoveAt(0); _prefixRunCache.Remove(evict); } _prefixRunCache[key] = entry; _prefixRunCacheOrder.Add(key); } private static PrefixRunCacheEntry BuildPrefixRunCache(in ReadOnlySpan qOneOffsets, int prefixCount) { // TODO: Remove per-call allocations int[] starts = new int[prefixCount]; int[] lengths = new int[prefixCount]; int runCount = 0; int prev = qOneOffsets[0]; int currentStart = prev; int currentLength = 1; for (int i = 1; i < prefixCount; i++) { int offset = qOneOffsets[i]; if (offset == prev + 1) { currentLength++; } else { starts[runCount] = currentStart; lengths[runCount] = currentLength; runCount++; currentStart = offset; currentLength = 1; } prev = offset; } starts[runCount] = currentStart; lengths[runCount] = currentLength; runCount++; int[] cachedStarts = new int[runCount]; int[] cachedLengths = new int[runCount]; Array.Copy(starts, cachedStarts, runCount); Array.Copy(lengths, cachedLengths, runCount); return new PrefixRunCacheEntry( prefixCount, runCount, prev, currentStart, currentLength, cachedStarts, cachedLengths); } #if DETAILED_LOG private static bool TryGetTopDownCache(in TopDownCacheKey key, out TopDownCacheValue value) #else private static bool TryGetTopDownCache(in TopDownCacheKey key, out bool value) #endif { // TODO: Remove this check after Initialize if (_topDownCache == null) { value = default; return false; } return _topDownCache.TryGetValue(key, out value); } #if DETAILED_LOG private static void StoreTopDownCache(in TopDownCacheKey key, in TopDownCacheValue value) #else private static void StoreTopDownCache(in TopDownCacheKey key, in bool value) #endif { // TODO: Remove this check after Initialize _topDownCache ??= []; _topDownCacheOrder ??= new List(TopDownCacheCapacity + 1); if (_topDownCache.ContainsKey(key)) { _topDownCache[key] = value; return; } if (_topDownCacheOrder.Count >= TopDownCacheCapacity) { var evict = _topDownCacheOrder[0]; _topDownCacheOrder.RemoveAt(0); _topDownCache.Remove(evict); } _topDownCache[key] = value; _topDownCacheOrder.Add(key); } internal static ColumnBounds ComputeColumnBounds( ReadOnlySpan multiplicand, ReadOnlySpan multiplier, int columnIndex) { int forced = 0; int possible = 0; int multiplicandLength = multiplicand.Length; int multiplierLength = multiplier.Length; for (int i = 0; i <= columnIndex; i++) { bool a = i < multiplicandLength && multiplicand[i]; bool q = columnIndex - i < multiplierLength && multiplier[columnIndex - i]; if (!a && !q) { continue; } if (a && q) { forced++; } possible++; } return new ColumnBounds(forced, possible); } internal static bool TryPropagateCarry( ref CarryRange currentCarry, int forcedOnes, int possibleOnes, int requiredResultBit) { // sum = carryIn + ones long minSum = currentCarry.Min + forcedOnes; long maxSum = currentCarry.Max + possibleOnes; int parity = requiredResultBit & 1; // enforce sum parity: sum ≡ requiredResultBit (mod 2) minSum = AlignUpToParity(minSum, parity); maxSum = AlignDownToParity(maxSum, parity); if (minSum > maxSum) return false; // carryOut = (sum - requiredBit) / 2 long nextMin = (minSum - requiredResultBit) >> 1; long nextMax = (maxSum - requiredResultBit) >> 1; // // ones-range is [forcedOnes..possibleOnes], not [nextMin..nextMax] // byte onesMask4 = MaskFromIntervalMod4(forcedOnes, possibleOnes); // // IMPORTANT: NextCarryMask4 must account for requiredResultBit, otherwise it may over-prune for b=0. // byte nextMask4 = NextCarryMask4(currentCarry.Mod4Mask, onesMask4, requiredResultBit); // if (nextMask4 == 0) // return false; // // Optional safety: intersect with the interval’s own mod4 availability // byte intervalMask4 = MaskFromIntervalMod4(nextMin, nextMax); // nextMask4 &= intervalMask4; // if (nextMask4 == 0) // return false; currentCarry = new CarryRange(nextMin, nextMax, byte.MaxValue); return true; } /// /// Mini ping-pong (Option B): compute a conservative allowed carry range after walking /// top-down for a given number of columns starting at . /// /// Semantics: /// - We treat "carryOut at column = startColumn+1" as the starting set (seed). /// - Each step processes one column (carryOut -> carryIn) moving downward. /// - After columns, the resulting set corresponds to carry-in /// at column (startColumn - steps + 1). /// private static bool TryHighBitCarryRangeToColumnJ( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int pLong, int maxAllowedA, int startColumn, // e.g. p-1, then J1-1, ... int steps, // number of columns to process downward long startCarryOutMin, // seed carryOut range long startCarryOutMax, out long allowedMin, out long allowedMax) { allowedMin = 0; allowedMax = long.MaxValue; if (steps <= 0) { // No processing: carry-in == seed (conservative) allowedMin = startCarryOutMin < 0 ? 0 : startCarryOutMin; allowedMax = startCarryOutMax < 0 ? 0 : startCarryOutMax; if (allowedMax < allowedMin) (allowedMin, allowedMax) = (0, long.MaxValue); return true; } // Clamp startColumn to legal high column range. If outside, be conservative. int maxHigh = pLong - 1; if (startColumn > maxHigh || startColumn < 0) { allowedMin = 0; allowedMax = long.MaxValue; return true; } // Determine last processed column. int endColumn = startColumn - (steps - 1); if (endColumn < 0) { // Processing would go past column 0. We can clamp steps to reach column 0 exactly. steps = startColumn + 1; endColumn = 0; if (steps <= 0) { allowedMin = 0; allowedMax = 0; return true; } } // Two fixed interval buffers on stack: each interval is [lo, hi], disjoint and sorted. Span lo0 = stackalloc long[64]; Span hi0 = stackalloc long[64]; Span lo1 = stackalloc long[64]; Span hi1 = stackalloc long[64]; Span outLo = lo0, outHi = hi0; Span inLo = lo1, inHi = hi1; int outCount = 1; outLo[0] = startCarryOutMin; outHi[0] = startCarryOutMax; int inCount = 0; bool flip = false; [MethodImpl(MethodImplOptions.AggressiveInlining)] static bool AddInterval(ref int count, Span loArr, Span hiArr, long lo, long hi) { if (hi < lo) return true; if (lo < 0) lo = 0; if (hi < 0) return true; if (count == 0) { loArr[0] = lo; hiArr[0] = hi; count = 1; return true; } long lastHi = hiArr[count - 1]; if (lo <= lastHi + 1) { if (hi > lastHi) hiArr[count - 1] = hi; return true; } if ((uint)count >= (uint)loArr.Length) { // Too many intervals -> disable this filter conservatively. count = 0; // signal overflow / disabled return true; // do NOT fail the candidate } loArr[count] = lo; hiArr[count] = hi; count++; return true; } for (int step = 0; step < steps; step++) { int column = startColumn - step; // n = count of q-offsets t such that 0 <= column - t <= maxAllowedA // <=> column - maxAllowedA <= t <= column int lowerT; long lowerTLong = (long)column - (long)maxAllowedA; if (lowerTLong <= 0) lowerT = 0; else if (lowerTLong >= int.MaxValue) lowerT = int.MaxValue; else lowerT = (int)lowerTLong; int upperT = column <= 0 ? 0 : column; FindBounds(qOneOffsets, lowerT, upperT, out int left, out int right); int n = right - left; if (n < 0) return false; inCount = 0; for (int i = 0; i < outCount; i++) { long a = outLo[i]; long b = outHi[i]; if (b < 0) continue; if (a < 0) a = 0; if (b < a) continue; const long TwicePlusOneOverflowLimit = (long.MaxValue - 1) >> 1; // If carryOut is too large for safe 2*x+1 in long, we must remain conservative. if (a > TwicePlusOneOverflowLimit || b > TwicePlusOneOverflowLimit) { allowedMin = 0; allowedMax = long.MaxValue; return true; } long rMin = (a << 1) + 1; long rMax = (b << 1) + 1; // carryIn ∈ [ max(0, r-n), r ] if (n <= 0) { // n==0 => ones==0 deterministycznie, bit wyniku = 1 => carryIn musi być nieparzyste. long lo02 = rMin; long hi02 = rMax; if (!IntervalHasParity(lo02, hi02, 1)) continue; if ((lo02 & 1L) == 0) lo02++; if (!AddInterval(ref inCount, inLo, inHi, lo02, hi02) || inCount == 0) { allowedMin = 0; allowedMax = long.MaxValue; return true; } continue; } long lo = rMin - n; if (lo < 0) lo = 0; if (!AddInterval(ref inCount, inLo, inHi, lo, rMax) || inCount == 0) { allowedMin = 0; allowedMax = long.MaxValue; return true; } } if (inCount <= 0) { return false; } // Swap buffers for next step. outCount = inCount; flip = !flip; if (!flip) { outLo = lo0; outHi = hi0; inLo = lo1; inHi = hi1; } else { outLo = lo1; outHi = hi1; inLo = lo0; inHi = hi0; } } // After steps, 'out' corresponds to carry-in at the last processed column (endColumn). allowedMin = outLo[0]; allowedMax = outHi[outCount - 1]; #if DETAILED_LOG Console.WriteLine($"[TD] steps={steps} endRange=[{allowedMin},{allowedMax}] intervals={outCount}"); #endif return true; } private static bool TryModularPrefilter( ReadOnlySpan qOneOffsets, long qHash, ulong p, long maxAllowedA) { if (maxAllowedA < 0 || maxAllowedA > ModLocalPrefilterMaxABits) { return true; } int bitCount = (int)(maxAllowedA + 1); if (bitCount <= 0 || bitCount >= 62) { return true; } long maxAValue = (1L << bitCount) - 1L; int modLocalPrefilterModuliLength = ModLocalPrefilterModuli.Length; for (int i = 0; i < modLocalPrefilterModuliLength; i++) { int mod = ModLocalPrefilterModuli[i]; int qMod = GetQModCached(qOneOffsets, qHash, mod); int target = ModPow2Minus1(p, mod); if (qMod == 0) { if (target != 0) { return false; } continue; } int inverse = ModPow(qMod, mod - 2, mod); long x = (long)target * inverse % mod; if (x > maxAValue) { return false; } } return true; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void FindBounds(ReadOnlySpan values, int lowerTarget, int upperTarget, out int lower, out int upper) { // Find lower bound for lowerTarget int lo = 0; int valuesLength = values.Length; int hi = valuesLength; int mid; while (lo < hi) { mid = lo + ((hi - lo) >> 1); if (values[mid] < lowerTarget) { lo = mid + 1; } else { hi = mid; } } lower = lo; // Find upper bound for upperTarget, starting from the lower bound's result. // Since upperTarget >= lowerTarget is assumed, we can start the search for the upper bound // from the position of the lower bound, effectively reducing the search space. hi = valuesLength; while (lo < hi) { mid = lo + ((hi - lo) >> 1); if (values[mid] <= upperTarget) { lo = mid + 1; } else { hi = mid; } } upper = lo; } /// /// Extremely small, allocation-free prefilter over the last few columns (near p-1). /// /// It uses only the fact that the product must have no carry beyond column p-1 (i.e. carry-out from /// column p-1 is zero) and that bits 0..p-1 of M_p are all ones. /// /// This is intentionally conservative (never rejects a valid divisor). In practice it can still /// remove some q with "awkward" high-bit geometry where the required carry values near the top /// become impossible given the availability of contributions implied by qOneOffsets and maxAllowedA. /// private static bool TryHighBitTopCarryPrefilter( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, long pLong, long maxAllowedA) { // We process last K columns: (p-1), (p-2), ... int columns = HighBitCarryPrefilterColumns; long maxPossibleColumns = pLong - 1; if (columns > maxPossibleColumns) columns = (int)maxPossibleColumns; // Two fixed interval buffers on stack. // Each interval is [lo, hi], disjoint and sorted. // 32 is intentionally small: in practice this stays tiny (often 1-2 intervals). Span lo0 = stackalloc long[64]; Span hi0 = stackalloc long[64]; Span lo1 = stackalloc long[64]; Span hi1 = stackalloc long[64]; // "out" = carryOut set for current column (start at top with {0}) Span outLo = lo0, outHi = hi0; Span inLo = lo1, inHi = hi1; int outCount = 1; outLo[0] = 0; outHi[0] = 0; int inCount = 0; bool flip = false; [MethodImpl(MethodImplOptions.AggressiveInlining)] static void AddInterval(ref int count, Span loArr, Span hiArr, long lo, long hi) { if (hi < lo) return; if (lo < 0) lo = 0; if (hi < 0) return; if (count == 0) { loArr[0] = lo; hiArr[0] = hi; count = 1; return; } // Merge into last interval if overlapping/adjacent. long lastHi = hiArr[count - 1]; if (lo <= lastHi + 1) { if (hi > lastHi) hiArr[count - 1] = hi; return; } // Append if ((uint)count >= (uint)loArr.Length) { // If this ever happens, we conservatively stop filtering (do not reject). // Returning "true" keeps correctness (prefilter is optional). count = 0; return; } loArr[count] = lo; hiArr[count] = hi; count++; } for (int step = 0; step < columns; step++) { long column = maxPossibleColumns - step; if (column < 0) break; // n = count of q-offsets t such that 0 <= column - t <= maxAllowedA // <=> column - maxAllowedA <= t <= column long lowerTLong = column - maxAllowedA; int lowerT = lowerTLong <= 0 ? 0 : (lowerTLong >= int.MaxValue ? int.MaxValue : (int)lowerTLong); int upperT = column <= 0 ? 0 : (column >= int.MaxValue ? int.MaxValue : (int)column); FindBounds(qOneOffsets, lowerT, upperT, out int left, out int right); int n = right - left; if (n < 0) { // Console.WriteLine("n < 0 - excluding in pre-filter"); return false; } // Build carryIn set into inLo/inHi inCount = 0; // For M_p, bit in columns [0..p-1] is 1: // carryIn + ones = 1 + 2*carryOut // carryIn = (2*carryOut + 1) - ones, with ones in [0..min(n, 2*carryOut+1)]. // // That yields carryIn interval: // r = 2*carryOut + 1 // carryIn ∈ [ r - min(n, r), r ] = [ max(0, r-n), r ]. // // Union over carryOut intervals is merged cheaply. for (int i = 0; i < outCount; i++) { long a = outLo[i]; long b = outHi[i]; if (b < 0) continue; if (a < 0) a = 0; if (b < a) continue; const long TwicePlusOneOverflowLimit = (long.MaxValue - 1) >> 1; // When carryOut is already too big, top-down becomes useless in long. // We want to be conservative to avoid rejecting correct divisors => stop filtering with this filter. if (a > TwicePlusOneOverflowLimit || b > TwicePlusOneOverflowLimit) { return true; } // r ranges from rMin..rMax over carryOut in [a..b] long rMin = (a << 1) + 1; long rMax = (b << 1) + 1; if (n <= 0) { // ones=0 only -> carryIn=r (odd sequence). Conservatively interval [rMin..rMax]. AddInterval(ref inCount, inLo, inHi, rMin, rMax); if (inCount == 0) return true; // overflow of intervals -> stop filtering continue; } if (n == 1) { // ones in {0,1} -> carryIn in {r, r-1}. Conservatively [rMin-1 .. rMax]. AddInterval(ref inCount, inLo, inHi, rMin - 1, rMax); if (inCount == 0) return true; continue; } // n >= 2 : union becomes dense enough that interval bounds are safe conservative approximations. if (rMax <= n) { // For all r<=n, lower bound is 0. AddInterval(ref inCount, inLo, inHi, 0, rMax); if (inCount == 0) return true; } else if (rMin > n) { // Always r>n => [r-n .. r] AddInterval(ref inCount, inLo, inHi, rMin - n, rMax); if (inCount == 0) return true; } else { // Crossing: part uses [0..r], later [r-n..r]. For n>=2 union is [0..rMax]. AddInterval(ref inCount, inLo, inHi, 0, rMax); if (inCount == 0) return true; } } if (inCount <= 0) { // Console.WriteLine("inCount <= 0 - excluding in pre-filter"); return false; } // Prepare for next column: carryOut := carryIn (shift down one column) outCount = inCount; // Toggle which buffers are used for out vs in. // Avoid tuple swap / temp Span to satisfy Span lifetime rules. flip = !flip; if (!flip) { outLo = lo0; outHi = hi0; inLo = lo1; inHi = hi1; } else { outLo = lo1; outHi = hi1; inLo = lo0; inHi = hi0; } } return true; } private static int[] GetPow2DiffLut(int mod) { _pow2DiffModCache ??= new Dictionary(); if (_pow2DiffModCache.TryGetValue(mod, out int[]? lut)) { return lut; } int[] values = new int[256]; values[0] = mod > 1 ? 1 : 0; for (int i = 1; i < values.Length; i++) { values[i] = (values[i - 1] << 1) % mod; } _pow2DiffModCache[mod] = values; return values; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ComputeForcedRemainingAndClean_Prefix( in ulong[] qMask, in ulong[] oneWin, in ulong[] zeroWin, int words, ulong lastMask, out int forced, out int remaining) { if (words == 0) { forced = 0; remaining = 0; return; } int i = 0; int wordsMinusOne = words - 1; var vForced = Vector.Zero; var vRemaining = Vector.Zero; if (Vector.IsHardwareAccelerated) { // Console.WriteLine("SIMD accelerator ComputeForcedRemainingAndClean_Prefix"); int limit = wordsMinusOne - (wordsMinusOne % Vector.Count); for (; i < limit; i += Vector.Count) { var vq = new Vector(qMask, i); var vOne = new Vector(oneWin, i); var vZero = new Vector(zeroWin, i); var vKnown = vOne | vZero; var vForcedBits = vOne & vq; var vUnknownBits = vq & ~vKnown; vForced += SimdPopCount(vForcedBits); vRemaining += SimdPopCount(vUnknownBits); } } // else // { // Console.WriteLine("No acceleration for ComputeForcedRemainingAndClean_Prefix"); // } forced = (int)Vector.Sum(vForced); remaining = (int)Vector.Sum(vRemaining); for (; i < wordsMinusOne; i++) { ulong unknownBits = qMask[i]; ulong one = oneWin[i]; ulong known = one | zeroWin[i]; one &= unknownBits; unknownBits &= ~known; forced += PopCount(one); remaining += PopCount(unknownBits); } if (words > 0) { ulong unknownBits = qMask[words - 1] & lastMask; ulong one = oneWin[words - 1]; ulong known = one | zeroWin[words - 1]; one &= unknownBits; unknownBits &= ~known; forced += PopCount(one); remaining += PopCount(unknownBits); } } private static int ComputeQMod(ReadOnlySpan qOneOffsets, int mod) { int count = qOneOffsets.Length; if (count == 0) { return 0; } int[] diffLut = GetPow2DiffLut(mod); int diffLutLength = diffLut.Length; long sum = 0; int prevOffset = qOneOffsets[0]; long pow = PowMod2((ulong)prevOffset, mod); sum += pow; for (int i = 1; i < count; i++) { int offset = qOneOffsets[i]; int diff = offset - prevOffset; if (diff <= 0) { prevOffset = offset; continue; } if (diff < diffLutLength) { pow = pow * diffLut[diff] % mod; } else { pow = pow * PowMod2((ulong)diff, mod) % mod; } sum += pow; if (sum >= mod) { sum %= mod; } prevOffset = offset; } return (int)(sum % mod); } private static bool ExistsCongruentInRange(BigInteger minA, BigInteger maxA, int residue, int mod) { if (minA > maxA) { return false; } int minMod = (int)((minA % mod + mod) % mod); int delta = residue - minMod; if (delta < 0) { delta += mod; } BigInteger candidate = minA + delta; return candidate <= maxA; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int NormalizeMod(int value, int mod) { int normalized = value % mod; return normalized < 0 ? normalized + mod : normalized; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool IsResidueReachable(ulong mask, int residue) { return ((mask >> residue) & 1UL) != 0UL; } private static bool TryPruneWithKnownBitsModResidue( ReadOnlySpan qOneOffsets, long qHash, ulong p, int mod, long maxAllowedA, HashSet knownOneIndices, HashSet knownZeroIndices) { if (maxAllowedA < 0) { return true; } int maxIndex = maxAllowedA > int.MaxValue ? int.MaxValue : (int)maxAllowedA; int bitsNeeded = 0; int bitsSource = mod - 1; while (bitsSource > 0) { bitsNeeded++; bitsSource >>= 1; } if (bitsNeeded <= 0) { bitsNeeded = 1; } int qMod = GetQModCached(qOneOffsets, qHash, mod); int target = ModPow2Minus1(p, mod); if (qMod == 0) { return target != 0; } int inverse = ModPow(qMod, mod - 2, mod); int targetA = (int)((long)target * inverse % mod); int fixedResidue = 0; foreach (int index in knownOneIndices) { if ((uint)index > (uint)maxIndex) { continue; } fixedResidue += PowMod2((ulong)index, mod); if (fixedResidue >= mod) { fixedResidue %= mod; } } targetA = NormalizeMod(targetA - fixedResidue, mod); List knownIndices = new(knownOneIndices.Count + knownZeroIndices.Count); foreach (int index in knownOneIndices) { if ((uint)index <= (uint)maxIndex) { knownIndices.Add(index); } } foreach (int index in knownZeroIndices) { if ((uint)index <= (uint)maxIndex) { knownIndices.Add(index); } } knownIndices.Sort(); ulong reachableMask = 1UL; ulong moduloMask = (1UL << mod) - 1UL; int previousKnown = -1; int knownCount = knownIndices.Count; for (int i = 0; i <= knownCount; i++) { int nextKnown = i == knownCount ? maxIndex + 1 : knownIndices[i]; if (i > 0 && nextKnown == knownIndices[i - 1]) { continue; } int runStart = previousKnown + 1; int runEnd = nextKnown - 1; int runLength = runEnd - runStart + 1; if (runLength >= bitsNeeded) { return false; } for (int bitIndex = runStart; bitIndex <= runEnd; bitIndex++) { int shift = PowMod2((ulong)bitIndex, mod); ulong shiftedMask; if (shift == 0) { shiftedMask = reachableMask; } else { shiftedMask = ((reachableMask << shift) | (reachableMask >> (mod - shift))) & moduloMask; } reachableMask |= shiftedMask; if (IsResidueReachable(reachableMask, targetA)) { return false; } } previousKnown = nextKnown; } return !IsResidueReachable(reachableMask, targetA); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void GetQPrefixMask(int startColumn, int qBitLen, int qWordCount, out int words, out ulong lastMask) { // We want only offsets t <= startColumn. // start column is inclusive for maxT if (startColumn >= qBitLen - 1) { words = qWordCount; lastMask = ulong.MaxValue; // caller can still & lastWordMask separately return; } if (startColumn < 0) { words = 0; lastMask = 0; return; } qBitLen = startColumn >> 6; words = qBitLen + 1; qWordCount = startColumn & 63; lastMask = qWordCount == 63 ? ulong.MaxValue : ((1UL << (qWordCount + 1)) - 1UL); } private static void ComputeARangeFromTopConstraints( long maxAllowedA, bool top0KnownOne, bool top0KnownZero, bool top1KnownOne, long zeroTailStart, out BigInteger minA, out BigInteger maxA) { _dynamicModChecks++; if (zeroTailStart <= maxAllowedA + 1) { zeroTailStart--; } // zeroTailStart is effectiveTop from here. We're reusing variable to limit registry pressure. if (zeroTailStart < 0) { minA = BigInteger.Zero; maxA = BigInteger.Zero; return; } int shift = zeroTailStart + 1L > int.MaxValue ? int.MaxValue : (int)(zeroTailStart + 1L); maxA = Pow2Provider.BigIntegerMinusOne(shift); minA = BigInteger.Zero; if (zeroTailStart == maxAllowedA) { if (top0KnownOne) { shift = maxAllowedA > int.MaxValue ? int.MaxValue : (int)maxAllowedA; } else if (top0KnownZero && top1KnownOne && maxAllowedA - 1L >= 0) { // zeroTailStart is index from here. We're reusing the variable to limit registry pressure. zeroTailStart = maxAllowedA - 1L; shift = zeroTailStart > int.MaxValue ? int.MaxValue : (int)zeroTailStart; } minA = Pow2Provider.BigInteger(shift); } } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ComputeGuaranteedTopPrefixFromInterval( in BigInteger qExact, int qBitLen, int maxAllowedA, out int knownBits, out ulong knownValue) { knownBits = 1; knownValue = 1UL; if (maxAllowedA <= 0) { return; } int windowBits = TopPrefixProbeBits; int maxBits = maxAllowedA + 1; if (windowBits > maxBits) { windowBits = maxBits; } if (windowBits <= 1) { return; } BigInteger scaled = BigInteger.One << (qBitLen + windowBits - 1); BigInteger prefixHighBig = scaled / qExact; if (prefixHighBig <= BigInteger.Zero) { return; } // prefixHighBig fits in windowBits (<= TopPrefixProbeBits), so ulong is sufficient. ulong prefixHigh = (ulong)prefixHighBig; ulong prefixLow = prefixHigh - 1UL; ulong xor = prefixHigh ^ prefixLow; int differingBits = xor == 0UL ? 0 : 64 - BitOperations.LeadingZeroCount(xor); int commonBits = windowBits - differingBits; if (commonBits < 1) { // MSB of a is still guaranteed by construction. return; } knownBits = commonBits; int shift = windowBits - commonBits; knownValue = shift == 0 ? prefixHigh : (prefixHigh >> shift); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TryGetTopPrefixBitValue( long maxAllowedA, int topPrefixKnownBits, ulong topPrefixValue, long bitIndex, out bool bitValue) { bitValue = false; long offsetFromTop = maxAllowedA - bitIndex; if (offsetFromTop < 0 || offsetFromTop >= topPrefixKnownBits) { return false; } int shift = topPrefixKnownBits - 1 - (int)offsetFromTop; bitValue = ((topPrefixValue >> shift) & 1UL) != 0UL; return true; } private static bool TryDynamicPrefixConsistencyPruneLarge( long maxAllowedA, bool top0KnownOne, bool top0KnownZero, bool top1KnownOne, long zeroTailStart, int topPrefixKnownBits, ulong topPrefixValue) { if (topPrefixKnownBits <= 0) { return true; } if (TryGetTopPrefixBitValue(maxAllowedA, topPrefixKnownBits, topPrefixValue, maxAllowedA, out bool top0PrefixValue)) { if ((top0KnownOne && !top0PrefixValue) || (top0KnownZero && top0PrefixValue)) { return false; } } long topIndex1 = maxAllowedA - 1; if (topIndex1 >= 0 && top1KnownOne && TryGetTopPrefixBitValue(maxAllowedA, topPrefixKnownBits, topPrefixValue, topIndex1, out bool top1PrefixValue) && !top1PrefixValue) { return false; } if (zeroTailStart <= maxAllowedA) { long prefixMinIndex = maxAllowedA - topPrefixKnownBits + 1; long start = zeroTailStart > prefixMinIndex ? zeroTailStart : prefixMinIndex; for (long bitIndex = start; bitIndex <= maxAllowedA; bitIndex++) { if (TryGetTopPrefixBitValue(maxAllowedA, topPrefixKnownBits, topPrefixValue, bitIndex, out bool bitValue) && bitValue) { return false; } } } return true; } private static bool IsCleanInQ_Prefix(ulong[] qMask, ulong[] oneWin, ulong[] zeroWin, int words, ulong lastMask) { for (int i = 0; i < words; i++) { ulong q = qMask[i]; if (i == words - 1) q &= lastMask; if ((q & (oneWin[i] | zeroWin[i])) != 0) return false; } return true; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TryPruneWithQTimesAInterval( in BigInteger qExact, in BigInteger minA, in BigInteger maxA, ulong p) { if (minA > BigInteger.Zero) { long minProductBits = (qExact * minA).GetBitLength(); if ((ulong)minProductBits > p) { return true; } } BigInteger maxProduct = qExact * maxA; long maxProductBits = maxProduct.GetBitLength(); if ((ulong)maxProductBits < p) { return true; } if ((ulong)maxProductBits == p) { BigInteger plusOne = maxProduct + BigInteger.One; if (!plusOne.IsPowerOfTwo) { return true; } } return false; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TrySetGlobalKnownBit( int aIndex, bool valueOne, HashSet knownOneIndices, HashSet knownZeroIndices, ref BigInteger globalMinA, ref BigInteger globalMaxA) { if (aIndex < 0) { return true; } if (valueOne) { if (knownZeroIndices.Contains(aIndex)) { return false; } if (knownOneIndices.Add(aIndex)) { globalMinA += Pow2Provider.BigInteger(aIndex); } } else { if (knownOneIndices.Contains(aIndex)) { return false; } if (knownZeroIndices.Add(aIndex)) { globalMaxA -= Pow2Provider.BigInteger(aIndex); } } return true; } private static bool TryDynamicModularRangePrune( ReadOnlySpan qOneOffsets, long qHash, in BigInteger qExact, ulong p, in BigInteger globalKnownMinA, in BigInteger globalKnownMaxA, HashSet globalKnownOneIndices, HashSet globalKnownZeroIndices, long maxAllowedA, bool top0KnownOne, bool top0KnownZero, bool top1KnownOne, long zeroTailStart, int topPrefixKnownBits, ulong topPrefixValue) { _dynamicModChecks++; long effectiveTop = maxAllowedA; if (zeroTailStart <= effectiveTop + 1) { effectiveTop = zeroTailStart - 1; } if (globalKnownMinA > globalKnownMaxA) { _dynamicModFailures++; return false; } if (effectiveTop > 4096) { if (!TryDynamicPrefixConsistencyPruneLarge(maxAllowedA, top0KnownOne, top0KnownZero, top1KnownOne, zeroTailStart, topPrefixKnownBits, topPrefixValue)) { _dynamicModFailures++; #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopPrefixRangeConflict); #endif return false; } if (TryPruneWithQTimesAInterval(qExact, globalKnownMinA, globalKnownMaxA, p)) { _dynamicModFailures++; #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithQTimesAInterval); #endif return false; } return true; } ComputeARangeFromTopConstraints( maxAllowedA, top0KnownOne, top0KnownZero, top1KnownOne, zeroTailStart, out BigInteger minA, out BigInteger maxA); if (topPrefixKnownBits > 0) { int shift = (int)(maxAllowedA + 1L - topPrefixKnownBits); BigInteger prefixBase = new(topPrefixValue); BigInteger prefixMin = prefixBase << shift; BigInteger prefixMax = ((prefixBase + BigInteger.One) << shift) - BigInteger.One; if (prefixMin > minA) { minA = prefixMin; } if (prefixMax < maxA) { maxA = prefixMax; } if (minA > maxA) { _dynamicModFailures++; #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopPrefixRangeConflict); #endif return false; } } if (globalKnownMinA > minA) { minA = globalKnownMinA; } if (globalKnownMaxA < maxA) { maxA = globalKnownMaxA; } if (minA > maxA) { _dynamicModFailures++; return false; } if (TryPruneWithQTimesAInterval(qExact, minA, maxA, p)) { _dynamicModFailures++; #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithQTimesAInterval); #endif return false; } int knownBitsModResidueModuliLength = KnownBitsModResidueModuli.Length; for (int i = 0; i < knownBitsModResidueModuliLength; i++) { int mod = KnownBitsModResidueModuli[i]; if (TryPruneWithKnownBitsModResidue(qOneOffsets, qHash, p, mod, maxAllowedA, globalKnownOneIndices, globalKnownZeroIndices)) { _dynamicModFailures++; #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithKnownBitsModResidue); #endif return false; } } int modLocalPrefilterModuliLength = ModLocalPrefilterModuli.Length; for (int i = 0; i < modLocalPrefilterModuliLength; i++) { int mod = ModLocalPrefilterModuli[i]; int qMod = GetQModCached(qOneOffsets, qHash, mod); int target = ModPow2Minus1(p, mod); if (qMod == 0) { if (target != 0) { return false; } continue; } int inverse = ModPow(qMod, mod - 2, mod); int aResidue = (int)((long)target * inverse % mod); if (!ExistsCongruentInRange(minA, maxA, aResidue, mod)) { _dynamicModFailures++; return false; } } return true; } private static int ModPow2Minus1(ulong p, int mod) { int pow = PowMod2(p, mod); int result = pow - 1; return result < 0 ? result + mod : result; } private static int PowMod2(ulong exp, int mod) { long result = 1; long baseVal = 2; ulong value = exp; while (value > 0) { if ((value & 1) != 0) { result = (result * baseVal) % mod; } baseVal = (baseVal * baseVal) % mod; value >>= 1; } return (int)result; } private static int ModPow(int value, int exp, int mod) { long result = 1 % mod; long baseVal = value % mod; int e = exp; while (e > 0) { if ((e & 1) != 0) { result = (result * baseVal) % mod; } baseVal = (baseVal * baseVal) % mod; e >>= 1; } return (int)result; } private enum CarryDropPruneReason : byte { None = 0, ZeroDrop = 1, Dominance67ThenZeroDrop = 2, Dominance10AndMoreDrop = 3, Dominance15161718Drop = 4, DominanceReached10AndStayedAbove5 = 5, Dominance1011121314 = 6, HighDensityWithShortLowStreaks = 7, } private static readonly BigInteger _testCase = BigInteger.Parse("337429718205079"); // #if WITH_PRUNE_STATS [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void TraceCarryDropHeuristicPrune( int pLong, BigInteger qExact, CarryDropPruneReason pruneReason, bool reached10, bool reachedBelow5, int reached5, int zeroCount) { if (pLong == 138000953 && qExact.Equals(_testCase)) { Console.WriteLine($"CarryDrop target reached: p={pLong}, q={qExact}, reason={pruneReason}, reached10={reached10}, reachedBelow5={reachedBelow5}, reached5={reached5}, zeroCount={zeroCount}"); } } // #endif private static bool TryCarryDropHeuristicPrune( BigInteger qExact, ReadOnlySpan qOneOffsets, int pLong, int maxAllowedA, ulong[] aBitsBuffer, int[] activeAIndicesBuffer, out CarryDropPruneReason pruneReason) { pruneReason = CarryDropPruneReason.None; if (!EnableCarryDropHeuristic || maxAllowedA <= 0 || qOneOffsets[0] != 0) { return false; } int bitLimit = maxAllowedA + 1; if (bitLimit > pLong) { bitLimit = pLong; } if (bitLimit > CarryDropHeuristicMaxBits) { bitLimit = CarryDropHeuristicMaxBits; } if (bitLimit <= 1) { return false; } int wordCount = (bitLimit + 63) >> 6; int qOneOffsetsLength = qOneOffsets.Length; if ((uint)wordCount > (uint)aBitsBuffer.Length || (uint)qOneOffsetsLength > (uint)activeAIndicesBuffer.Length) { return false; } Array.Clear(aBitsBuffer, 0, wordCount); long carry = 0; bool hadPositiveCarry = false; int upperOffsetExclusive = 1; int totalCarry2Count = 0, totalCarry6Count = 0, totalCarry7Count = 0, totalCarry8Count = 0, totalCarry9Count = 0, totalCarry10Count = 0, totalCarry11Count = 0, totalCarry12Count = 0, totalCarry13Count = 0, totalCarry14Count = 0, totalCarry15Count = 0, totalCarry01234Count = 0, totalCarryFirst32Count = 0, totalCarryAbove10 = 0, totalCarryBelow10 = 0, longestStreakBelow10 = 0, shortestStreakBelow10 = int.MaxValue, lastStreakBelow10 = 0, lowestCarryAfter10 = int.MaxValue; // , totalCarry16Count = 0, totalCarry17Count = 0, totalCarry18Count = 0; int totalCarryOtherPositive67Count = 0, totalCarryOtherPositive8910OtherCount = 0, totalCarryOtherPositive1011121314OtherCount = 0, totalCarryOtherPositive15161718OtherCount = 0; bool hadZeroBefore = false; bool reached10 = false; bool reachedBelow5 = false; int reached5 = 0; int zeroCount = 0; bool restoredTo10 = false; ushort first10 = 0; ushort totalTransitionsOver10 = 0; long previousCarryOut = 0; int activeContributionCount = 0; for (ushort blockStart = 0; blockStart < bitLimit; blockStart += CarryDropHeuristicBlockBits) { int blockEnd = blockStart + CarryDropHeuristicBlockBits; if (blockEnd > bitLimit) { blockEnd = bitLimit; } for (ushort col = blockStart; col < blockEnd; col++) { // Optimization: maintain rolling a-indices for active q offsets. // Each active a-index advances by +1 per column, so we avoid recomputing (col - offset) // in the hot loop while preserving the exact same contribution set and carry semantics. while (upperOffsetExclusive < qOneOffsetsLength && qOneOffsets[upperOffsetExclusive] <= col) { activeAIndicesBuffer[activeContributionCount] = col - qOneOffsets[upperOffsetExclusive]; activeContributionCount++; upperOffsetExclusive++; } long sumWithoutCurrent = carry; for (int i = 0; i < activeContributionCount; i++) { int aIndex = activeAIndicesBuffer[i]; ulong word = aBitsBuffer[aIndex >> 6]; sumWithoutCurrent += (long)((word >> (aIndex & 63)) & 1UL); activeAIndicesBuffer[i] = aIndex + 1; } int aBit = (1 - ((int)sumWithoutCurrent & 1)) & 1; if (aBit != 0) { aBitsBuffer[col >> 6] |= 1UL << (col & 63); } long carryOut = (sumWithoutCurrent + aBit) >> 1; if (col < 32) { totalCarryFirst32Count += (int)carryOut; } if ((previousCarryOut < 10 && carryOut >= 10) || (previousCarryOut >= 10 && carryOut < 10)) { totalTransitionsOver10++; } previousCarryOut = carryOut; if (carryOut == 6) { totalCarry6Count++; } else if (carryOut == 7) { totalCarry7Count++; } else if (carryOut > 0) { totalCarryOtherPositive67Count++; } if (carryOut == 8) { totalCarry8Count++; } else if (carryOut == 9) { totalCarry9Count++; } else if (carryOut == 10) { if (!reached10) { reached10 = true; first10 = col; } totalCarry10Count++; } else if (carryOut == 11) { totalCarry11Count++; } else if (carryOut == 12) { totalCarry12Count++; } else if (carryOut == 13) { totalCarry13Count++; } else if (carryOut == 14) { totalCarry14Count++; } else if (carryOut > 0) { totalCarryOtherPositive8910OtherCount++; totalCarryOtherPositive1011121314OtherCount++; } if (reached10) { if (carryOut == 5) { reached5++; restoredTo10 = false; } else if (carryOut < 5) { reachedBelow5 = true; } else if (reached5 >= 1 && carryOut >= 10) { restoredTo10 = true; } if (carryOut < 10) { totalCarryBelow10++; lastStreakBelow10++; if (carryOut < lowestCarryAfter10) { lowestCarryAfter10 = (int)carryOut; } } else if (carryOut >= 10) { totalCarryAbove10++; if (lastStreakBelow10 > longestStreakBelow10) { longestStreakBelow10 = lastStreakBelow10; } if (lastStreakBelow10 < shortestStreakBelow10) { shortestStreakBelow10 = lastStreakBelow10; } if (lastStreakBelow10 > 0) { lastStreakBelow10 = 0; } } } if (carryOut <= 4) { totalCarry01234Count++; if (carryOut == 2) { totalCarry2Count++; } } if (carryOut >= 14) // if (carryOut == 15) { totalCarry15Count++; } // else if (carryOut == 16) // { // totalCarry16Count++; // } // else if (carryOut == 17) // { // totalCarry17Count++; // } // else if (carryOut == 18) // { // totalCarry18Count++; // } else if (carryOut > 0) { totalCarryOtherPositive15161718OtherCount++; } if (hadPositiveCarry && carryOut == 0 && col < maxAllowedA) { if (hadZeroBefore) { pruneReason = CarryDropPruneReason.ZeroDrop; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } hadZeroBefore = true; } if (carryOut == 0) { zeroCount++; } if (carryOut > 0) { hadPositiveCarry = true; } carry = carryOut; } } if (reached10 && zeroCount <= 5 && !(reached5 == 0 && zeroCount <= 4) && !reachedBelow5 && reached5 <= 1 && (reached5 == 0 || restoredTo10) ) { pruneReason = CarryDropPruneReason.DominanceReached10AndStayedAbove5; TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #if WITH_PRUNE_STATS #endif return true; } if (totalCarry12Count > totalCarry10Count && totalCarry12Count > totalCarry11Count && totalCarry12Count > totalCarry13Count && totalCarry12Count > totalCarry14Count && totalCarry11Count > 100 && totalCarry13Count > 100 && (first10 < 34 || totalCarry01234Count < 21) && totalCarryFirst32Count > 118 ) { pruneReason = CarryDropPruneReason.Dominance1011121314; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } var totalCarry67Count = totalCarry6Count + totalCarry7Count; if (totalCarry67Count > totalCarryOtherPositive67Count || (totalCarry67Count > totalCarry8Count && zeroCount < first10 && lowestCarryAfter10 >= 3 && (totalCarry7Count > totalCarry8Count || (totalCarry7Count & 1) == 1))) { pruneReason = CarryDropPruneReason.Dominance67ThenZeroDrop; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } if (totalCarry15Count /*+ totalCarry16Count + totalCarry17Count + totalCarry18Count */> totalCarryOtherPositive15161718OtherCount // totalCarry16Count + totalCarry17Count >= totalCarry15Count + totalCarry18Count) // totalCarry16Count >= 100 && totalCarry17Count >= 100 // /* totalCarry15Count >= 100 &&*/ totalCarry16Count >= 150 && totalCarry17Count >= 150 /*&& totalCarry18Count >= 100*/ ) { pruneReason = CarryDropPruneReason.Dominance15161718Drop; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } if (totalCarry8Count + totalCarry9Count + totalCarry10Count > totalCarryOtherPositive8910OtherCount && totalCarry8Count >= 150 && totalCarry9Count >= 150 && totalCarry10Count >= 150 && (totalCarry10Count >= totalCarry9Count || totalCarry8Count >= totalCarry9Count || totalCarry10Count < 190)) { pruneReason = CarryDropPruneReason.Dominance10AndMoreDrop; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } bool totalCarry8CountIsEven = (totalCarry8Count & 1) == 0; bool totalCarry9CountIsEven = (totalCarry9Count & 1) == 0; if (totalCarry8CountIsEven && totalCarry9CountIsEven && ((totalCarry10Count & 1) == 0 || (totalCarry11Count > totalCarry10Count && totalCarry8Count > 14 && totalCarry11Count > totalCarry9Count && totalCarry11Count > totalCarry8Count)) ) { pruneReason = CarryDropPruneReason.Dominance10AndMoreDrop; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } bool totalCarry10CountIsEven = (totalCarry10Count & 1) == 0; bool totalCarry11CountIsEven = (totalCarry11Count & 1) == 0; if ((!totalCarry8CountIsEven && !totalCarry10CountIsEven && totalCarry8Count + totalCarry9Count < totalCarry10Count + totalCarry11Count && totalCarryFirst32Count > 118) || (totalCarry8CountIsEven && !totalCarry10CountIsEven && totalCarry8Count < 188 && totalCarry8Count + totalCarry9Count > totalCarry10Count + totalCarry11Count) || (((!totalCarry8CountIsEven && !totalCarry11CountIsEven) || (totalCarry8CountIsEven && totalCarry8Count < 66) || (totalCarry8CountIsEven && totalCarry10CountIsEven && totalCarry11CountIsEven && totalCarry8Count < 72) || (!totalCarry11CountIsEven && (!totalCarry8CountIsEven || !totalCarry10CountIsEven) && totalCarry11Count > 16)) && totalCarry9Count + totalCarry11Count > totalCarry10Count + totalCarry12Count)) { pruneReason = CarryDropPruneReason.Dominance10AndMoreDrop; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } if (totalCarry7Count + totalCarry8Count > totalCarry9Count + totalCarry10Count && totalCarry11Count > 100) { pruneReason = CarryDropPruneReason.Dominance10AndMoreDrop; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } if (totalCarryAbove10 >= 840 && shortestStreakBelow10 <= 3 && ( (lastStreakBelow10 == 0 && first10 > 34) || (totalCarryBelow10 < 100 && longestStreakBelow10 <= 12 && totalCarry11Count < 150 && lastStreakBelow10 < 12) ) ) { pruneReason = CarryDropPruneReason.HighDensityWithShortLowStreaks; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } if (reached10 && !reachedBelow5 && ( ( reached5 == 0 && ( totalTransitionsOver10 < 35 || ( lastStreakBelow10 <= 89 && ( zeroCount > 4 || lastStreakBelow10 > 0 || first10 < 34 || totalCarry2Count < 7 || totalCarry11Count < 156 || totalCarry10Count > totalCarry11Count /* || longestStreakBelow10 < 34*/ ) ) || ( totalTransitionsOver10 > 100 && lastStreakBelow10 > 100 ) ) ) || ( reached5 >= 8 && totalTransitionsOver10 > 35 && ( lastStreakBelow10 >= 5 || lastStreakBelow10 == 0 || (lastStreakBelow10 < 5 && totalTransitionsOver10 > 100) || (lastStreakBelow10 < 5 && totalCarryBelow10 < 186) || (totalTransitionsOver10 < 100) ) ) || ( reached5 <= 2 && totalTransitionsOver10 > 35 && (lastStreakBelow10 <= 28 || lastStreakBelow10 > 35 || longestStreakBelow10 < 49 || totalCarryBelow10 > 500) ) || ( reached5 >= 3 && reached5 <= 7 && ( totalTransitionsOver10 > 50 || (totalTransitionsOver10 < 50 && first10 < 33) ) && ( first10 < 58 || first10 > 100 || totalTransitionsOver10 > 150 || (totalTransitionsOver10 > 100 && longestStreakBelow10 < 49) || (totalTransitionsOver10 > 100 && lastStreakBelow10 > 0) || (totalTransitionsOver10 > 100 && totalCarryBelow10 > 400 && zeroCount < 5) || (totalCarryBelow10 < 320) ) ) ) ) { pruneReason = CarryDropPruneReason.DominanceReached10AndStayedAbove5; #if WITH_PRUNE_STATS TraceCarryDropHeuristicPrune(pLong, qExact, pruneReason, reached10, reachedBelow5, reached5, zeroCount); #endif return true; } else if (reached10 && !reachedBelow5 && reached5 > 0) { Console.WriteLine($"Something excluded this filter for p={pLong}, q={qExact}, reached5={reached5}, totalTransitionsOver10={totalTransitionsOver10}, lastStreakBelow10={lastStreakBelow10}, longestStreakBelow10={longestStreakBelow10}, zeroCount={zeroCount}, first10={first10}, totalCarryBelow10={totalCarryBelow10}, totalCarry2={totalCarry2Count}"); } return false; } // internal static bool TryCheckDivisibilityFromDivisor(ulong q, ulong p, out bool divides, out ContradictionReason reason) // { // divides = false; // if (q == 0UL) // { // reason = ContradictionReason.ParityUnreachable; // return true; // } // Span oneOffsets = stackalloc int[64]; // int count = 0; // ulong remaining = q; // int baseIndex = 0; // while (remaining != 0UL) // { // int tz = BitOperations.TrailingZeroCount(remaining); // baseIndex += tz; // oneOffsets[count++] = baseIndex; // remaining >>= tz + 1; // baseIndex++; // } // return TryCheckDivisibilityFromOneOffsets(oneOffsets.Slice(0, count), p, out divides, out reason); // } [ThreadStatic] private static bool _initialized; public static void Initialize() { if (_initialized) { return; } // _privateWorkSet = new(); // _stableUnknownDelta8 = new(); _initialized = true; // Start with a small size and grow on demand. _stableUnknownDelta8ByUnknown = new int[64][]; _privateWorkSet = new PrivateWorkSet(); } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int ComputeForcedSumFromKnownPrefix( ReadOnlySpan qOneOffsets, int col, ulong[] knownOne0) { int forced = 0; // qOneOffsets są rosnące; przerywamy gdy t>col (bo wtedy aIndex<0 => wkład 0) for (int i = 0; i < qOneOffsets.Length; i++) { int t = qOneOffsets[i]; if (t > col) break; int aIndex = col - t; // w prefiksie gwarantujemy 0 <= aIndex <= maxFixed ulong word = knownOne0[aIndex >> 6]; forced += (int)((word >> (aIndex & 63)) & 1UL); } return forced; } /// /// Full column-wise propagation for a given set of q one-bit offsets (LSB=0) and exponent p. /// Returns true if the solver could decide locally (divides or contradiction). /// #if DETAILED_LOG internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, BigInteger qExact, BigInteger k, bool isPowerOfTwo, out bool divides, out ContradictionReason reason) #else internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, BigInteger qExact, BigInteger k, bool isPowerOfTwo, out bool divides) #endif { if (!_initialized) { Initialize(); } PrivateWorkSet workSet = _privateWorkSet; divides = false; // Dbg(p, k, $"[TGT] Enter solver. p={p} k={k} isPowerOfTwo={isPowerOfTwo}"); #if DETAILED_LOG _lastTopDownFailure = null; _statStableAdvanceCalls = 0; _statStableAdvanceCols = 0; _statStableAdvanceColsBatched = 0; _statScalarCols = 0; _statDeltaRequests = 0; _statDeltaHits = 0; _statDeltaBuilds = 0; _statStableAdvanceBigShift = 0; _statCleanChecks = 0; _statCleanFalse = 0; _statFirstDirtyColumn = -1; #endif _dynamicModChecks = 0; _dynamicModFailures = 0; // These conditions will never be true on EvenPerfectBitScanner execution path. // if (p == 0 || qOneOffsets.Length == 0) // { // reason = ContradictionReason.ParityUnreachable; // return true; // } // While the scanner supports p >= 31, it's a tiny subset of all checked exponents. No need to support it. // Exact small-case evaluation to short-circuit both success and failure deterministically. // int highestBit = qOneOffsets[^1]; // if (p <= 64 && highestBit < 63) // { // ulong q = 0; // foreach (int bit in qOneOffsets) // { // q |= 1UL << bit; // } // ulong mersenne = (1UL << (int)p) - 1UL; // divides = q != 0 && (mersenne % q == 0); // reason = divides ? ContradictionReason.None : ContradictionReason.ParityUnreachable; // return true; // } int qOneOffsetsLength = qOneOffsets.Length; int qBitLen = qOneOffsets[qOneOffsetsLength - 1] + 1; if ((ulong)qBitLen >= p) { #if DETAILED_LOG reason = ContradictionReason.TruncatedLength; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTruncatedLength1); #endif return true; } // q from EvenPerfectBitScanner will always have the correct 2kp+1 form. // if (qOneOffsets[0] != 0) // { // reason = ContradictionReason.ParityUnreachable; // return true; // } int pLong = (int)p; int maxAllowedA = pLong - qBitLen; int carryDropBitLimit = maxAllowedA + 1; if (carryDropBitLimit > pLong) { carryDropBitLimit = pLong; } if (carryDropBitLimit > CarryDropHeuristicMaxBits) { carryDropBitLimit = CarryDropHeuristicMaxBits; } int carryDropWordCount = (carryDropBitLimit + 63) >> 6; workSet.EnsureCapacity(carryDropWordCount, 0, qOneOffsetsLength); ulong[] knownOne1 = workSet.KnownOne1; int[] carryDropActiveAIndices = workSet.CarryDropActiveAIndices; Array.Clear(knownOne1, 0, carryDropWordCount); #if false Array.Clear(qMaskWords, 0, qWordCount); #endif #if false int segmentBase = 0; #endif if (maxAllowedA < 0) { #if DETAILED_LOG reason = ContradictionReason.TruncatedLength; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTruncatedLength2); #endif return true; } if (TryCarryDropHeuristicPrune(qExact, qOneOffsets, pLong, maxAllowedA, knownOne1, carryDropActiveAIndices, out CarryDropPruneReason carryDropPruneReason)) { divides = false; #if WITH_PRUNE_STATS if (carryDropPruneReason == CarryDropPruneReason.DominanceReached10AndStayedAbove5) { IncrementPruneCounter(ref _prunedWithCarryDominanceReached10AndStayedAbove5); } else if (carryDropPruneReason == CarryDropPruneReason.Dominance1011121314) { IncrementPruneCounter(ref _prunedWithCarryDominance1011121314); } else if (carryDropPruneReason == CarryDropPruneReason.Dominance67ThenZeroDrop) { IncrementPruneCounter(ref _prunedWithCarryDominance67); } else if (carryDropPruneReason == CarryDropPruneReason.Dominance10AndMoreDrop) { IncrementPruneCounter(ref _prunedWithCarryDominance8910); } else if (carryDropPruneReason == CarryDropPruneReason.Dominance15161718Drop) { IncrementPruneCounter(ref _prunedWithCarryDominance15161718); } else if (carryDropPruneReason == CarryDropPruneReason.ZeroDrop) { IncrementPruneCounter(ref _prunedWithCarryDropHeuristic); } else if (carryDropPruneReason == CarryDropPruneReason.HighDensityWithShortLowStreaks) { IncrementPruneCounter(ref _prunedWithHighDensityWithShortLowStreaksCarry); } #endif #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS if (pLong == 138000953 && qExact.Equals(_testCase)) { Console.WriteLine($"Validating heuristic prune: p={pLong}, q={qExact}, reason={carryDropPruneReason}"); } #endif ValidateHeuristicPruneAgainstPowmod2OrThrow(qExact, pLong, $"CarryDropHeuristic:{carryDropPruneReason}"); return true; } #if false // Force low bits of 'a' from a ≡ -q^{-1} (mod 2^ForcedALowBits). // This is a necessary condition for q | (2^p - 1) when p >= ForcedALowBits. BigInteger mask = Pow2Provider.BigIntegerMinusOne(ForcedALowBits); 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; } byte[] qMod2kBytes = workSet.QMod2kBytes; Span qPrefix = stackalloc ulong[16]; int i, lastBit, t, index, forced, word; #if false ulong wordValue; #endif ulong[] bitMask64 = BitMask64; if (usedPow2Transform) { for (i = 0; i < qOneOffsetsLength; i++) { t = qOneOffsets[i]; #if false qMaskWords[t >> 6] |= bitMask64[t & 63]; #endif if ((uint)t < 1024u) { qPrefix[t >> 6] |= bitMask64[t & 63]; } } } else { Array.Clear(qMod2kBytes, 0, qMod2kBytes.Length); int currentByte = -1; byte currentValue = 0; for (i = 0; i < qOneOffsetsLength; i++) { t = qOneOffsets[i]; if (t < ForcedALowBits) { // lastBit stand for lastByte here. We're reusing variable to limit registry pressure. lastBit = t >> 3; if (lastBit != currentByte) { if (currentByte >= 0) { qMod2kBytes[currentByte] = currentValue; } currentByte = lastBit; currentValue = 0; } currentValue |= (byte)(1 << (t & 7)); } #if false qMaskWords[t >> 6] |= bitMask64[t & 63]; #endif if ((uint)t < 1024u) { qPrefix[t >> 6] |= bitMask64[t & 63]; } } if (currentByte >= 0) { qMod2kBytes[currentByte] = currentValue; } qMod2k = new BigInteger(qMod2kBytes, isUnsigned: true, isBigEndian: false); } ulong[] reversePrefix = workSet.LastQReverse; GetReversePrefix1024(workSet, qPrefix, reversePrefix); #if false ulong[] lastPrefix = workSet.LastQPrefix; for (i = 0; i < 16; i++) { lastPrefix[i] = qPrefix[i]; } workSet.HasLastQPrefix1024 = true; #endif // Compute inv = q^{-1} mod 2^ForcedALowBits using Hensel/Newton iteration in Z/2^kZ. BigInteger[] invByKMod = workSet.InvByKMod; bool[] hasInvByKMod = workSet.HasInvByKMod; int kMod = (int)(k % ForcedALowBits); BigInteger inv; if (hasInvByKMod[kMod]) { inv = invByKMod[kMod]; } else if (workSet.HasLastInvMod2k) { inv = workSet.LastInvMod2k; } else { inv = qMod2k; } BigInteger qTimesInvMod2k = (qMod2k * inv) & mask; inv = (inv * (BigIntegerNumbers.Two - qTimesInvMod2k)) & mask; qTimesInvMod2k = (qMod2k * inv) & mask; inv = (inv * (BigIntegerNumbers.Two - qTimesInvMod2k)) & mask; qTimesInvMod2k = (qMod2k * inv) & mask; inv = (inv * (BigIntegerNumbers.Two - qTimesInvMod2k)) & mask; qTimesInvMod2k = (qMod2k * inv) & mask; if (qTimesInvMod2k != BigInteger.One) { inv = qMod2k; for (i = 0; i < 5; i++) { qTimesInvMod2k = (qMod2k * inv) & mask; inv = (inv * (BigIntegerNumbers.Two - qTimesInvMod2k)) & mask; qTimesInvMod2k = (qMod2k * inv) & mask; inv = (inv * (BigIntegerNumbers.Two - qTimesInvMod2k)) & mask; } } BigInteger aLow = (-inv) & mask; if (isPowerOfTwo) { workSet.LastPow2QMod2k = qMod2k; workSet.LastPow2K = k; workSet.HasLastPow2QMod2k = true; } workSet.LastInvMod2k = inv; workSet.HasLastInvMod2k = true; invByKMod[kMod] = inv; hasInvByKMod[kMod] = true; // a ≡ -aLow (mod 2^ForcedALowBits) // Apply constraints for bits i in [0..min(ForcedALowBits-1, maxAllowedA)]. Span aLowBytes = stackalloc byte[ForcedALowBits >> 3]; if (!aLow.TryWriteBytes(aLowBytes, out int aLowBytesWritten, isUnsigned: true, isBigEndian: false)) { aLowBytesWritten = 0; } int maxKnownA; int maxFixed = maxAllowedA < (ForcedALowBits - 1) ? maxAllowedA : (ForcedALowBits - 1); if (maxFixed >= 0) { int numWordsToProcess = (maxFixed >> 6) + 1; Array.Clear(knownOne0, 0, numWordsToProcess); int bytesToCopy = Math.Min(aLowBytesWritten, numWordsToProcess * 8); if (bytesToCopy > 0) { System.Runtime.InteropServices.MemoryMarshal.Cast(aLowBytes.Slice(0, bytesToCopy)) .CopyTo(new Span(knownOne0, 0, (bytesToCopy + 7) / 8)); } i = 0; if (Vector.IsHardwareAccelerated) { int limit = numWordsToProcess - (numWordsToProcess % Vector.Count); for (; i < limit; i += Vector.Count) { var vKnownOne = new Vector(knownOne0, i); (~vKnownOne).CopyTo(knownZero0, i); } } for (; i < numWordsToProcess; i++) { knownZero0[i] = ~knownOne0[i]; } int lastWordIndex = maxFixed >> 6; int lastBitInWord = maxFixed & 63; if (lastBitInWord < 63) { ulong finalMask = (1UL << (lastBitInWord + 1)) - 1UL; knownOne0[lastWordIndex] &= finalMask; knownZero0[lastWordIndex] &= finalMask; } } // We have learned bits up to maxFixed in the current segment. maxKnownA = -1; if (maxFixed > maxKnownA) { maxKnownA = maxFixed; } CarryRange carryLow = CarryRange.Zero; // ------------------------------------------------------------ // Cheap low-column prefilter (prefix-style): // For the first N columns, all a-indices that contribute to the column are within // the forced low-bit region (aLow). When every contributing bit is known (remaining==0), // we can propagate carry deterministically and reject immediately on contradiction. // This is significantly cheaper than entering the full row-aware window DP. // ------------------------------------------------------------ // maxFixed now serves as preMax. maxFixed = 512 - 1; // 256..512 is a good trade-off; keep <= ForcedALowBits if (maxFixed > maxKnownA) { maxFixed = maxKnownA; } if (maxFixed > maxAllowedA) { maxFixed = maxAllowedA; } // Do NOT clamp to maxOffset: for columns above maxOffset the full set of q=1 positions // still contributes, and aLow still fixes the corresponding a-bits. // Only run if we have at least a few columns to check. if (false && maxFixed >= 8) { Span prefixMask = stackalloc ulong[16]; // index is column here. We're reusing the variable to limit registry pressure. for (index = 0; index <= maxFixed; index++) { // Required result bit for M_p in columns [0..p-1] is 1. // Compute forced ones exactly: sum of a_{col - t} for all q_t=1 with t<=col. forced = 0; if (index <= 1023) { int shift = 1023 - index; int wordShift = shift >> 6; int bitShift = shift & 63; for (word = 0; word < 16; word++) { int src = word + wordShift; ulong value = 0; if (src < 16) { value = reversePrefix[src] >> bitShift; if (bitShift != 0 && src + 1 < 16) { value |= reversePrefix[src + 1] << (64 - bitShift); } } prefixMask[word] = value; } forced = 0; for (word = 0; word < 16; word++) { forced += PopCount(knownOne0[word] & prefixMask[word]); } } else { for (i = 0; i < qOneOffsetsLength; i++) { t = qOneOffsets[i]; if (t > index) { break; } int aIndex = index - t; // For this prefilter we only run where aIndex is guaranteed within [0..maxFixed]. word = aIndex >> 6; ulong m64 = bitMask64[aIndex & 63]; if ((knownOne0[word] & m64) != 0) forced++; } } // remaining==0 here, so possible==forced. Propagate carry range deterministically. if (!TryPropagateCarry(ref carryLow, forced, forced, 1)) { // reason = ContradictionReason.ParityUnreachable; divides = false; #if DETAILED_LOG Console.WriteLine($"[bit-contradiction] Pruned with {nameof(TryPropagateCarry)}, reversed prefix filter"); #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryPropagateCarry1); #endif ValidateHeuristicPruneAgainstPowmod2OrThrow(qExact, pLong, "TryPropagateCarry1"); return true; } } } #endif #if false // Disabled intentionally: all filters/calculations except TryPropagateCarry1, CarryDropHeuristic, Powmod2Tracked. // Original block retained below for easy restore. // ----------------------------- // Force MSB of 'a' to be 1 (last row must be active => a has exact bit-length maxAllowedA+1) // ----------------------------- #endif // Active filters only: CarryDropHeuristic (above), Powmod2 zero-before-last-block, Powmod2Tracked (optimized). if (TryGetPowmod2BlockStep(qExact, pLong, out int blockBits, out BigInteger blockResidue)) { divides = ComputePowmod2TrackedOptimizedWithZeroBeforeLastBlock( pLong, qExact, blockBits, blockResidue, out bool prunedWithZeroBeforeLastBlock); #if WITH_PRUNE_STATS if (!divides) { if (prunedWithZeroBeforeLastBlock) { IncrementPruneCounter(ref _prunedWithPowmod2ZeroBeforeLastBlock); } else { IncrementPruneCounter(ref _prunedWithPowmod2TrackedOptimized); } } #endif } else { divides = BigInteger.ModPow(BigIntegerNumbers.Two, pLong, qExact) == BigInteger.One; #if WITH_PRUNE_STATS if (!divides) { IncrementPruneCounter(ref _prunedWithPowmod2Tracked); } #endif } #if DETAILED_LOG reason = ContradictionReason.None; PrintStats(); #endif return true; } #if DETAILED_LOG internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, BigInteger k, bool isPowerOfTwo, out bool divides, out ContradictionReason reason) => TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, BuildQExact(qOneOffsets), k, isPowerOfTwo, out divides, out reason); internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, out bool divides) => TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, BuildQExact(qOneOffsets), BigInteger.Zero, false, out divides, out _); #endif #if !DETAILED_LOG internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, BigInteger k, bool isPowerOfTwo, out bool divides) => TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, BuildQExact(qOneOffsets), k, isPowerOfTwo, out divides); internal static bool TryCheckDivisibilityFromOneOffsets( ReadOnlySpan qOneOffsets, ulong p, out bool divides) => TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, BuildQExact(qOneOffsets), BigInteger.Zero, false, out divides); #endif [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TryGetPowmod2BlockStep(in BigInteger modulus, int exponent, out int blockBits, out BigInteger blockResidue) { // blockBits = 0; // blockResidue = BigInteger.Zero; // if (modulus <= BigInteger.One) // { // return false; // } // long blockBitsLong = modulus.GetBitLength(); // if (blockBitsLong <= 0 || blockBitsLong > int.MaxValue) // { // return false; // } // blockBits = (int)blockBitsLong; // if (exponent > blockBits * Powmod2AdaptiveMinRatio) // { // int adaptiveBits = (exponent + Powmod2AdaptiveTargetBlocks - 1) / Powmod2AdaptiveTargetBlocks; // if (adaptiveBits > blockBits) // { // blockBits = adaptiveBits; // } blockBits = 262_144; // } blockResidue = Pow2Provider.BigInteger(blockBits) % modulus; // blockResidue = BigInteger.ModPow(BigIntegerNumbers.Two, blockBits, modulus); return true; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool ComputePowmod2TrackedOptimizedWithZeroBeforeLastBlock( int exponent, in BigInteger modulus, int blockBits, in BigInteger blockResidue, out bool prunedWithZeroBeforeLastBlock) { prunedWithZeroBeforeLastBlock = false; int fullBlocks = exponent / blockBits; int fullBlocksMinusOne = fullBlocks - 1; int tailBits = exponent - (fullBlocks * blockBits); BigInteger residue = BigInteger.One; for (int i = 0; i < fullBlocks; i++) { residue *= blockResidue; if (residue >= modulus) { residue %= modulus; } if (i < fullBlocksMinusOne && residue.IsZero) { prunedWithZeroBeforeLastBlock = true; return false; } } if (tailBits > 0) { BigInteger tailResidue = Pow2Provider.BigInteger(tailBits) % modulus; // BigInteger tailResidue = BigInteger.ModPow(BigIntegerNumbers.Two, tailBits, modulus); residue *= tailResidue; if (residue >= modulus) { residue %= modulus; } } return residue == BigInteger.One; } [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static bool ShouldPruneWithFinalPow2Shift(BigInteger prefixResidue, BigInteger modulus, int tailShift) { if (tailShift <= 0) { return false; } if (prefixResidue == BigInteger.One) { return false; } return (prefixResidue << tailShift) < modulus; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void AdvancePow2Residue( ref BigInteger residue, int delta, in BigInteger modulus, ref int cachedDelta, ref BigInteger cachedPow2Delta) { if (delta <= 0) { return; } if (cachedDelta != delta) { cachedPow2Delta = BigInteger.ModPow(BigIntegerNumbers.Two, delta, modulus); cachedDelta = delta; } residue = (residue * cachedPow2Delta) % modulus; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static BigInteger BuildQExact(ReadOnlySpan qOneOffsets) { BigInteger q = BigInteger.Zero; for (int i = 0; i < qOneOffsets.Length; i++) { q |= BigInteger.One << qOneOffsets[i]; } return q; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong BuildQLow64(ReadOnlySpan qOneOffsets) { ulong q = 0; for (int i = 0; i < qOneOffsets.Length; i++) { int bit = qOneOffsets[i]; if ((uint)bit >= 64u) break; q |= 1UL << bit; } return q; } /// /// Modular inverse of odd q modulo 2^64. /// Uses Newton/Hensel lifting: x_{k+1} = x_k * (2 - q*x_k) mod 2^64. /// [MethodImpl(MethodImplOptions.AggressiveInlining)] private static ulong ModInversePow2_64(ulong qOdd) { // qOdd must be odd. // Start with x ≡ 1 (mod 2). A common fast seed is x = qOdd. // Each iteration doubles correct bits. 6 iterations are enough for 64 bits: 1->2->4->8->16->32->64. unchecked { qOdd *= 2UL - qOdd * qOdd; qOdd *= 2UL - qOdd * qOdd; qOdd *= 2UL - qOdd * qOdd; qOdd *= 2UL - qOdd * qOdd; qOdd *= 2UL - qOdd * qOdd; qOdd *= 2UL - qOdd * qOdd; } return qOdd; } [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)] private static bool TryAdvanceZeroColumns(ref CarryRange carry, int columnCount) { const int Step = 30; long carryMin = carry.Min; long carryMax = carry.Max; int add; while (columnCount >= Step) { add = (1 << Step) - 1; carryMin = (carryMin + add) >> Step; carryMax >>= Step; if (carryMin > carryMax) { return false; } columnCount -= Step; } if (columnCount > 0) { add = (1 << columnCount) - 1; carryMin = (carryMin + add) >> columnCount; carryMax >>= columnCount; if (carryMin > carryMax) { return false; } } carry = new CarryRange(carryMin, carryMax); return true; } // [ThreadStatic] // private static Dictionary _stableUnknownDelta8; // Indexed by "unknown" (usually qOneOffsetsLength). Elements are lazy initialized. private static volatile int[][]? _stableUnknownDelta8ByUnknown; private static int[] GetStableUnknownDelta(int unknown) { #if DETAILED_LOG _statDeltaRequests++; #endif // Lazily allocate the outer array. Keep it volatile to publish safely. int[][] cache = _stableUnknownDelta8ByUnknown!; // Ensure capacity for this "unknown". int cacheLength = cache.Length, parity; if ((uint)unknown >= (uint)cacheLength) { // Grow to at least unknown+1, doubling for amortized O(1). parity means newLength here - we're reusing variable to limit registry pressure. parity = cacheLength; while (parity <= unknown) parity <<= 1; var newCache = new int[parity][]; Array.Copy(cache, newCache, cacheLength); // Publish the grown cache (racy publish is OK; contents are immutable after creation). _stableUnknownDelta8ByUnknown = newCache; cache = newCache; } // Fast path: already computed. var cached = cache[unknown]; if (cached != null) { #if DETAILED_LOG _statDeltaHits++; #endif return cached; } // Compute. const int PrefixLength = 1 << ColumnsAtOnce; int[] delta = new int[PrefixLength]; for (int low = 0; low < PrefixLength; low++) { long value = low; // 16 iterations correspond to processing 16 columns at once. It was slower with: 32, 128, 1024 and rolled. for (int i = 0; i < ColumnsAtOnce; i += 4) { parity = (int)(value & 1L); int sMax = (((unknown ^ parity) & 1) != 0) ? unknown : (unknown - 1); value = (value + sMax - 1) >> 1; parity = (int)(value & 1L); sMax = (((unknown ^ parity) & 1) != 0) ? unknown : (unknown - 1); value = (value + sMax - 1) >> 1; parity = (int)(value & 1L); sMax = (((unknown ^ parity) & 1) != 0) ? unknown : (unknown - 1); value = (value + sMax - 1) >> 1; parity = (int)(value & 1L); sMax = (((unknown ^ parity) & 1) != 0) ? unknown : (unknown - 1); value = (value + sMax - 1) >> 1; } delta[low] = (int)((value << ColumnsAtOnce) - low); } #if DETAILED_LOG _statDeltaBuilds++; #endif // Publish into cache. If a race computes it twice, that's acceptable (same deterministic result). cache[unknown] = delta; return delta; } // private static int[] GetStableUnknownDelta8(int unknown) // { // var stableUnknownDelta8Cache = _stableUnknownDelta8; // if (stableUnknownDelta8Cache.TryGetValue(unknown, out int[]? cached)) // { // return cached!; // } // const int PrefixLength = 131072; // int[] delta = new int[PrefixLength]; // for (int low = 0; low < PrefixLength; low++) // { // long value = low; // for (int i = 0; i < 8; i++) // { // int sMax = (((unknown ^ (int)value) & 1) != 0) ? unknown : unknown - 1; // value = (value + sMax - 1) >> 1; // } // delta[low] = (int)((value << 8) - low); // } // stableUnknownDelta8Cache[unknown] = delta; // return delta; // } const int ColumnsAtOnce = 20; #if DETAILED_LOG private static bool TryAdvanceStableUnknown(ref CarryRange carry, int unknown, int columnCount, out ContradictionReason reason) #else private static bool TryAdvanceStableUnknown(ref CarryRange carry, int unknown, int columnCount) #endif { #if DETAILED_LOG _statStableAdvanceCalls++; _statStableAdvanceCols += columnCount; if (columnCount >= 64) _statStableAdvanceBigShift++; #endif if (columnCount <= 0) { #if DETAILED_LOG reason = ContradictionReason.None; _statStableAdvanceColsBatched += columnCount; #endif return true; } if (unknown <= 0) { if (!TryAdvanceZeroColumns(ref carry, columnCount)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } #if DETAILED_LOG _statStableAdvanceColsBatched += columnCount; reason = ContradictionReason.None; #endif return true; } long carryMin = columnCount >= 63 ? 0 : (carry.Min >> columnCount); long carryMax = carry.Max; int remaining = columnCount; int[] delta = GetStableUnknownDelta(unknown); const long CarryDecimalMask = (1L << ColumnsAtOnce) - 1L; int carryMask; while (remaining >= ColumnsAtOnce) { carryMask = (int)(carryMax & CarryDecimalMask); carryMax = (carryMax + delta[carryMask]) >> ColumnsAtOnce; remaining -= ColumnsAtOnce; } while (remaining > 0) { // carryMask is sMax here. We're reusing variable to limit registry pressure. carryMask = (((unknown ^ (int)carryMax) & 1) != 0) ? unknown : unknown - 1; carryMax = (carryMax + carryMask - 1) >> 1; remaining--; } if (carryMax < carryMin) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } #if DETAILED_LOG _statStableAdvanceColsBatched += columnCount; reason = ContradictionReason.None; #endif carry = new CarryRange(carryMin, carryMax); return true; } // ----------------------------- // Row-aware bottom-up block processing // ----------------------------- [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void ShiftLeft1InPlace(ulong[] wordsOne, ulong[] wordsZero, int wordCount, ulong lastWordMask) { if (wordCount > 0) { ulong carryOne = 0, carryZero = 0; for (int i = 0; i < wordCount; i++) { ulong w = wordsOne[i]; wordsOne[i] = (w << 1) | carryOne; carryOne = w >> 63; w = wordsZero[i]; wordsZero[i] = (w << 1) | carryZero; carryZero = w >> 63; } wordCount--; wordsOne[wordCount] &= lastWordMask; wordsZero[wordCount] &= lastWordMask; } } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool AnyAndNonZero(ulong[] a, ulong[] b, int wordCount) { for (int i = 0; i < wordCount; i++) { if ((a[i] & b[i]) != 0) return true; } return false; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int PopCountAnd(ulong[] a, ulong[] b, int wordCount) { int c = 0; for (int i = 0; i < wordCount; i++) { c += BitOperations.PopCount(a[i] & b[i]); } return c; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int PopCountUnknownInQ(ulong[] qMask, ulong[] oneWin, ulong[] zeroWin, int wordCount) { int c = 0; for (int i = 0; i < wordCount; i++) { ulong known = oneWin[i] | zeroWin[i]; ulong unkQ = qMask[i] & ~known; c += BitOperations.PopCount(unkQ); } return c; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TryFindSingleUnknownInQ(ulong[] qMask, ulong[] oneWin, ulong[] zeroWin, int wordCount, out int offset) { // Returns true if exactly one unknown bit among positions where q has 1. offset = -1; int found = 0; for (int w = 0; w < wordCount; w++) { ulong known = oneWin[w] | zeroWin[w]; ulong unk = qMask[w] & ~known; if (unk == 0) { continue; } int pc = BitOperations.PopCount(unk); found += pc; if (found > 1) { return false; } int bit = BitOperations.TrailingZeroCount(unk); offset = (w << 6) + bit; } return found == 1; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int GetAKnownStateRowAware( int aIndex, int segmentBase, int windowSize, in ulong[] knownOne0, in ulong[] knownZero0, in ulong[] knownOne1, in ulong[] knownZero1) { // 0 unknown, 1 one, 2 zero if (aIndex < 0) { return 2; } int rel, w; ulong m; if (aIndex >= segmentBase && aIndex < segmentBase + windowSize) { rel = aIndex - segmentBase; w = rel >> 6; m = 1UL << (rel & 63); rel = (knownOne0[w] & m) != 0 ? 1 : (knownZero0[w] & m) != 0 ? 2 : 0; return rel; } w = segmentBase - windowSize; if (aIndex >= w && aIndex < segmentBase) { rel = aIndex - w; w = rel >> 6; m = 1UL << (rel & 63); rel = (knownOne1[w] & m) != 0 ? 1 : (knownZero1[w] & m) != 0 ? 2 : 0; return rel; } return 0; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void SetAKnownRowAware( int aIndex, bool valueOne, int segmentBase, int windowSize, ulong[] knownOne0, ulong[] knownZero0, ulong[] knownOne1, ulong[] knownZero1) { if (aIndex < 0) return; int rel, w; ulong m; if (aIndex >= segmentBase && aIndex < segmentBase + windowSize) { rel = aIndex - segmentBase; w = rel >> 6; m = 1UL << (rel & 63); if (valueOne) knownOne0[w] |= m; else knownZero0[w] |= m; return; } // windowSize is seg1Start here. We're reusing the variable to limit registry pressure. windowSize = segmentBase - windowSize; if (aIndex >= windowSize && aIndex < segmentBase) { rel = aIndex - windowSize; w = rel >> 6; m = 1UL << (rel & 63); if (valueOne) knownOne1[w] |= m; else knownZero1[w] |= m; } } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TryFindSingleUnknownInQ_Prefix( in ulong[] qMask, in ulong[] oneWin, in ulong[] zeroWin, int words, ulong lastMask, out int offset) { offset = -1; int found = 0; for (int w = 0; w < words; w++) { ulong q = qMask[w]; if (w == words - 1) q &= lastMask; ulong unknown = oneWin[w] | zeroWin[w]; unknown = q & ~unknown; if (unknown == 0) continue; found += PopCount(unknown); if (found > 1) return false; int bit = BitOperations.TrailingZeroCount(unknown); offset = (w << 6) + bit; } return found == 1; } #if DETAILED_LOG private static bool TryProcessBottomUpBlockRowAware( ReadOnlySpan qOneOffsets, in ulong[] qMaskWords, int qWordCount, ulong lastWordMask, int maxAllowedA, int startColumn, int columnCount, ref CarryRange carry, ref int maxKnownA, ref int segmentBase, ref ulong[] knownOne0, ref ulong[] knownOne1, ref ulong[] knownZero0, ref ulong[] knownZero1, int windowSize, int wordCount, ulong[] aOneWin, ulong[] aZeroWin, ref int windowColumn, ulong prime, long qHash, in BigInteger qExact, int topIndex0, int topIndex1, ref bool top0KnownOne, ref bool top1KnownOne, ref bool top0KnownZero, ref bool top1KnownZero, ref int zeroTailStart, int topPrefixKnownBits, ulong topPrefixValue, HashSet globalKnownOneIndices, HashSet globalKnownZeroIndices, ref BigInteger globalKnownMinA, ref BigInteger globalKnownMaxA, out ContradictionReason reason) #else private static bool TryProcessBottomUpBlockRowAware( ReadOnlySpan qOneOffsets, in ulong[] qMaskWords, int qWordCount, ulong lastWordMask, int maxAllowedA, int startColumn, int columnCount, ref CarryRange carry, ref int maxKnownA, ref int segmentBase, ref ulong[] knownOne0, ref ulong[] knownOne1, ref ulong[] knownZero0, ref ulong[] knownZero1, int windowSize, int wordCount, ulong[] aOneWin, ulong[] aZeroWin, ref int windowColumn, ulong prime, long qHash, in BigInteger qExact, int topIndex0, int topIndex1, ref bool top0KnownOne, ref bool top1KnownOne, ref bool top0KnownZero, ref bool top1KnownZero, ref int zeroTailStart, int topPrefixKnownBits, ulong topPrefixValue, HashSet globalKnownOneIndices, HashSet globalKnownZeroIndices, ref BigInteger globalKnownMinA, ref BigInteger globalKnownMaxA) #endif { // columnCount is endColumn from here. We're reusing the variable to limit registry pressure. columnCount += startColumn; int stableUnknown = qOneOffsets.Length; int maxOffset = qOneOffsets[stableUnknown - 1]; ulong[] temp; // Ensure window is aligned. ArgumentOutOfRangeException.ThrowIfNotEqual(windowColumn, startColumn); int maxOffsetPlusOne = maxOffset + 1, qWordCountMinusOne = qWordCount - 1, endColumnMinusOne = columnCount - 1; while (startColumn < columnCount) { int aIndex, remaining, state, step; bool oldTop0KnownOne = top0KnownOne; bool oldTop1KnownOne = top1KnownOne; bool oldTop0KnownZero = top0KnownZero; bool oldTop1KnownZero = top1KnownZero; int oldZeroTailStart = zeroTailStart; bool topChanged = false; // Segment rollover based on current column (because we may set aIndex within [startColumn-maxOffset..startColumn]). int bound = (startColumn / windowSize) * windowSize; if (bound != segmentBase) { if (bound == segmentBase + windowSize) { temp = knownOne0; knownOne0 = knownOne1; knownOne1 = temp; temp = knownZero0; knownZero0 = knownZero1; knownZero1 = temp; Array.Clear(knownOne0, 0, wordCount); Array.Clear(knownZero0, 0, wordCount); } else { // Console.WriteLine("Rebuilding window"); Array.Clear(knownOne0, 0, wordCount); Array.Clear(knownOne1, 0, wordCount); Array.Clear(knownZero0, 0, wordCount); Array.Clear(knownZero1, 0, wordCount); // Also clear window because it may reference bits outside cache. Array.Clear(aOneWin, 0, qWordCount); Array.Clear(aZeroWin, 0, qWordCount); // Build window for current startColumn: bit t corresponds to a_{startColumn - t} for (step = 0; step <= maxOffset; step++) { aIndex = startColumn - step; state = aIndex < 0 || aIndex > maxAllowedA ? 2 : GetAKnownStateRowAware(aIndex, bound, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); if (state == 1) { aOneWin[step >> 6] |= 1UL << (step & 63); } else if (state == 2) { aZeroWin[step >> 6] |= 1UL << (step & 63); } } if (qWordCount > 0) { aOneWin[qWordCountMinusOne] &= lastWordMask; aZeroWin[qWordCountMinusOne] &= lastWordMask; } } segmentBase = bound; } int qWordsEff = qWordCount; ulong qLastEff = lastWordMask; if (startColumn < maxOffset) { GetQPrefixMask(startColumn, maxOffsetPlusOne, qWordCount, out qWordsEff, out qLastEff); } // StableUnknown batching is valid only if the window contains no known bits in the q-active positions. if (startColumn >= maxOffset && startColumn <= maxAllowedA) { remaining = endColumnMinusOne; if (remaining > maxAllowedA) { remaining = maxAllowedA; } remaining = remaining - startColumn + 1; if (remaining > 0) { // Additionally, ensure current window has no known bits (should be true if maxKnownA condition holds). // Conservative batching: ignore row-aware window state and assume all q positions are unknown. // This cannot create false negatives; it may only reduce pruning power. #if DETAILED_LOG _statBatchEntry++; _statBatchCols += remaining; #endif while (remaining >= TailCarryBatchColumns) { #if DETAILED_LOG if (!TryAdvanceStableUnknown(ref carry, stableUnknown, TailCarryBatchColumns, out var propagateReason)) #else if (!TryAdvanceStableUnknown(ref carry, stableUnknown, TailCarryBatchColumns)) #endif { #if DETAILED_LOG reason = propagateReason; Console.WriteLine($"[bit-contradiction] Pruned with {nameof(TryAdvanceStableUnknown)} 1, remaining={remaining}"); #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryAdvanceStableUnknown); #endif return false; } startColumn += TailCarryBatchColumns; windowColumn = startColumn; remaining -= TailCarryBatchColumns; } if (remaining > 0) { #if DETAILED_LOG if (!TryAdvanceStableUnknown(ref carry, stableUnknown, remaining, out var propagateReason)) #else if (!TryAdvanceStableUnknown(ref carry, stableUnknown, remaining)) #endif { #if DETAILED_LOG reason = propagateReason; Console.WriteLine($"[bit-contradiction] Pruned with {nameof(TryAdvanceStableUnknown)} 2, remaining={remaining}"); #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryAdvanceStableUnknown); #endif return false; } startColumn += remaining; windowColumn = startColumn; } // Reset window to unknown (no known bits) because we skipped updating it across columns. Array.Clear(aOneWin, 0, qWordCount); Array.Clear(aZeroWin, 0, qWordCount); continue; } } // Mark out-of-range (negative) a indices as zero in this column. // if (startColumn < maxOffset) // { // int firstOut = startColumn + 1; // for (int t = firstOut; t <= maxOffset; t++) // { // int w = t >> 6; // int b = t & 63; // aZeroWin[w] |= 1UL << b; // } // if (qWordCount > 0) aZeroWin[qWordCount - 1] &= lastWordMask; // } ComputeForcedRemainingAndClean_Prefix(qMaskWords, aOneWin, aZeroWin, qWordsEff, qLastEff, out int forced, out remaining); // OPTION-A - doesn't work and skips valid q candidates // // possible = forced + remaining (later you already do remaining += forced) // int possible = forced + remaining; // // --- ADD THIS PRUNE --- // if (!CouldStillReachTopCarryZero_UpperBound(carry, startColumn, endColumnMinusOne, possible)) // { // #if DETAILED_LOG // reason = ContradictionReason.ParityUnreachable; // #endif // // Console.WriteLine("[bit-contradiction] Pruned with CouldStillReachTopCarryZero_UpperBound"); // #if WITH_PRUNE_STATS // Atomic.Add(ref _prunedWithCouldStillReachTopCarryZero_UpperBound, 1UL); // #endif // return false; // } // // remaining = possible; // OPTION-B - doesn't work with steps > 1. // int possible = forced + remaining; // // --- ADD: best-case carry check to see if we can still reach 0 before we pass endColumnMinusOne --- // long best = carry.Min; // if (best > 0) // { // // We only need a quick check; simulate a few steps or until it hits 0. // int stepsLeft = Math.Min(1, endColumnMinusOne - startColumn + 1); // long c = best; // for (int s = 0; s < stepsLeft && c > 0; s++) // { // c = BestCaseCarryOutOneStep(c, forced, possible, 1); // if (c < 0) // { // #if DETAILED_LOG // reason = ContradictionReason.ParityUnreachable; // #endif // Console.WriteLine("[bit-contradiction] Pruned with BestCaseCarryOutOneStep"); // #if WITH_PRUNE_STATS // Atomic.Add(ref _prunedWithBestCaseCarryOutOneStep, 1UL); // #endif // return false; // } // // For subsequent columns, forced/possible will change, so this is still conservative (weak). // // But even this single-step parity feasibility is a useful prune. // // break; // cheap mode: only 1 step; remove break for stronger prune (but slower) // } // // If even best-case can't reduce at all due to parity infeasibility, prune already handled. // // Full multi-step would require using future columns' forced/possible, which we don't have here cheaply. // } // If exactly one unknown contribution among q=1 positions AND carry parity is fixed, // we can force the corresponding a bit to satisfy the required result bit parity. if (remaining == 1 && ((carry.Min ^ carry.Max) & 1L) == 0) { if (TryFindSingleUnknownInQ_Prefix(qMaskWords, aOneWin, aZeroWin, qWordsEff, qLastEff, out step)) { aIndex = startColumn - step; // Required result bit for M_p in columns [0..p-1] is 1. const int requiredBit = 1; // We need (carryIn + forced + x) % 2 == requiredBit, where x ∈ {0,1}. int b = (int)(carry.Min & 1L); // fixed parity because Min/Max same parity // We're reusing qWordsEff as w from here to limit registry pressure qWordsEff = (requiredBit ^ b ^ (forced & 1)) & 1; bool needOne = qWordsEff != 0; // Check contradiction with already-known value. state = GetAKnownStateRowAware(aIndex, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); if (needOne) { if (state == 2) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; Console.WriteLine($"[bit-contradiction] Pruned with {nameof(TryFindSingleUnknownInQ_Prefix)}, needOne={needOne}, aIndex={aIndex}, state=2"); #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryFindSingleUnknownInQ_Prefix); #endif return false; } SetAKnownRowAware(aIndex, true, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); if (!TrySetGlobalKnownBit(aIndex, true, globalKnownOneIndices, globalKnownZeroIndices, ref globalKnownMinA, ref globalKnownMaxA)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } if (aIndex == topIndex0) { if (top0KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } if (!top0KnownOne) { top0KnownOne = true; topChanged = true; } } else if (aIndex == topIndex1) { if (top1KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } if (!top1KnownOne) { top1KnownOne = true; topChanged = true; } } // Reflect in window at offset t qWordsEff = step >> 6; b = step & 63; aOneWin[qWordsEff] |= 1UL << b; forced++; } else { if (state == 1) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; Console.WriteLine($"[bit-contradiction] Pruned with {nameof(TryFindSingleUnknownInQ_Prefix)}, needOne={needOne}, aIndex={aIndex}, state=1"); #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryFindSingleUnknownInQ_Prefix); #endif return false; } SetAKnownRowAware(aIndex, false, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); if (!TrySetGlobalKnownBit(aIndex, false, globalKnownOneIndices, globalKnownZeroIndices, ref globalKnownMinA, ref globalKnownMaxA)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } if (aIndex == topIndex0) { if (top0KnownOne) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } if (!top0KnownZero) { top0KnownZero = true; topChanged = true; } } else if (aIndex == topIndex1) { if (top1KnownOne) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } if (!top1KnownZero) { top1KnownZero = true; topChanged = true; } } { UpdateZeroTailFromTop(ref zeroTailStart, topIndex0, topIndex1, top0KnownZero, top1KnownZero); } qWordsEff = step >> 6; b = step & 63; aZeroWin[qWordsEff] |= 1UL << b; } remaining = 0; if (aIndex > maxKnownA) { maxKnownA = aIndex; } } } topChanged = topChanged || oldTop0KnownOne != top0KnownOne || oldTop1KnownOne != top1KnownOne || oldTop0KnownZero != top0KnownZero || oldTop1KnownZero != top1KnownZero || oldZeroTailStart != zeroTailStart; bool runDynamicModPrune = topChanged || ((startColumn & DynamicModPeriodicMask) == 0); if (runDynamicModPrune) { if (!TryDynamicModularRangePrune(qOneOffsets, qHash, qExact, prime, globalKnownMinA, globalKnownMaxA, globalKnownOneIndices, globalKnownZeroIndices, maxAllowedA, top0KnownOne, top0KnownZero, top1KnownOne, zeroTailStart, topPrefixKnownBits, topPrefixValue)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryDynamicModularRangePrune); #endif return false; } } if ((top0KnownZero && top0KnownOne) || (top1KnownZero && top1KnownOne)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopStateConflict); #endif return false; } if (topIndex1 < 0) { if (top0KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopBothZeroConflict); #endif return false; } } else if (top0KnownZero && top1KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopBothZeroConflict); #endif return false; } if (remaining == 0) { int requiredCarryParity = 1 ^ (forced & 1); long carryMin = AlignUpToParity(carry.Min, requiredCarryParity); long carryMax = AlignDownToParity(carry.Max, requiredCarryParity); if (carryMin > carryMax) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithCarryParityPrecheck); #endif return false; } carry = new CarryRange(carryMin, carryMax); } long minSum = carry.Min + forced; long maxSum = carry.Max + forced + remaining; const int requiredParity = 1; if (minSum == maxSum) { if ((minSum & 1L) != requiredParity) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithCarryParityPrecheck); #endif return false; } carry = CarryRange.Single((minSum - requiredParity) >> 1); } else { if ((minSum & 1L) != requiredParity && remaining == 0) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithCarryParityPrecheck); #endif return false; } if (!TryPropagateCarry(ref carry, forced, forced + remaining, requiredParity)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; Console.WriteLine($"[bit-contradiction] Pruned with {nameof(TryPropagateCarry)}"); #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryPropagateCarry2); #endif return false; } } // --- SAFE TOP-CARRY PRUNE (uses true remaining columns to p-1) --- // We never reach the case when colsLeftToTop is small. We eliminate candidates very quickly, // with colsLeftToTop ~138 mln. // int colsLeftToTop = topColumnExclusive - (startColumn + 1); // +1 bo kolumnę startColumn już właśnie przetworzyłeś // if (carry.Min > 0 && colsLeftToTop > 0) // { // if (colsLeftToTop < 63) // { // Console.WriteLine("[bit-contradiction] Checking TOP-CARRY PRUNE"); // if ((carry.Min >> colsLeftToTop) > 0) // { // #if DETAILED_LOG // reason = ContradictionReason.ParityUnreachable; // Console.WriteLine("[bit-contradiction] Pruned with TOP-CARRY PRUNE"); // #endif // #if WITH_PRUNE_STATS // Atomic.Add(ref _prunedWithTopCarry, 1UL); // #endif // return false; // } // } // else // { // Console.WriteLine($"[bit-contradiction] TOP-CARRY PRUNE too far away from top {colsLeftToTop}"); // } // // colsLeftToTop >= 63 => shift dałby 0, więc nigdy nie tnij (konserwatywnie) // } // Advance to next column: shift window, then insert status of a_{startColumn+1} into bit0. // remaining is nextA from here. We're reusing variable to limit registry pressure. remaining = startColumn + 1; ShiftLeft1InPlace(aOneWin, aZeroWin, qWordCount, lastWordMask); // Insert bit0 for nextA if (remaining > maxAllowedA) { // state = 2; aZeroWin[0] |= 1UL; } else { state = GetAKnownStateRowAware(remaining, segmentBase, windowSize, knownOne0, knownZero0, knownOne1, knownZero1); if (state == 1) { aOneWin[0] |= 1UL; } else if (state == 2) { aZeroWin[0] |= 1UL; } } startColumn++; windowColumn = startColumn; #if DETAILED_LOG _statScalarCols++; #endif } #if DETAILED_LOG reason = ContradictionReason.None; #endif return true; } private static bool TryProcessBottomUpBlock( ulong prime, ReadOnlySpan qOneOffsets, in int[] runStarts, in int[] runLengths, int runCount, int maxAllowedA, int startColumn, int columnCount, ref CarryRange carry, ref int maxKnownA, ref int segmentBase, ref int runStartIdx, ref int runEndIdx, ulong[] knownOne0, ulong[] knownOne1, ulong[] knownZero0, ulong[] knownZero1, int windowSize, int topIndex0, int topIndex1, ref bool top0KnownOne, ref bool top1KnownOne, ref bool top0KnownZero, ref bool top1KnownZero, ref int zeroTailStart, int topPrefixKnownBits, ulong topPrefixValue #if DETAILED_LOG , out ContradictionReason reason #endif ) { int endColumn = startColumn + columnCount; int qOneOffsetsLength = qOneOffsets.Length; int maxOffset = qOneOffsets[qOneOffsetsLength - 1]; int stableUnknown = qOneOffsetsLength; long qHash = ComputeQOffsetsHash(qOneOffsets); BigInteger qExact = BuildQExact(qOneOffsets); HashSet globalKnownOneIndices = new(); HashSet globalKnownZeroIndices = new(); BigInteger globalKnownMinA = BigInteger.Zero; BigInteger globalKnownMaxA = Pow2Provider.BigIntegerMinusOne(maxAllowedA + 1); int bound, runStart, runEnd; ulong[] temp; // Span temp; while (startColumn < endColumn) { int remaining; if (startColumn >= maxOffset && startColumn <= maxAllowedA && maxKnownA < startColumn - maxOffset) { remaining = endColumn - 1; if (remaining > maxAllowedA) { remaining = maxAllowedA; } remaining = remaining - startColumn + 1; if (remaining >= 16) { runStartIdx = 0; runEndIdx = runCount; while (remaining >= TailCarryBatchColumns) { #if DETAILED_LOG if (!TryAdvanceStableUnknown(ref carry, stableUnknown, TailCarryBatchColumns, out var propagateReason)) #else if (!TryAdvanceStableUnknown(ref carry, stableUnknown, TailCarryBatchColumns)) #endif { #if DETAILED_LOG reason = propagateReason; #endif return false; } startColumn += TailCarryBatchColumns; remaining -= TailCarryBatchColumns; } if (remaining > 0) { #if DETAILED_LOG if (!TryAdvanceStableUnknown(ref carry, stableUnknown, remaining, out var propagateReason)) #else if (!TryAdvanceStableUnknown(ref carry, stableUnknown, remaining)) #endif { #if DETAILED_LOG reason = propagateReason; #endif return false; } remaining = 0; } bound = (startColumn / windowSize) * windowSize; if (bound != segmentBase) { if (bound == segmentBase + windowSize) { // rollover forward by one segment: swap and clear only the new "current" segment temp = knownOne0; knownOne0 = knownOne1; knownOne1 = temp; temp = knownZero0; knownZero0 = knownZero1; knownZero1 = temp; Array.Clear(knownOne0); Array.Clear(knownZero0); } else { // jumped more than one segment: we genuinely lost locality; clear both Array.Clear(knownOne0); Array.Clear(knownOne1); Array.Clear(knownZero0); Array.Clear(knownZero1); } segmentBase = bound; } continue; } } int forced = 0; remaining = 0; int chosenIndex = -1; bool oldTop0KnownOne = top0KnownOne; bool oldTop1KnownOne = top1KnownOne; bool oldTop0KnownZero = top0KnownZero; bool oldTop1KnownZero = top1KnownZero; int oldZeroTailStart = zeroTailStart; bool topChanged = false; bound = (startColumn / windowSize) * windowSize; if (bound != segmentBase) { if (bound == segmentBase + windowSize) { temp = knownOne0; knownOne0 = knownOne1; knownOne1 = temp; temp = knownZero0; knownZero0 = knownZero1; knownZero1 = temp; // knownOne0.Clear(); // knownZero0.Clear(); Array.Clear(knownOne0); Array.Clear(knownZero0); } else { // knownOne0.Clear(); // knownOne1.Clear(); // knownZero0.Clear(); // knownZero1.Clear(); Array.Clear(knownOne0); Array.Clear(knownOne1); Array.Clear(knownZero0); Array.Clear(knownZero1); } segmentBase = bound; } bound = startColumn - maxAllowedA; while (runEndIdx < runCount && runStarts[runEndIdx] <= startColumn) { runEndIdx++; } while (runStartIdx < runEndIdx) { runStart = runStarts[runStartIdx]; runEnd = runStart + runLengths[runStartIdx] - 1; if (runEnd < bound) { runStartIdx++; continue; } break; } int known0Len, runIndex; int seg0Start = segmentBase; for (runIndex = runStartIdx; runIndex < runEndIdx; runIndex++) { runStart = runStarts[runIndex]; if (runStart > startColumn) { break; } runEnd = runStart + runLengths[runIndex] - 1; if (runEnd > startColumn) { runEnd = startColumn; } if (runStart < bound) { runStart = bound; } if (runStart > runEnd) { continue; } int aStart = startColumn - runEnd; runEnd = startColumn - runStart; if (aStart > maxAllowedA) { continue; } if (runEnd > maxAllowedA) { runEnd = maxAllowedA; } int rangeLength = runEnd - aStart + 1; if (rangeLength <= 0) { continue; } int seg0End = seg0Start + windowSize - 1; int seg1Start = seg0Start - windowSize; int seg1End = seg0Start - 1; int overlapStart = aStart > seg0Start ? aStart : seg0Start; int overlap0End = runEnd < seg0End ? runEnd : seg0End; // seg1End is now overlap1End to limit registry pressure seg1End = runEnd < seg1End ? runEnd : seg1End; known0Len = 0; int ones; if (overlapStart <= overlap0End) { runStart = overlapStart - seg0Start; // overlap0Start is now len to limit registry pressure overlapStart = overlap0End - overlapStart + 1; ones = CountBitsInRange(knownOne0, runStart, overlapStart); // overlap0End is now zeros to limit registry pressure overlap0End = CountBitsInRange(knownZero0, runStart, overlapStart); forced += ones; remaining += overlapStart - ones - overlap0End; known0Len = overlapStart; if (chosenIndex == -1 && overlapStart > ones + overlap0End && TryFindFirstUnknown(knownOne0, knownZero0, runStart, overlapStart, out int idx)) { chosenIndex = seg0Start + idx; } } // overlap0End is now knownTotal to limit registry pressure overlap0End = 0; overlapStart = aStart > seg1Start ? aStart : seg1Start; if (overlapStart <= seg1End) { runStart = overlapStart - seg1Start; // overlap1Start is now len to limit registry pressure overlapStart = seg1End - overlapStart + 1; ones = CountBitsInRange(knownOne1, runStart, overlapStart); // seg1End is now zeros to limit registry pressure seg1End = CountBitsInRange(knownZero1, runStart, overlapStart); forced += ones; remaining += overlapStart - ones - seg1End; overlap0End = overlapStart; if (chosenIndex == -1 && overlapStart > ones + seg1End && TryFindFirstUnknown(knownOne1, knownZero1, runStart, overlapStart, out int idx)) { chosenIndex = seg1Start + idx; } } overlap0End = known0Len + overlap0End; if (overlap0End < rangeLength) { remaining += rangeLength - overlap0End; } if (zeroTailStart <= runEnd) { // overlap0End is now tailStart to limit registry pressure overlap0End = aStart > zeroTailStart ? aStart : zeroTailStart; if (overlap0End <= runEnd) { rangeLength = runEnd - overlap0End + 1; if (rangeLength > 0) { remaining -= rangeLength; if (remaining < 0) { remaining = 0; } } } } } long minSum = carry.Min; long maxSum = carry.Max; ulong mask; if (remaining == 1 && ((minSum ^ maxSum) & 1) == 0) { // known0Len is now requiredUnknownParity to limit registry pressure known0Len = 1 ^ ((int)minSum & 1) ^ (forced & 1); if (known0Len != 0) { if (chosenIndex >= 0) { if (chosenIndex >= seg0Start) { bound = chosenIndex - seg0Start; // runIndex is now word to limit registry pressure runIndex = bound >> 6; mask = 1UL << (bound & 63); knownOne0[runIndex] |= mask; knownZero0[runIndex] &= ~mask; } else { bound = chosenIndex - (seg0Start - windowSize); // runIndex is now word to limit registry pressure runIndex = bound >> 6; mask = 1UL << (bound & 63); knownOne1[runIndex] |= mask; knownZero1[runIndex] &= ~mask; } forced++; remaining = 0; if (chosenIndex > maxKnownA) { maxKnownA = chosenIndex; } if (chosenIndex == topIndex0) { if (top0KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } top0KnownOne = true; topChanged = true; } else if (chosenIndex == topIndex1) { if (top1KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } top1KnownOne = true; topChanged = true; } } } else { // if (_debugEnabled && _parityZeroLogCount < 8) // { // _parityZeroLogCount++; // Console.WriteLine($"[bit-contradiction] unknown==1 parity forces zero at column={column} carry=[{minSum},{maxSum}] forced={forced} chosenIndex={chosenIndex}"); // } // Treat the single unknown as 0 for parity. remaining = 0; if (chosenIndex == topIndex0) { if (top0KnownOne) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } top0KnownZero = true; } else if (chosenIndex == topIndex1) { if (top1KnownOne) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } top1KnownZero = true; } UpdateZeroTailFromTop(ref zeroTailStart, topIndex0, topIndex1, top0KnownZero, top1KnownZero); } } topChanged = topChanged || oldTop0KnownOne != top0KnownOne || oldTop1KnownOne != top1KnownOne || oldTop0KnownZero != top0KnownZero || oldTop1KnownZero != top1KnownZero || oldZeroTailStart != zeroTailStart; bool runDynamicModPrune = topChanged || ((startColumn & DynamicModPeriodicMask) == 0); if (runDynamicModPrune) { if (!TryDynamicModularRangePrune(qOneOffsets, qHash, qExact, prime, globalKnownMinA, globalKnownMaxA, globalKnownOneIndices, globalKnownZeroIndices, maxAllowedA, top0KnownOne, top0KnownZero, top1KnownOne, zeroTailStart, topPrefixKnownBits, topPrefixValue)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTryDynamicModularRangePrune); #endif return false; } } if ((top0KnownZero && top0KnownOne) || (top1KnownZero && top1KnownOne)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopStateConflict); #endif return false; } if (topIndex1 < 0) { if (top0KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopBothZeroConflict); #endif return false; } } else if (top0KnownZero && top1KnownZero) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithTopBothZeroConflict); #endif return false; } if (remaining == 0) { // known0Len is now requiredCarryParity to limit registry pressure known0Len = 1 ^ (forced & 1); long carryMin = AlignUpToParity(carry.Min, known0Len); long carryMax = AlignDownToParity(carry.Max, known0Len); if (carryMin > carryMax) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithCarryParityPrecheck); #endif return false; } carry = new CarryRange(carryMin, carryMax); minSum = carryMin; maxSum = carryMax; } minSum += forced; maxSum += forced + remaining; const int requiredParity = 1; if (minSum == maxSum) { if ((minSum & 1) != requiredParity) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithCarryParityPrecheck); #endif return false; } minSum = (minSum - requiredParity) >> 1; carry = CarryRange.Single(minSum); } else { if ((minSum & 1) != requiredParity && remaining == 0) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif #if WITH_PRUNE_STATS IncrementPruneCounter(ref _prunedWithCarryParityPrecheck); #endif return false; } if (!TryPropagateCarry( ref carry, forced, forced + remaining, requiredParity)) { #if DETAILED_LOG reason = ContradictionReason.ParityUnreachable; #endif return false; } } startColumn++; #if DETAILED_LOG _statScalarCols++; #endif } #if DETAILED_LOG reason = ContradictionReason.None; #endif return true; } #if DETAILED_LOG private static bool TryProcessTopDownBlock( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int maxAllowedA, int startHighColumn, int columnCount, CarryRange carryOut, ref int windowStart, ref int windowEnd, out CarryRange nextCarryOut, out TopDownPruneFailure? failure) #else private static bool TryProcessTopDownBlock( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int maxAllowedA, int startHighColumn, int columnCount, CarryRange carryOut, ref int windowStart, ref int windowEnd, out CarryRange nextCarryOut) #endif { #if DETAILED_LOG failure = null; #endif long carryOutMin = carryOut.Min, carryOutMax = carryOut.Max; int value; int endColumn = startHighColumn - columnCount + 1; int column = startHighColumn; while (column >= endColumn) { // value is used as high column index value = column; while (windowEnd > 0 && qOneOffsets[windowEnd - 1] > value) { windowEnd--; } // value is used as low column index value = column - maxAllowedA; if (value < 0) { value = 0; } while (windowStart > 0 && qOneOffsets[windowStart - 1] >= value) { windowStart--; } if (windowStart > windowEnd) { windowStart = windowEnd; } // value is unknown bits count here value = windowEnd - windowStart; if (value == 0) { int lastEnd = windowEnd > 0 ? qOneOffsets[windowEnd - 1] : int.MinValue; int lastStart = windowStart > 0 ? qOneOffsets[windowStart - 1] + maxAllowedA : int.MinValue; int chunkEnd = lastEnd > lastStart ? lastEnd : lastStart; if (chunkEnd < endColumn) { chunkEnd = endColumn; } int steps = column - chunkEnd + 1; if (steps > 1) { if (!AdvanceTopDownUnknown0(ref carryOutMin, ref carryOutMax, steps)) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, carryOutMin, carryOutMax, value); #endif nextCarryOut = default; return false; } column = chunkEnd - 1; continue; } } else if (value == 1) { int lastEnd = windowEnd > 0 ? qOneOffsets[windowEnd - 1] : int.MinValue; int lastStart = windowStart > 0 ? qOneOffsets[windowStart - 1] + maxAllowedA : int.MinValue; int chunkEnd = lastEnd > lastStart ? lastEnd : lastStart; if (chunkEnd < endColumn) { chunkEnd = endColumn; } int steps = (int)(column - chunkEnd + 1); if (steps > 1) { if (!AdvanceTopDownUnknown1(ref carryOutMin, ref carryOutMax, steps)) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, carryOutMin, carryOutMax, value); #endif nextCarryOut = default; return false; } column = chunkEnd - 1; continue; } } // carryOutMin now becomes nextCarryMin to limit registry pressure carryOutMin = (carryOutMin << 1) + 1 - value; // carryOutMax now becomes nextCarryMax to limit registry pressure carryOutMax = (carryOutMax << 1) + 1; if (carryOutMin < 0) { carryOutMin = 0; } if (carryOutMax < 0) { carryOutMax = 0; } if (carryOutMin > carryOutMax) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, carryOutMin, carryOutMax, value); #endif nextCarryOut = default; return false; } column--; } nextCarryOut = new CarryRange(carryOutMin, carryOutMax); return true; } private static long AlignUpToParity(long value, int parity) => (value & 1L) == parity ? value : value + 1L; private static long AlignDownToParity(long value, int parity) => (value & 1L) == parity ? value : value - 1L; #if DETAILED_LOG private static bool TryTopDownCarryPrune( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int pLong, int maxAllowedA, out TopDownPruneFailure? failure) #else private static bool TryTopDownCarryPrune( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int pLong, int maxAllowedA) #endif { #if DETAILED_LOG failure = null; #endif int windowStart = qOneOffsetsLength; int windowEnd = qOneOffsetsLength; long carryOutMin = 0, carryOutMax = 0; int value; int column = pLong - 1; while (column >= 0) { // value is used as high column index value = column; while (windowEnd > 0 && qOneOffsets[windowEnd - 1] > value) { windowEnd--; } // value is used as low column index value = column - maxAllowedA; if (value < 0) { value = 0; } while (windowStart > 0 && qOneOffsets[windowStart - 1] >= value) { windowStart--; } if (windowStart > windowEnd) { windowStart = windowEnd; } // value is used as unknown bits count here value = windowEnd - windowStart; if (value == 0) { int lastEnd = windowEnd > 0 ? qOneOffsets[windowEnd - 1] : int.MinValue; int lastStart = windowStart > 0 ? qOneOffsets[windowStart - 1] + maxAllowedA : int.MinValue; int chunkEnd = lastEnd > lastStart ? lastEnd : lastStart; if (chunkEnd < 0) { chunkEnd = 0; } int steps = (int)(column - chunkEnd + 1); if (steps > 1) { if (!AdvanceTopDownUnknown0(ref carryOutMin, ref carryOutMax, steps)) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, carryOutMin, carryOutMax, value); #endif return false; } column = chunkEnd - 1; continue; } } else if (value == 1) { int lastEnd = windowEnd > 0 ? qOneOffsets[windowEnd - 1] : int.MinValue; int lastStart = windowStart > 0 ? qOneOffsets[windowStart - 1] + maxAllowedA : int.MinValue; int chunkEnd = lastEnd > lastStart ? lastEnd : lastStart; if (chunkEnd < 0) { chunkEnd = 0; } int steps = (int)(column - chunkEnd + 1); if (steps > 1) { if (!AdvanceTopDownUnknown1(ref carryOutMin, ref carryOutMax, steps)) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, carryOutMin, carryOutMax, value); #endif return false; } column = chunkEnd - 1; continue; } } // carryOutMin now becomes nextCarryMin to limit registry pressure carryOutMin = (carryOutMin << 1) + 1 - value; // carryOutMax now becomes nextCarryMax to limit registry pressure carryOutMax = (carryOutMax << 1) + 1; if (carryOutMin < 0) { carryOutMin = 0; } if (carryOutMax < 0) { carryOutMax = 0; } if (carryOutMin > carryOutMax) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, carryOutMin, carryOutMax, value); #endif return false; } column--; } if (carryOutMin > 0) { #if DETAILED_LOG failure = new TopDownPruneFailure(0, carryOutMin, carryOutMax, 0); #endif return false; } return true; } #if DETAILED_LOG private static bool TryTopDownBorrowPrefilter( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int pLong, int maxAllowedA, out TopDownPruneFailure? failure) #else private static bool TryTopDownBorrowPrefilter( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int pLong, int maxAllowedA) #endif { #if DETAILED_LOG failure = null; #endif CarryRange borrow = CarryRange.Zero; int windowStart = qOneOffsetsLength; int windowEnd = qOneOffsetsLength; int column = pLong - 1; int endColumn = column - (TopDownBorrowColumns - 1); if (endColumn < 0) { endColumn = 0; } while (column >= endColumn) { int value = column; while (windowEnd > 0 && qOneOffsets[windowEnd - 1] > value) { windowEnd--; } value = column - maxAllowedA; if (value < 0) { value = 0; } while (windowStart > 0 && qOneOffsets[windowStart - 1] >= value) { windowStart--; } if (windowStart > windowEnd) { windowStart = windowEnd; } int unknown = windowEnd - windowStart; if (!TryPropagateCarry(ref borrow, 0, unknown, 1)) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, borrow.Min, borrow.Max, unknown); #endif return false; } column--; } return true; } private static bool ShouldRunTopDownBorrowPrefilter( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, long pLong, long maxAllowedA) { if (qOneOffsetsLength < TopDownBorrowMinUnknown) { return false; } long column = pLong - 1; long start = column - maxAllowedA; if (start < 0) { start = 0; } FindBounds(qOneOffsets, (int)start, (int)column, out int lower, out int upper); return upper - lower >= TopDownBorrowMinUnknown; } private static bool AdvanceTopDownUnknown1(ref long carryOutMin, ref long carryOutMax, int steps) { while (steps > 0) { int step = steps > 30 ? 30 : steps; int pow = 1 << step; if (carryOutMin > (int.MaxValue >> step) || carryOutMax > ((int.MaxValue - (pow - 1)) >> step)) { carryOutMin = 0; carryOutMax = int.MaxValue; return true; } carryOutMin <<= step; carryOutMax = (carryOutMax << step) + (pow - 1); steps -= step; } return true; } private static bool AdvanceTopDownUnknown0(ref long carryOutMin, ref long carryOutMax, int steps) { while (steps > 0) { int step = steps > 30 ? 30 : steps; int add = (1 << step) - 1; if (carryOutMin > ((int.MaxValue - add) >> step) || carryOutMax > ((int.MaxValue - add) >> step)) { carryOutMin = 0; carryOutMax = int.MaxValue; return true; } carryOutMin = (carryOutMin << step) + add; carryOutMax = (carryOutMax << step) + add; steps -= step; } return true; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int CountBitsInRange(in ulong[] bits, int start, int length) { if (length <= 0) { return 0; } length += start - 1; int firstWord = start >> 6; int lastWord = length >> 6; start &= 63; length &= 63; ulong mask; if (firstWord == lastWord) { mask = (length == 63 ? ulong.MaxValue : ((1UL << (length + 1)) - 1UL)) & (ulong.MaxValue << start); return PopCount(bits[firstWord] & mask); } mask = ulong.MaxValue << start; // start is count from here. We're reusing the variable to limit registry pressure. start = PopCount(bits[firstWord] & mask); firstWord++; for (; firstWord < lastWord; firstWord++) { start += PopCount(bits[firstWord]); } mask = length == 63 ? ulong.MaxValue : ((1UL << (length + 1)) - 1UL); start += PopCount(bits[lastWord] & mask); return start; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static bool TryFindFirstUnknown(in ulong[] knownOne, in ulong[] knownZero, int start, int length, out int index) { if (length <= 0) { index = -1; return false; } int end = start + length - 1; int firstWord = start >> 6; int lastWord = end >> 6; int startBit = start & 63; int endBit = end & 63; ulong firstMask = ulong.MaxValue << startBit; ulong lastMask = endBit == 63 ? ulong.MaxValue : ((1UL << (endBit + 1)) - 1UL); for (int word = firstWord; word <= lastWord; word++) { ulong mask = word == firstWord ? firstMask : ulong.MaxValue; if (word == lastWord) { mask &= lastMask; } ulong known = (knownOne[word] | knownZero[word]) & mask; ulong unknown = mask & ~known; if (unknown != 0) { int bit = BitOperations.TrailingZeroCount(unknown); index = (word << 6) + bit; return true; } } index = -1; return false; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static void UpdateZeroTailFromTop(ref int zeroTailStart, int topIndex0, int topIndex1, bool top0KnownZero, bool top1KnownZero) { if (!top0KnownZero) { return; } if (zeroTailStart > topIndex0) { zeroTailStart = topIndex0; } if (top1KnownZero && topIndex1 == topIndex0 - 1 && zeroTailStart > topIndex1) { zeroTailStart = topIndex1; } } #if DETAILED_LOG private static bool TryRunHighBitAndBorrowPrefiltersCombined( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int pLong, int maxAllowedA, out TopDownPruneFailure? failure) #else private static bool TryRunHighBitAndBorrowPrefiltersCombined( ReadOnlySpan qOneOffsets, int qOneOffsetsLength, int pLong, int maxAllowedA) #endif { #if DETAILED_LOG failure = null; #endif int column = pLong - 1, columnPlusOne = column + 1; int highBitColumns = HighBitCarryPrefilterColumns; if (highBitColumns > columnPlusOne) { highBitColumns = columnPlusOne; } int borrowColumns = 0; bool runBorrowPrefilter = qOneOffsetsLength >= TopDownBorrowMinUnknown; if (runBorrowPrefilter) { borrowColumns = TopDownBorrowColumns; if (borrowColumns > columnPlusOne) { borrowColumns = columnPlusOne; } } int sequentialColumns = highBitColumns > borrowColumns ? highBitColumns : borrowColumns; // Given the above, sequentialColumns will never <= 0, because highBitColumns is always > 0. // if (sequentialColumns <= 0) // { // return true; // } Span lo0 = stackalloc long[64]; Span hi0 = stackalloc long[64]; Span lo1 = stackalloc long[64]; Span hi1 = stackalloc long[64]; Span outLo = lo0, outHi = hi0; Span inLo = lo1, inHi = hi1; int outCount = 1; outLo[0] = 0; outHi[0] = 0; int inCount = 0; bool flip = false; CarryRange borrow = CarryRange.Zero; [MethodImpl(MethodImplOptions.AggressiveInlining)] static void AddInterval(ref int count, Span loArr, Span hiArr, long lo, long hi) { if (hi < lo) return; if (lo < 0) lo = 0; if (hi < 0) return; if (count == 0) { loArr[0] = lo; hiArr[0] = hi; count = 1; return; } long lastHi = hiArr[count - 1]; if (lo <= lastHi + 1) { if (hi > lastHi) hiArr[count - 1] = hi; return; } if ((uint)count >= (uint)loArr.Length) { count = 0; return; } loArr[count] = lo; hiArr[count] = hi; count++; } int step = 0; int upperT; int lowerT = column - maxAllowedA; if (lowerT < 0) { lowerT = 0; } FindBounds(qOneOffsets, lowerT, column, out lowerT, out upperT); int n = upperT - lowerT; if (n < 0) { return false; } if (runBorrowPrefilter && n < TopDownBorrowMinUnknown) { runBorrowPrefilter = false; if (highBitColumns < sequentialColumns) { sequentialColumns = highBitColumns; } } while (true) { if (step < highBitColumns) { // Build carryIn set into inLo/inHi inCount = 0; for (int i = 0; i < outCount; i++) { long rMin = outLo[i]; long rMax = outHi[i]; if (rMax < 0) continue; if (rMin < 0) rMin = 0; if (rMax < rMin) continue; const long TwicePlusOneOverflowLimit = (long.MaxValue - 1) >> 1; if (rMin > TwicePlusOneOverflowLimit || rMax > TwicePlusOneOverflowLimit) { return true; } rMin = (rMin << 1) + 1; rMax = (rMax << 1) + 1; if (n <= 0) { AddInterval(ref inCount, inLo, inHi, rMin, rMax); if (inCount == 0) return true; continue; } if (n == 1) { AddInterval(ref inCount, inLo, inHi, rMin - 1, rMax); if (inCount == 0) return true; continue; } if (rMax <= n) { AddInterval(ref inCount, inLo, inHi, 0, rMax); if (inCount == 0) return true; } else if (rMin > n) { AddInterval(ref inCount, inLo, inHi, rMin - n, rMax); if (inCount == 0) return true; } else { AddInterval(ref inCount, inLo, inHi, 0, rMax); if (inCount == 0) return true; } } if (inCount <= 0) { return false; } // Swap buffers for next step. outCount = inCount; flip = !flip; if (!flip) { outLo = lo0; outHi = hi0; inLo = lo1; inHi = hi1; } else { outLo = lo1; outHi = hi1; inLo = lo0; inHi = hi0; } } if (runBorrowPrefilter && step < borrowColumns) { if (!TryPropagateCarry(ref borrow, 0, n, 1)) { #if DETAILED_LOG failure = new TopDownPruneFailure(column, borrow.Min, borrow.Max, n); #endif return false; } } step++; column--; if (step >= sequentialColumns || column < 0) { break; } lowerT = column - maxAllowedA; if (lowerT < 0) { lowerT = 0; } FindBounds(qOneOffsets, lowerT, column, out lowerT, out upperT); n = upperT - lowerT; if (n < 0) { return false; } } return true; } }