Files
even-perfect-numbers-scanner/PerfectNumbers.Core/ByDivisor/BitContradictionSolverWithAMultiplier.cs
T

5718 lines
186 KiB
C#
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
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 = [];
#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)
{
if (KnownOne1.Length < wordCount)
{
EnsureUlong(ref KnownOne1, wordCount);
}
#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);
#if false
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static void EnsureInt(ref int[] arr, int needed) => Array.Resize(ref arr, needed);
#endif
}
/// <summary>
/// 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.
/// </summary>
[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;
#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}"
);
}
#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<ulong> prefix, Span<ulong> 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<ulong> prefix, Span<ulong> 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<TopDownCacheKey, TopDownCacheValue>? _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<TopDownCacheKey, bool>? _topDownCache;
#endif
[ThreadStatic]
private static long _dynamicModChecks;
[ThreadStatic]
private static long _dynamicModFailures;
[ThreadStatic]
private static Dictionary<QModCacheKey, int>? _qModCache;
[ThreadStatic]
private static List<TopDownCacheKey>? _topDownCacheOrder;
[ThreadStatic]
private static Dictionary<long, PrefixRunCacheEntry>? _prefixRunCache;
[ThreadStatic]
private static List<long>? _prefixRunCacheOrder;
[ThreadStatic]
private static Dictionary<int, int[]>? _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<ulong> SimdPopCount(Vector<ulong> value)
{
// A classic parallel popcount algorithm (SWAR)
var m1 = new Vector<ulong>(0x5555555555555555); // 0101...
var m2 = new Vector<ulong>(0x3333333333333333); // 0011...
var m4 = new Vector<ulong>(0x0F0F0F0F0F0F0F0F); // 00001111...
var h01 = new Vector<ulong>(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<ulong>.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<ulong>.Zero;
var vRemaining = Vector<ulong>.Zero;
var vContradictions = Vector<ulong>.Zero;
int limit = wordCount - (wordCount % Vector<ulong>.Count);
for (i = 0; i < limit; i += Vector<ulong>.Count)
{
var vq = new Vector<ulong>(qMask, i);
var vOne = new Vector<ulong>(oneWin, i);
var vZero = new Vector<ulong>(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<ulong>.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<int> 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<int> qOneOffsets, long qHash, int mod)
{
Dictionary<QModCacheKey, int>? 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<int> 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<long>(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<int> 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<TopDownCacheKey>(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<bool> multiplicand,
ReadOnlySpan<bool> 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;
}
/// <summary>
/// Mini ping-pong (Option B): compute a conservative allowed carry range after walking
/// top-down for a given number of columns starting at <paramref name="startColumn"/>.
///
/// Semantics:
/// - We treat "carryOut at column = startColumn+1" as the starting set (seed).
/// - Each step processes one column (carryOut -> carryIn) moving downward.
/// - After <paramref name="steps"/> columns, the resulting set corresponds to carry-in
/// at column (startColumn - steps + 1).
/// </summary>
private static bool TryHighBitCarryRangeToColumnJ(
ReadOnlySpan<int> 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<long> lo0 = stackalloc long[64];
Span<long> hi0 = stackalloc long[64];
Span<long> lo1 = stackalloc long[64];
Span<long> hi1 = stackalloc long[64];
Span<long> outLo = lo0, outHi = hi0;
Span<long> 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<long> loArr, Span<long> 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<int> 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<int> 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;
}
/// <summary>
/// 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.
/// </summary>
private static bool TryHighBitTopCarryPrefilter(
ReadOnlySpan<int> 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<long> lo0 = stackalloc long[64];
Span<long> hi0 = stackalloc long[64];
Span<long> lo1 = stackalloc long[64];
Span<long> hi1 = stackalloc long[64];
// "out" = carryOut set for current column (start at top with {0})
Span<long> outLo = lo0, outHi = hi0;
Span<long> 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<long> loArr, Span<long> 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<int, int[]>();
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<ulong>.Zero;
var vRemaining = Vector<ulong>.Zero;
if (Vector.IsHardwareAccelerated)
{
// Console.WriteLine("SIMD accelerator ComputeForcedRemainingAndClean_Prefix");
int limit = wordsMinusOne - (wordsMinusOne % Vector<ulong>.Count);
for (; i < limit; i += Vector<ulong>.Count)
{
var vq = new Vector<ulong>(qMask, i);
var vOne = new Vector<ulong>(oneWin, i);
var vZero = new Vector<ulong>(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<int> 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<int> qOneOffsets,
long qHash,
ulong p,
int mod,
long maxAllowedA,
HashSet<int> knownOneIndices,
HashSet<int> 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<int> 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<int> knownOneIndices,
HashSet<int> 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<int> qOneOffsets,
long qHash,
in BigInteger qExact,
ulong p,
in BigInteger globalKnownMinA,
in BigInteger globalKnownMaxA,
HashSet<int> globalKnownOneIndices,
HashSet<int> 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,
}
private static bool TryCarryDropHeuristicPrune(
ReadOnlySpan<int> qOneOffsets,
int pLong,
int maxAllowedA,
ulong[] aBitsBuffer,
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;
if ((uint)wordCount > (uint)aBitsBuffer.Length)
{
return false;
}
Array.Clear(aBitsBuffer, 0, wordCount);
long carry = 0;
bool hadPositiveCarry = false;
int upperOffsetExclusive = 1;
int qOneOffsetsLength = qOneOffsets.Length;
int totalCarry67Count = 0,
totalCarry8Count = 0,
totalCarry9Count = 0,
totalCarry10Count = 0,
totalCarry11Count = 0,
totalCarry12Count = 0,
totalCarry13Count = 0,
totalCarry14Count = 0,
totalCarry15Count = 0;//, 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;
for (int blockStart = 0; blockStart < bitLimit; blockStart += CarryDropHeuristicBlockBits)
{
int blockEnd = blockStart + CarryDropHeuristicBlockBits;
if (blockEnd > bitLimit)
{
blockEnd = bitLimit;
}
int blockCarry67Count = 0;
int blockCarryOtherPositive67Count = 0;
for (int col = blockStart; col < blockEnd; col++)
{
while (upperOffsetExclusive < qOneOffsetsLength && qOneOffsets[upperOffsetExclusive] <= col)
{
upperOffsetExclusive++;
}
long sumWithoutCurrent = carry;
for (int i = 1; i < upperOffsetExclusive; i++)
{
int aIndex = col - qOneOffsets[i];
ulong word = aBitsBuffer[aIndex >> 6];
sumWithoutCurrent += (long)((word >> (aIndex & 63)) & 1UL);
}
int aBit = (1 - ((int)sumWithoutCurrent & 1)) & 1;
if (aBit != 0)
{
aBitsBuffer[col >> 6] |= 1UL << (col & 63);
}
long carryOut = (sumWithoutCurrent + aBit) >> 1;
if (carryOut == 6 || carryOut == 7)
{
blockCarry67Count++;
}
else if (carryOut > 0)
{
blockCarryOtherPositive67Count++;
}
// if (carryOut == 8)
// {
// totalCarry8Count++;
// }
// else if (carryOut == 9)
// {
// totalCarry9Count++;
// }
// else
if (carryOut == 10)
{
if (!reached10)
{
reached10 = true;
}
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++;
}
else if (carryOut < 5)
{
reachedBelow5 = true;
}
}
if (carryOut >= 15)
// 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;
return true;
}
hadZeroBefore = true;
}
if (carryOut == 0)
{
zeroCount++;
}
if (carryOut > 0)
{
hadPositiveCarry = true;
}
carry = carryOut;
}
totalCarry67Count += blockCarry67Count;
totalCarryOtherPositive67Count += blockCarryOtherPositive67Count;
}
if (reached10 && !reachedBelow5 && reached5 <= 1 && zeroCount <= 5)
{
pruneReason = CarryDropPruneReason.DominanceReached10AndStayedAbove5;
return true;
}
if (totalCarry12Count > totalCarry10Count && totalCarry12Count > totalCarry11Count &&
totalCarry12Count > totalCarry13Count && totalCarry12Count > totalCarry14Count &&
totalCarry11Count > 100 && totalCarry13Count > 100
)
{
pruneReason = CarryDropPruneReason.Dominance1011121314;
return true;
}
if (totalCarry67Count > totalCarryOtherPositive67Count)
{
pruneReason = CarryDropPruneReason.Dominance67ThenZeroDrop;
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;
return true;
}
// if (totalCarry8Count + totalCarry9Count + totalCarry10Count > totalCarryOtherPositive8910OtherCount &&
// totalCarry8Count >= 150 && totalCarry9Count >= 150 && totalCarry10Count >= 150)
// {
// pruneReason = CarryDropPruneReason.Dominance10AndMoreDrop;
// return true;
// }
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<int> 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<int> 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;
}
/// <summary>
/// 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).
/// </summary>
#if DETAILED_LOG
internal static bool TryCheckDivisibilityFromOneOffsets(
ReadOnlySpan<int> qOneOffsets,
ulong p,
BigInteger qExact,
BigInteger k,
bool isPowerOfTwo,
out bool divides,
out ContradictionReason reason)
#else
internal static bool TryCheckDivisibilityFromOneOffsets(
ReadOnlySpan<int> 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);
ulong[] knownOne1 = workSet.KnownOne1;
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(qOneOffsets, pLong, maxAllowedA, knownOne1, 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);
}
#endif
#if DETAILED_LOG
reason = ContradictionReason.ParityUnreachable;
#endif
// #if WITH_PRUNE_STATS
ValidateHeuristicPruneAgainstPowmod2OrThrow(qExact, pLong, $"CarryDropHeuristic:{carryDropPruneReason}");
// #endif
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<ulong> 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<byte> 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<byte, ulong>(aLowBytes.Slice(0, bytesToCopy))
.CopyTo(new Span<ulong>(knownOne0, 0, (bytesToCopy + 7) / 8));
}
i = 0;
if (Vector.IsHardwareAccelerated)
{
int limit = numWordsToProcess - (numWordsToProcess % Vector<ulong>.Count);
for (; i < limit; i += Vector<ulong>.Count)
{
var vKnownOne = new Vector<ulong>(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<ulong> 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);
ValidateHeuristicPruneAgainstPowmod2OrThrow(qExact, pLong, "TryPropagateCarry1");
#endif
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<int> 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<int> 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<int> qOneOffsets,
ulong p,
BigInteger k,
bool isPowerOfTwo,
out bool divides)
=> TryCheckDivisibilityFromOneOffsets(qOneOffsets, p, BuildQExact(qOneOffsets), k, isPowerOfTwo, out divides);
internal static bool TryCheckDivisibilityFromOneOffsets(
ReadOnlySpan<int> 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<int> 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<int> 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;
}
/// <summary>
/// Modular inverse of odd q modulo 2^64.
/// Uses Newton/Hensel lifting: x_{k+1} = x_k * (2 - q*x_k) mod 2^64.
/// </summary>
[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<int, int[]> _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<int> 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<int> globalKnownOneIndices,
HashSet<int> globalKnownZeroIndices,
ref BigInteger globalKnownMinA,
ref BigInteger globalKnownMaxA,
out ContradictionReason reason)
#else
private static bool TryProcessBottomUpBlockRowAware(
ReadOnlySpan<int> 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<int> globalKnownOneIndices,
HashSet<int> 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<int> 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<int> globalKnownOneIndices = new();
HashSet<int> globalKnownZeroIndices = new();
BigInteger globalKnownMinA = BigInteger.Zero;
BigInteger globalKnownMaxA = Pow2Provider.BigIntegerMinusOne(maxAllowedA + 1);
int bound, runStart, runEnd;
ulong[] temp;
// Span<ulong> 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<int> 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<int> 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<int> qOneOffsets,
int qOneOffsetsLength,
int pLong,
int maxAllowedA,
out TopDownPruneFailure? failure)
#else
private static bool TryTopDownCarryPrune(
ReadOnlySpan<int> 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<int> qOneOffsets,
int qOneOffsetsLength,
int pLong,
int maxAllowedA,
out TopDownPruneFailure? failure)
#else
private static bool TryTopDownBorrowPrefilter(
ReadOnlySpan<int> 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<int> 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<int> qOneOffsets,
int qOneOffsetsLength,
int pLong,
int maxAllowedA,
out TopDownPruneFailure? failure)
#else
private static bool TryRunHighBitAndBorrowPrefiltersCombined(
ReadOnlySpan<int> 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<long> lo0 = stackalloc long[64];
Span<long> hi0 = stackalloc long[64];
Span<long> lo1 = stackalloc long[64];
Span<long> hi1 = stackalloc long[64];
Span<long> outLo = lo0, outHi = hi0;
Span<long> 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<long> loArr, Span<long> 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;
}
}