"""Reproduces the kernel computation behind Proposition 4.1 of the note
"The partition exponent in the graded pieces of the Claude 3-SAT algorithm is
optimal for pieces of that shape".

For one-part monomials a = t_n, b = t_m on the Fock space with modes S, the operators
of eq. (47) are t_n[r] = c_n(r) nabla_{-r} for |r| in S, and exchanging the factors in
Q_s(t_n, t_m) changes it by

    sum_{i=0}^{H} (-1)^i C(H, i) g(s - m + i),   g(r) = -r c_n(r) c_m(-r) [|r| in S],

H = n + m. Lemma 9(a) for all s says the H-th finite difference of g vanishes, so g is
a polynomial of degree < H vanishing off +-S. This script computes, for |r| <= R, the
dimension of the space of sequences g supported on +-S with vanishing H-th difference
on every window inside [-R, R] (exact rational elimination).

Expected: S = all positive integers gives H - 1 (polynomials of degree < H with
g(0) = 0); the positive integers without 5 give H - 3 for H >= 4 (roots at 0, +-5);
squares, cubes and odd numbers give 0, so c_n(r) c_m(-r) = 0 for every r.

Usage: python3 kernel_check.py
"""
from fractions import Fraction
from math import comb


def rank(rows):
    M = [[Fraction(x) for x in row] for row in rows]
    r = 0
    ncols = len(M[0]) if M else 0
    for c in range(ncols):
        piv = next((i for i in range(r, len(M)) if M[i][c] != 0), None)
        if piv is None:
            continue
        M[r], M[piv] = M[piv], M[r]
        for i in range(len(M)):
            if i != r and M[i][c] != 0:
                f = M[i][c] / M[r][c]
                M[i] = [x - f * y for x, y in zip(M[i], M[r])]
        r += 1
    return r


def kernel_dim(S, H, R):
    support = [r for r in range(-R, R + 1) if r != 0 and abs(r) in S]
    col = {r: j for j, r in enumerate(support)}
    rows = []
    for s in range(-R, R - H + 1):
        row = [0] * len(support)
        for i in range(H + 1):
            if s + i in col:
                row[col[s + i]] += (-1) ** i * comb(H, i)
        rows.append(row)
    return len(support) - rank(rows)


if __name__ == "__main__":
    R = 40
    families = {
        "all positive integers": set(range(1, R + 1)),
        "all but 5": set(range(1, R + 1)) - {5},
        "squares": {j * j for j in range(1, 7)},
        "cubes": {j ** 3 for j in range(1, 4)},
        "odd numbers": {j for j in range(1, R + 1) if j % 2},
    }
    Hs = list(range(2, 13, 2))
    print(f"kernel dimension, |r| <= {R}")
    print(f"{'S':>22} " + " ".join(f"H={H:<3d}" for H in Hs))
    for name, S in families.items():
        print(f"{name:>22} " + " ".join(f"{kernel_dim(S, H, R):<5d}" for H in Hs))