#!/usr/bin/env python3
"""Exact replay of the finite identities (Proposition "Finite identities", prop:finite) from native-data.json.

usage:  python3 finite-identities.py native-data.json [--base-only]

Everything is exact: numbers lie in Q(sqrt-1, sqrt3), stored as integer 5-tuples
(a, b, c, d, N) = (a + b i + c s + d i s)/N with i^2 = -1, s^2 = 3.  The parameters enter as
Laurent polynomials in U_j = xi_j / xi_*j (j = 0..3), on the torus |U_j| = 1, where the
complex conjugate of U^e is U^-e.  Two Laurent polynomials agree on the torus (S^1)^4 iff
they agree coefficient by coefficient, so each identity below is checked identically in xi
(and, with --base-only, only at xi = xi_*, i.e. U = 1).

Built only from the paper: the octonion rule (eq:octonionrule), the coordinate map and
eq:literal for P(xi), the residual rows of Appendix C, and the written-out H, T of Appendix C
(typed below).  From the data file: the 24 covectors, the transport recipe of h, the 24
matrices V^(j) and the redundant seed columns.  Standard library only.

Also: the file's formula for V_xi u is compared, exactly and identically in
xi, with its definition as the sum of the four family velocities (eq:anchored-family with
g = h = exp(t n_i(u_i)), n_i (fraktur n) from the table of Appendix C, transported by eq:anchorbinding with
G_i, L_i, R_i, zeta_i, s_i, phi_i, psi_i, iota_i of Appendix A).
"""
import json, sys, time, hashlib
from math import gcd

# ---------------------------------------------------------------- the field Q(i, sqrt3)
def kn(a, b, c, d, N):
    if N < 0: a, b, c, d, N = -a, -b, -c, -d, -N
    g = gcd(gcd(gcd(a, b), gcd(c, d)), N)
    if g > 1: a, b, c, d, N = a // g, b // g, c // g, d // g, N // g
    return (a, b, c, d, N)
K0 = (0, 0, 0, 0, 1); K1 = (1, 0, 0, 0, 1); KI = (0, 1, 0, 0, 1)
def kadd(x, y):
    a, b, c, d, N = x; e, f, g, h, M = y
    if N == M: return kn(a + e, b + f, c + g, d + h, N)
    return kn(a * M + e * N, b * M + f * N, c * M + g * N, d * M + h * N, N * M)
def kmul(x, y):
    a, b, c, d, N = x; e, f, g, h, M = y
    return kn(a*e - b*f + 3*c*g - 3*d*h, a*f + b*e + 3*c*h + 3*d*g,
              a*g + c*e - b*h - d*f, a*h + d*e + b*g + c*f, N * M)
def kneg(x): return (-x[0], -x[1], -x[2], -x[3], x[4])
def kconj(x): return (x[0], -x[1], x[2], -x[3], x[4])
def kiszero(x): return x[0] == 0 and x[1] == 0 and x[2] == 0 and x[3] == 0
def kinv_unit_gauss(x):   # inverse of a unit Gaussian rational (a + b i)/N, |.| = 1: its conjugate
    assert x[2] == 0 and x[3] == 0 and x[0] * x[0] + x[1] * x[1] == x[4] * x[4]
    return kconj(x)
def kstr(x):
    a, b, c, d, N = x; return f'({a}+{b}i+{c}s+{d}is)/{N}'

# ---------------------------------------------------------------- Laurent polynomials over K
Z4 = (0, 0, 0, 0)
def L(c, e=Z4): return {} if kiszero(c) else {e: c}
def ladd(p, q):
    r = dict(p)
    for e, c in q.items():
        if e in r:
            s = kadd(r[e], c)
            if kiszero(s): del r[e]
            else: r[e] = s
        else: r[e] = c
    return r
def lacc(r, p, q):        # r += p*q in place
    for e1, c1 in p.items():
        for e2, c2 in q.items():
            e = (e1[0]+e2[0], e1[1]+e2[1], e1[2]+e2[2], e1[3]+e2[3]); c = kmul(c1, c2)
            if e in r:
                s = kadd(r[e], c)
                if kiszero(s): del r[e]
                else: r[e] = s
            else: r[e] = c
def lmul(p, q):
    r = {}; lacc(r, p, q); return r
def lscale(p, k): return {} if kiszero(k) else {e: kmul(c, k) for e, c in p.items()}
def lconj(p): return {(-e[0], -e[1], -e[2], -e[3]): kconj(c) for e, c in p.items()}  # on the torus
def lre(p): return lscale(ladd(p, lconj(p)), (1, 0, 0, 0, 2))
def lim(p): return lscale(ladd(p, lscale(lconj(p), (-1, 0, 0, 0, 1))), (0, -1, 0, 0, 2))  # (p - pbar)/(2i)
BASE_ONLY = False
def lspec(p):             # restrict to U = 1 when --base-only
    if not BASE_ONLY: return p
    s = K0
    for c in p.values(): s = kadd(s, c)
    return L(s)

# ---------------------------------------------------------------- octonions (eq:octonionrule)
# x = c + z1 j + z2 l + z3 n, real coordinates (Re c, Im c, Re z1, Im z1, Re z2, Im z2, Re z3, -Im z3)
def to_cz(x):
    return complex(x[0], x[1]), [complex(x[2], x[3]), complex(x[4], x[5]), complex(x[6], -x[7])]
def from_cz(c, z):
    return [c.real, c.imag, z[0].real, z[0].imag, z[1].real, z[1].imag, z[2].real, -z[2].imag]
def omul_rule(x, y):
    c, z = to_cz(x); d, w = to_cz(y)
    dot = lambda p, q: p[0]*q[0] + p[1]*q[1] + p[2]*q[2]
    cr = [z[1]*w[2] - z[2]*w[1], z[2]*w[0] - z[0]*w[2], z[0]*w[1] - z[1]*w[0]]
    return from_cz(c*d - dot(z, [q.conjugate() for q in w]),
                   [c*w[r] + d.conjugate()*z[r] + cr[r].conjugate() for r in range(3)])
OT = {}                    # (a, b) -> list of (k, int): e_a e_b = sum OT e_k
for a in range(8):
    for b in range(8):
        ea = [1 if m == a else 0 for m in range(8)]; eb = [1 if m == b else 0 for m in range(8)]
        v = omul_rule(ea, eb)
        assert all(abs(t - round(t)) < 1e-12 for t in v)
        OT[a, b] = [(k, int(round(t))) for k, t in enumerate(v) if round(t) != 0]
assert OT[1, 2] == [(3, 1)] and OT[1, 4] == [(5, 1)] and OT[2, 4] == [(6, 1)] and OT[1, 6] == [(7, -1)]
OCONJ = [1, -1, -1, -1, -1, -1, -1, -1]

# ---------------------------------------------------------------- Jordan product on 27 coordinates
# coordinates (alpha_1, alpha_2, alpha_3, o01, o02, o12); y has o01 at (0,1), o02 at (0,2), o12 at (1,2)
OFF = {(0, 1): 3, (0, 2): 11, (1, 2): 19}
def entry_terms(r, s):     # matrix entry (r,s) as list of (coordinate index, basis m, sign)
    if r == s: return [(r, 0, 1)]
    if r < s: return [(OFF[r, s] + m, m, 1) for m in range(8)]
    return [(OFF[s, r] + m, m, OCONJ[m]) for m in range(8)]
def build_jordan():
    JT = {}                # (k, a, b) -> Fraction-free: stored as (num, 2) i.e. value num/2
    for r in range(3):
        for s in range(3):
            if r > s: continue
            outs = [(r, 0)] if r == s else [(OFF[r, s] + m, m) for m in range(8)]
            for t in range(3):
                for (ia, ma, sa) in entry_terms(r, t):
                    for (ib, mb, sb) in entry_terms(t, s):
                        for k_, val in OT[ma, mb]:
                            for (ko, mo) in outs:
                                if mo == k_:
                                    # XY + YX: X_{rt} Y_{ts} contributes to both orderings (a,b), (b,a)
                                    for key in ((ko, ia, ib), (ko, ib, ia)):
                                        JT[key] = JT.get(key, 0) + sa * sb * val
    # (XY+YX)/2 = sum_ab JT/2 x_a y_b ; JT now holds 2x the structure constant summed symmetrically
    out = {}
    for (k, a, b), v in JT.items():
        if v: out[k, a, b] = kn(v, 0, 0, 0, 2)
    return out
JT = build_jordan()
JTk = {}
for (k, a, b), v in JT.items(): JTk.setdefault(k, []).append((a, b, v))
JTab = {}
for (k, a, b), v in JT.items(): JTab.setdefault((a, b), []).append((k, v))
def jordan(x, y):          # x, y: dict coord -> Laurent
    out = {}
    for a, pa in x.items():
        for b, pb in y.items():
            for k, v in JTab.get((a, b), ()):
                r = out.setdefault(k, {}); lacc(r, lscale(pa, v), pb)
    return {k: p for k, p in out.items() if p}
TRG = {}
for (k, a, b), v in JT.items():
    if k < 3: TRG[a, b] = kadd(TRG.get((a, b), K0), v)
TRGa = {}
for (a, b), v in TRG.items():
    if not kiszero(v): TRGa.setdefault(a, []).append((b, v))
def trace_form(x, y):      # tr(x o y)
    r = {}
    for a, pa in x.items():
        for b, v in TRGa.get(a, ()):
            if b in y: lacc(r, lscale(pa, v), y[b])
    return r

# ---------------------------------------------------------------- the data file
D = json.load(open(sys.argv[1]))
BASE_ONLY = '--base-only' in sys.argv
def coef(A, B, C, Dd, N): return kn(A, B, C, Dd, N)
XS = [kn(a, b, 0, 0, d) for a, b, d in D['base_phase']['entries']]
assert [(x[0], x[1], x[4]) for x in XS] == [(3, 4, 5), (5, 12, 13), (15, 8, 17), (21, 20, 29)]
def xi_power(j, e):        # xi_j^e as a Laurent polynomial in U: xi_*j^e U_j^e
    c = K1; b = XS[j] if e > 0 else kinv_unit_gauss(XS[j])
    for _ in range(abs(e)): c = kmul(c, b)
    ex = [0, 0, 0, 0]; ex[j] = e
    return L(c, tuple(ex))
def xi_mono(es):
    r = L(K1)
    for j, e in enumerate(es):
        if e: r = lmul(r, xi_power(j, e))
    return lspec(r)
def u_mono(es):
    return lspec(L(K1, tuple(es)))

# ---------------------------------------------------------------- the configuration P(xi) (eq:literal)
OMEGA = kn(-1, 0, 0, 1, 2)              # (-1 + i s)/2
def opow(t):
    c = K1
    for _ in range(t % 3): c = kmul(c, OMEGA)
    return c
def oct_scalar_unit(cpoly, unit):       # octonion (c * unit), unit in {1, j, l, n}, c Laurent (complex)
    re, im = lre(cpoly), lim(cpoly)
    pos = {'1': (0, 1, 1), 'j': (2, 3, 1), 'l': (4, 5, 1), 'n': (6, 7, -1)}[unit]
    x = {}
    if re: x[pos[0]] = re
    if im: x[pos[1]] = lscale(im, (pos[2], 0, 0, 0, 1))
    return x
def omul_poly(x, y):
    out = {}
    for a, pa in x.items():
        for b, pb in y.items():
            for k, v in OT[a, b]:
                r = out.setdefault(k, {}); lacc(r, lscale(pa, (v, 0, 0, 0, 1)), pb)
    return {k: p for k, p in out.items() if p}
def oconj_poly(x): return {m: (p if m == 0 else lscale(p, (-1, 0, 0, 0, 1))) for m, p in x.items()}
def point(u, w):            # p(1,u,w) = v v^*/3: o01 = conj(u)/3, o02 = conj(w)/3, o12 = u conj(w)/3
    third = (1, 0, 0, 0, 3); X = {t: L(third) for t in range(3)}
    for base, o in ((3, oconj_poly(u)), (11, oconj_poly(w)), (19, omul_poly(u, oconj_poly(w)))):
        for m, p in o.items():
            q = lspec(lscale(p, third))
            if q: X[base + m] = q
    return X
def config():
    P = [None] * 39
    for t in range(3): P[t] = {t: L(K1)}
    for s in range(3):
        for t in range(3):
            P[3*(1+s)+t] = point(oct_scalar_unit(L(opow(t)), '1'), oct_scalar_unit(L(opow(s-t)), '1'))
            P[3*(4+s)+t] = point(oct_scalar_unit(L(opow(t)), 'j'), oct_scalar_unit(L(opow(s+t)), 'l'))
            P[3*(7+s)+t] = point(oct_scalar_unit(lscale(xi_mono((1, 0, 0, 0)), opow(t)), 'l'),
                                 oct_scalar_unit(lscale(xi_mono((0, 1, 0, 0)), opow(s+t)), 'n'))
            P[3*(10+s)+t] = point(oct_scalar_unit(lscale(xi_mono((0, 0, 1, 0)), opow(t)), 'n'),
                                  oct_scalar_unit(lscale(xi_mono((0, 0, 0, 1)), opow(s+t)), 'j'))
    return P

def pairrow(p, q): return 1443 + p*(77-p)//2 + q - p - 1
GAMMA = lambda p, q: K0 if p // 3 == q // 3 else (1, 0, 0, 0, 3)
def residual(X):            # F(X) as dict row -> Laurent (nonzero rows only)
    out = {}
    unit = {t: L(K1) for t in range(3)}
    for p in range(39):
        sq = jordan(X[p], X[p])
        for k in range(27):
            v = ladd(sq.get(k, {}), lscale(X[p].get(k, {}), (-1, 0, 0, 0, 1)))
            if v: out[27*p+k] = v
        v = ladd(ladd(ladd(X[p].get(0, {}), X[p].get(1, {})), X[p].get(2, {})), L((-1, 0, 0, 0, 1)))
        if v: out[1053+p] = v
    for f in range(13):
        for k in range(27):
            v = {}
            for t in range(3): v = ladd(v, X[3*f+t].get(k, {}))
            if k < 3: v = ladd(v, L((-1, 0, 0, 0, 1)))
            if v: out[1092+27*f+k] = v
    for p in range(39):
        for q in range(p+1, 39):
            v = ladd(trace_form(X[p], X[q]), L(kneg(GAMMA(p, q))))
            if v: out[pairrow(p, q)] = v
    return out

# ---------------------------------------------------------------- covectors
def h_covectors():
    rec = D['H_phase_recipe']; H = {}
    for c in D['H_base_covectors']:
        i, a = c['component']; lam = {}
        for r, *cf in c['entries']:
            v = coef(*cf)
            if c['conjugated']: v = kconj(v)
            if c['negated']: v = kneg(v)
            lam[r] = ladd(lam.get(r, {}), L(v))
        # E_i = tau_i(U)^-1 on planes (3+2j, 4+2j): chi_E = U^-d, chi'_E = U^d
        Emat = {}                         # (n, m) -> Laurent, the 27x27 block, identity elsewhere
        for t in range(3): Emat[t, t] = L(K1)
        for j in range(12):
            dd = rec['native_plane_exponents'][i][j]
            chi = u_mono([-x for x in dd]); chip = u_mono(dd)
            half = (1, 0, 0, 0, 2); ihalf = (0, 1, 0, 0, 2)
            x_, y_ = 3+2*j, 4+2*j
            Emat[x_, x_] = lscale(ladd(chi, chip), half)
            Emat[x_, y_] = lscale(ladd(chip, lscale(chi, (-1, 0, 0, 0, 1))), ihalf)
            Emat[y_, x_] = lscale(ladd(chi, lscale(chip, (-1, 0, 0, 0, 1))), ihalf)
            Emat[y_, y_] = lscale(ladd(chi, chip), half)
        cinv = u_mono([-x for x in rec['character_exponents'][i][a]])
        out = {}
        for r, p in lam.items():
            if r < 1053 or 1092 <= r < 1443:
                blk = r - (r % 27 if r < 1053 else (r - 1092) % 27); n = r - blk
                for m in range(27):
                    if (n, m) in Emat:
                        out[blk + m] = ladd(out.get(blk + m, {}), lmul(p, Emat[n, m]))
            else:
                out[r] = ladd(out.get(r, {}), p)
        H[i, a] = {r: lspec(lmul(p, cinv)) for r, p in out.items() if p and lspec(lmul(p, cinv))}
    return H
def t_covectors():
    T = {}
    for c in D['T_phase_covectors']:
        i, a = c['component']; lam = {}
        for r, A, B, C, Dd, N, *m in c['entries']:
            lam[r] = ladd(lam.get(r, {}), lscale(u_mono(m), coef(A, B, C, Dd, N)))
        T[i, a] = {r: lspec(p) for r, p in lam.items() if lspec(p)}
    return T

# ---------------------------------------------------------------- V_xi on the 24 real basis vectors of u
IND = [(0, 1, 2), (6, 9, 12), (8, 11, 14), (7, 10, 13)]
PAIRS = [(0, 3), (1, 5), (2, 4), (6, 15), (9, 21), (12, 18), (8, 16), (11, 22), (14, 19), (7, 17), (10, 23), (13, 20)]
def V_matrices():
    tg = D['tangent']
    assert [list(t) for t in IND] == tg['independent_indices']
    assert [list(p) for p in PAIRS] == [o['eigen_pair'] for o in tg['operator_pairs']]
    sub = tg['real_coordinate_substitution']
    assert all(sub[IND[b][k]] == [b, k, False, 1] for b in range(4) for k in range(3))
    assert all(sub[jp] == [sub[j][0], sub[j][1], True, -1] for j, jp in PAIRS)
    Vm = {}
    for o in tg['operator_pairs']:
        R = {}; S = {}
        for r, s, *c in o['real_difference']: R[r, s] = kadd(R.get((r, s), K0), coef(*c))
        for r, s, *c in o['imaginary_sum']: S[r, s] = kadd(S.get((r, s), K0), coef(*c))
        a, b = o['eigen_pair']; Ka = {}; Kb = {}
        for key in set(R) | set(S):
            rr = R.get(key, K0); ss = kmul(S.get(key, K0), KI)
            Ka[key] = kmul(kadd(rr, kneg(ss)), (1, 0, 0, 0, 2))
            Kb[key] = kmul(kadd(rr, ss), (-1, 0, 0, 0, 2))
        Vm[a] = Ka; Vm[b] = Kb
    return Vm
def apply_mat(M, X):        # M: (r,s) -> K ; X: coord -> Laurent
    out = {}
    for (r, s), k in M.items():
        if s in X: out[r] = ladd(out.get(r, {}), lscale(X[s], k))
    return {r: p for r, p in out.items() if p}
def velocity_basis(P, Vm):
    """v[alpha][p] = (V_xi e_alpha)_p, alpha = 6b + 2k + (0 Re, 1 Im) for u_{bk}; also max imaginary defect."""
    VP = {j: [apply_mat(Vm[j], P[p]) for p in range(39)] for j in range(24)}
    vs = []; imag_ok = True
    for b in range(4):
        for k in range(3):
            j = IND[b][k]; jp = [q for (a_, q) in PAIRS if a_ == j]
            jp = jp[0] if jp else [a_ for (a_, q) in PAIRS if q == j][0]
            for val in (K1, KI):                  # u_bk = 1, then u_bk = i
                ut = {j: val, jp: kneg(kconj(val))}
                v = []
                for p in range(39):
                    s = {}
                    for jj, c in ut.items():
                        for kk, poly in VP[jj][p].items():
                            s[kk] = ladd(s.get(kk, {}), lscale(poly, c))
                    re = {}
                    for kk, poly in s.items():
                        if lim(poly): imag_ok = False
                        r_ = lre(poly)
                        if r_: re[kk] = r_
                    v.append(re)
                vs.append(v)
    return vs, imag_ok, VP

# ---------------------------------------------------------------- the targets H, T (Appendix C, written out)
HT_TEXT = {
 'H0': [r'u30 ~u20 + u11 ~u31 - u22 ~u12', r'u10 ~u30 + u21 ~u11 - u32 ~u22', r'- u20 ~u10 - u31 ~u21 - u12 ~u32'],
 'H1': [r'- u01 ~u30 + ~u02 ~u20 + u31 u22', r'- u00 ~u31 - ~u01 ~u21 - u20 u32', r'~u00 ~u22 + u02 ~u32 - u30 u21'],
 'H2': [r'~u00 ~u30 - u02 ~u10 + u32 u11', r'u01 ~u11 - ~u02 ~u31 + u30 u12', r'- u00 ~u12 - ~u01 ~u32 - u10 u31'],
 'H3': [r'- u00 ~u20 - ~u01 ~u10 - u21 u12', r'- ~u00 ~u11 + u02 ~u21 + u10 u22', r'u01 ~u22 + ~u02 ~u12 - u20 u11'],
 'T0': [r'- u10 ~u32 + u30 ~u22 - u20 ~u12', r'- u10 ~u21 - u30 ~u11 - u20 ~u31', r'- u11 ~u22 + u31 ~u12 + u21 ~u32'],
 'T1': [r'u00 u21 + ~u01 u31 - ~u30 ~u22', r'- u00 ~u32 - ~u02 ~u22 - u20 u31', r'u01 u20 - ~u02 u30 + ~u32 ~u21'],
 'T2': [r'u00 u32 + ~u01 u12 + ~u30 ~u11', r'- u00 ~u10 + ~u02 ~u30 - u31 u12', r'u01 u31 - ~u02 u11 + ~u10 ~u32'],
 'T3': [r'u00 u10 + ~u01 u20 - ~u22 ~u11', r'- u00 ~u21 + ~u02 ~u11 - u20 u12', r'u01 u12 + ~u02 u22 + ~u10 ~u21'],
}
def parse_component(s):
    toks = s.replace('+', ' + ').replace('-', ' - ').split()
    terms = []; sign = 1; fac = []
    for tk in toks:
        if tk in '+-':
            if fac: terms.append((sign, fac)); fac = []
            sign = 1 if tk == '+' else -1
        else:
            bar = tk.startswith('~'); tk = tk.lstrip('~'); fac.append((int(tk[1]), int(tk[2]), bar))
    terms.append((sign, fac)); assert all(len(f) == 2 for _, f in terms); return terms
def target_coeffs(s):
    """Coefficients of the monomials r_alpha r_beta (alpha <= beta) of a written-out component."""
    out = {}
    for sign, ((i1, a1, b1), (i2, a2, b2)) in parse_component(s):
        f1 = [(6*i1+2*a1, K1), (6*i1+2*a1+1, (0, -1 if b1 else 1, 0, 0, 1))]
        f2 = [(6*i2+2*a2, K1), (6*i2+2*a2+1, (0, -1 if b2 else 1, 0, 0, 1))]
        for x, cx in f1:
            for y, cy in f2:
                key = (min(x, y), max(x, y)); c = kmul(kmul(cx, cy), (sign, 0, 0, 0, 1))
                out[key] = kadd(out.get(key, K0), c)
    return {k: v for k, v in out.items() if not kiszero(v)}

# ---------------------------------------------------------------- the checks
def contract_DF(lam, P):
    """lam DF_P as dict packed column -> Laurent (nonzero only)."""
    out = {}
    def add(col, p):
        if p: out[col] = ladd(out.get(col, {}), p)
    for r, l in lam.items():
        if r < 1053:
            p, k = divmod(r, 27)
            # d/dX_p of (X_p o X_p - X_p)_k = 2 (P_p o e_c)_k - delta_kc
            for (a, b, v) in JTk.get(k, ()):
                if a in P[p]: add(27*p + b, lmul(l, lscale(P[p][a], kmul(v, (2, 0, 0, 0, 1)))))
            add(27*p + k, lscale(l, (-1, 0, 0, 0, 1)))
        elif r < 1092:
            p = r - 1053
            for c in range(3): add(27*p + c, l)
        elif r < 1443:
            f, k = divmod(r - 1092, 27)
            for t in range(3): add(27*(3*f+t) + k, l)
        else:
            p, q = PAIR_OF[r]
            # d/dX_p tr(X_p o X_q) = tr(e_c o P_q)
            for (x, y) in ((p, q), (q, p)):
                for c in range(27):
                    g = trace_coeff(c, P[y])
                    if g: add(27*x + c, lmul(l, g))
    return {c: p for c, p in out.items() if p}
PAIR_OF = {pairrow(p, q): (p, q) for p in range(39) for q in range(p+1, 39)}
_TRC = {}
def trace_coeff(c, X):     # tr(e_c o X) as Laurent
    key = c
    if key not in _TRC:
        _TRC[key] = [(b, v) for (k, a, b), v in JT.items() if k < 3 and a == c]
    r = {}
    for b, v in _TRC[key]:
        if b in X: r = ladd(r, lscale(X[b], v))
    return r

def quad_coeffs(lam, vs):
    """C[alpha,beta] (alpha<=beta): coefficient of r_alpha r_beta in lam Q(V_xi u)."""
    C = {}
    def acc(key, p):
        if p: C[key] = ladd(C.get(key, {}), p)
    idrows = {}
    for r, l in lam.items():
        if r < 1053: idrows.setdefault(r // 27, {})[r % 27] = l
    for p, lk in idrows.items():
        M = {}                       # (a,b) -> Laurent, symmetric
        for k, l in lk.items():
            for (a, b, v) in JTk.get(k, ()):
                M[a, b] = ladd(M.get((a, b), {}), lscale(l, v))
        M = {ab: x for ab, x in M.items() if x}
        W = []
        for al in range(24):
            va = vs[al][p]; w = {}
            for (a, b), x in M.items():
                if a in va: w[b] = ladd(w.get(b, {}), lmul(x, va[a]))
            W.append(w)
        for al in range(24):
            for be in range(al, 24):
                vb = vs[be][p]; s = {}
                for b, x in W[al].items():
                    if b in vb: lacc(s, x, vb[b])
                acc((al, be), s if al == be else lscale(s, (2, 0, 0, 0, 1)))
    for r, l in lam.items():
        if r >= 1443:
            p, q = PAIR_OF[r]
            for al in range(24):
                for be in range(al, 24):
                    s = trace_form(vs[al][p], vs[be][q])
                    if al != be: s = ladd(s, trace_form(vs[be][p], vs[al][q]))
                    acc((al, be), lmul(l, s))
    return {k: v for k, v in C.items() if v}

def support_defect(lam, pts):
    bad = []
    idrows = {}
    for r, l in lam.items():
        if r < 1053: idrows.setdefault(r // 27, {})[r % 27] = l
    for p, lk in idrows.items():
        if p in pts: continue
        M = {}
        for k, l in lk.items():
            for (a, b, v) in JTk.get(k, ()):
                M[a, b] = ladd(M.get((a, b), {}), lscale(l, v))
        if any(M.values()): bad.append(('point', p))
    for r in lam:
        if r >= 1443:
            p, q = PAIR_OF[r]
            if p not in pts or q not in pts: bad.append(('pair', p, q))
    return bad

# ---------------------------------------------------------------- Section 6.3: V_xi u as a sum of family velocities
def c(n, d=1): return kn(n, 0, 0, 0, d)
def xm(*e): return xi_mono(e)
def neg(p): return lscale(p, c(-1))
W = lambda k: L(opow(k))
Z3 = [[{} for _ in range(3)] for _ in range(3)]
def mat(rows): return [[(e if isinstance(e, dict) else (L(c(e)) if e else {})) for e in r] for r in rows]
def mmul(A, B):
    C = [[{} for _ in range(3)] for _ in range(3)]
    for i in range(3):
        for j in range(3):
            for k in range(3):
                if A[i][k] and B[k][j]: lacc(C[i][j], A[i][k], B[k][j])
    return C
def madd(A, B): return [[ladd(A[i][j], B[i][j]) for j in range(3)] for i in range(3)]
def mscale(A, p): return [[lmul(A[i][j], p) for j in range(3)] for i in range(3)]
def diag(a, b, cc): return mat([[a, 0, 0], [0, b, 0], [0, 0, cc]])
I3 = diag(1, 1, 1)
LM = [mat([[0, 0, xm(0, 0, 1, 0)], [0, neg(xm(0, 0, 0, -1)), 0], [1, 0, 0]]),
      mat([[0, 0, neg(xm(-1, 0, 0, 0))], [0, xm(0, 1, 0, 0), 0], [-1, 0, 0]]),
      mat([[-1, 0, 0], [0, 0, neg(xm(0, 0, -1, -1))], [0, xm(0, 0, -1, 0), 0]]),
      mat([[0, -1, 0], [-1, 0, 0], [0, 0, -1]])]
RM = [diag(1, xm(0, 0, 0, 1), xm(0, 0, -1, 0)), diag(-1, neg(xm(0, -1, 0, 0)), neg(xm(1, 0, 0, 0))),
      diag(-1, neg(xm(0, 0, 1, 0)), neg(xm(0, 0, 1, 1))), I3]
DT = [(xm(1, 1, -1, 0), xm(-1, 0, 0, -1), xm(0, -1, 1, 1)), (xm(1, 0, 1, 1), neg(xm(0, 1, -1, 0)), neg(xm(-1, -1, 0, -1))),
      (neg(xm(0, 0, 1, 1)), xm(0, 0, -1, 0), neg(xm(0, 0, 0, -1))), (neg(xm(-1, 0, 0, 0)), neg(xm(0, -1, 0, 0)), xm(1, 1, 0, 0))]
GM = [I3, diag(1, -1, -1), diag(1, xm(1, 0, 0, 0), xm(0, -1, 0, 0)), diag(1, xm(0, 0, -1, 0), neg(xm(0, 0, 0, 1)))]
def zeta(i, a, b):
    k = [a*b, 0, a*b - a, (b if a == 0 else (-b if a == 1 else 0))][i]; return W(k)
def s(i, a, b):
    if i == 0: return (-b % 3, -a % 3)
    if i == 1: return ((-b) % 3, 2) if a == 2 else (b, a)
    if i == 2: return (a, b)
    return ((-b) % 3, a) if a == 2 else (b, a)
Xs = mat([[0, 0, 1], [1, 0, 0], [0, 1, 0]])          # X e_k = e_{k+1}
Zc = diag(1, W(1), W(2))
def mpow(A, n):
    R = I3
    for _ in range(n % 3): R = mmul(R, A)
    return R
def E(ab, d):
    cl, sh = ab; Dm = diag(*d) if sh == 2 else I3
    return mmul(mmul(mpow(Xs, sh), Dm), mpow(Zc, cl))
def Umat(i, a, b): return mscale(mmul(mmul(LM[i], E(s(i, a, b), DT[i])), RM[i]), zeta(i, a, b))
PHI = [list(range(13)), [0, 4, 5, 6, 1, 2, 3, 7, 8, 9, 10, 11, 12], [0, 7, 9, 8, 1, 10, 4, 3, 12, 5, 2, 11, 6], [0, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9]]
def psi(i, f, t):
    if i == 0: return t if f < 4 else -t % 3
    if i == 1:
        if f == 0: return t
        if f <= 6: return -t % 3
        if f <= 9: return (t - (f - 7)) % 3
        return ((f - 10) - t) % 3
    if i == 2: return -t % 3 if f in (5, 8, 11) else t
    return t if f in (0, 7, 8, 9) else -t % 3
TT = {1: ((6, 4, 2), ('-+++ +---', '--++ +--+', '++-+ -+++')),
      2: ((2, 6, 4), ('-+-- +++-', '+++- ---+', '-+-+ ++++')),
      3: ((4, 2, 6), ('++-- +--+', '--+- -+++', '-+++ -+++'))}
def oct_from_cz(cpoly, z):  # complex Laurent c, z1..z3 -> 8 real-coordinate Laurents
    out = [lre(cpoly), lim(cpoly), lre(z[0]), lim(z[0]), lre(z[1]), lim(z[1]), lre(z[2]), neg(lim(z[2]))]
    return out
def tits_to_vec(a, B, i):
    """(a,B) in the coordinates of quartet i -> 27 real coordinates of iota_i(a,B)."""
    cc = {1: a[1][2], 2: a[2][0], 3: a[0][1]}
    xs = {}
    for r in (1, 2, 3):
        z = [neg(lconj(B[r-1][k])) for k in range(3)]       # z_r = -conj(b_r)
        xs[r] = oct_from_cz(cc[r], z)
    if i > 0:
        mu, sg = TT[i]
        for r in (1, 2, 3):
            signs = sg[r-1].replace(' ', ''); new = [{} for _ in range(8)]
            for m in range(8):
                new[m ^ mu[r-1]] = ladd(new[m ^ mu[r-1]], lscale(xs[r][m], c(1 if signs[m] == '+' else -1)))
            xs[r] = new
    v = {}
    for t in range(3):
        if a[t][t]: v[t] = lre(a[t][t])
    o = {3: xs[3], 11: [(p if m == 0 else neg(p)) for m, p in enumerate(xs[2])], 19: xs[1]}
    for base, x in o.items():
        for m in range(8):
            if x[m]: v[base + m] = x[m]
    return {k: lspec(p) for k, p in v.items() if lspec(p)}
def lit_q(f, t):            # q_{f,t} in H_3(C), frames 0..3 of eq:literal
    if f == 0: return [[L(c(1)) if (r == s_ == t) else {} for s_ in range(3)] for r in range(3)]
    v = [L(K1), W(t), W(f - 1 - t)]
    return [[lscale(lmul(v[r], lconj(v[s_])), c(1, 3)) for s_ in range(3)] for r in range(3)]
def velocity_check(P, vs):
    # --- check eq:anchorbinding: P(xi)_{phi_i(f), psi_i(f,t)} = iota_i hatP_i(f,t), identically
    ok = True
    for i in range(4):
        G = GM[i]; Gs = [[lconj(G[s_][r]) for s_ in range(3)] for r in range(3)]
        for f in range(13):
            for t in range(3):
                if f < 4: a = mmul(mmul(G, lit_q(f, t)), Gs); B = Z3
                else:
                    ab = divmod(f - 4, 3); a = diag(L(c(1, 3)), L(c(1, 3)), L(c(1, 3)))
                    B = mscale(Umat(i, *ab), lscale(W(t), c(1, 3)))
                got = tits_to_vec(a, B, i); want = P[3*PHI[i][f] + psi(i, f, t)]
                if got != want: ok = False; print('binding mismatch', i, f, t)
    print('eq:anchorbinding reproduces P(xi) for all four quartets, 39 points each:', ok)
    # --- velocities
    TAB = [[('y0', 'x0-'), ('y2', 'x2'), ('y1', 'x1-')], [('y0', 'x0'), ('y1', 'x1-'), ('y2', 'x2')],
           [('y1', 'x1-'), ('y2', 'x2'), ('y0', 'x0-')], [('y2', 'x2'), ('y0', 'x0-'), ('y1', 'x1')]]
    PAIRS3 = [(0, 1), (0, 2), (1, 2)]
    def Nmat(i, xv, yv):        # x, y: lists of 3 rationals (Re, Im of u_i)
        N = [[{} for _ in range(3)] for _ in range(3)]
        for (r, s_), (im_, re_) in zip(PAIRS3, TAB[i]):
            val = lambda spec: (xv if spec[0] == 'x' else yv)[int(spec[1])] * (-1 if spec.endswith('-') else 1)
            ci, cr = val(im_), val(re_)
            # ci * sqrt-1 (E_rs + E_sr) + cr (E_rs - E_sr)
            N[r][s_] = ladd(N[r][s_], L(kn(cr, ci, 0, 0, 1))); N[s_][r] = ladd(N[s_][r], L(kn(-cr, ci, 0, 0, 1)))
        return N
    allok = True; nzero = 0
    for al in range(24):
        b, rem = divmod(al, 6); k, im = divmod(rem, 2)
        xv = [0, 0, 0]; yv = [0, 0, 0]
        (yv if im else xv)[k] = 1
        N = Nmat(b, xv, yv); vel = [dict() for _ in range(39)]
        for f in range(4, 13):
            ab = divmod(f - 4, 3); Um = Umat(b, *ab)
            comm = madd(mmul(N, Um), mscale(mmul(Um, N), L(c(-1))))
            for t in range(3):
                Bv = mscale(comm, lscale(W(t), c(1, 3)))
                vel[3*PHI[b][f] + psi(b, f, t)] = tits_to_vec(Z3, Bv, b)
        same = all(vel[p] == vs[al][p] for p in range(39))
        if not same:
            allok = False
            diffp = [p for p in range(39) if vel[p] != vs[al][p]]
            print(f'alpha {al} (u_{b}{k} {"Im" if im else "Re"}): differs at points {diffp[:8]}')
    print('V_xi e_alpha equals the transported family velocity with g=h=exp(t n_i(u_i)), all 24 alpha:', allok)
    return ok and allok

def main():
    t0 = time.time()
    mode = 'at xi = xi_* only (U = 1)' if BASE_ONLY else 'identically in xi (Laurent coefficients in U)'
    print('input sha256:', hashlib.sha256(open(sys.argv[1], 'rb').read()).hexdigest())
    print('mode:', mode)
    P = config()
    selfconj = all(lconj(p) == p for X in P for p in X.values())
    print('P(xi) real on the torus (each coordinate self-conjugate):', selfconj)
    Fp = residual(P)
    print('F(P(xi)) = 0 exactly:', not Fp, f'[{time.time()-t0:.0f}s]')
    Vm = V_matrices()
    vs, imag_ok, VP = velocity_basis(P, Vm)
    print('sum_j u~_j V^(j) P(xi) is real before taking Re (all 24 basis vectors):', imag_ok)
    # redundant seed columns: V^(j) P(xi_*) = column eigen_index[j]
    if True:
        tg = D['tangent']; den = tg['seed_columns_denominator']; ok = True
        for j in range(24):
            col = {}
            for pc, A, B, C_, Dd in tg['seed_columns'][tg['eigen_index'][j]]['entries']:
                col[pc] = kadd(col.get(pc, K0), kn(A, B, C_, Dd, den))
            got = {}
            for p in range(39):
                for kk, poly in VP[j][p].items():
                    s = K0
                    for c in poly.values(): s = kadd(s, c)       # evaluate at U = 1
                    if not kiszero(s): got[27*p+kk] = s
            ok &= (got == {k: v for k, v in col.items() if not kiszero(v)})
        print('seed columns equal V^(j) P(xi_*) for all 24 j (redundant check, at xi_* only):', ok)
    # DF V_xi u = 0 (first-order deformation)
    kerok = True
    for al in range(24):
        X = vs[al]
        lin = {}
        # DF_P v via residual polarization: F(P+v) - F(P) - Q(v) is linear; compute directly
        for p in range(39):
            pv = jordan(P[p], X[p])
            for k in set(pv) | set(X[p]):
                val = ladd(lscale(pv.get(k, {}), (2, 0, 0, 0, 1)), lscale(X[p].get(k, {}), (-1, 0, 0, 0, 1)))
                if val: kerok = False
            if ladd(ladd(X[p].get(0, {}), X[p].get(1, {})), X[p].get(2, {})): kerok = False
        for f in range(13):
            for k in range(27):
                s = {}
                for t in range(3): s = ladd(s, X[3*f+t].get(k, {}))
                if s: kerok = False
        for p in range(39):
            for q in range(p+1, 39):
                if ladd(trace_form(P[p], X[q]), trace_form(X[p], P[q])): kerok = False
    print('DF_P (V_xi e_alpha) = 0 for all 24 real basis vectors:', kerok, f'[{time.time()-t0:.0f}s]')
    velok = velocity_check(P, vs)
    Hc = h_covectors(); Tc = t_covectors()
    # E_i P(xi) = P(xi_*) on quartet i, i.e. tau_i(U) P(xi_*) = P(xi)
    ANCH = D['H_phase_recipe']['anchor_frames']
    assert ANCH == [[0, 1, 2, 3], [0, 4, 5, 6], [0, 7, 8, 9], [0, 10, 11, 12]]
    allok = True
    for name, cov in (('h', Hc), ('t', Tc)):
        for (i, a), lam in sorted(cov.items()):
            eq = contract_DF(lam, P)
            C = quad_coeffs(lam, vs)
            tgt = target_coeffs(HT_TEXT[('H' if name == 'h' else 'T') + str(i)][a])
            nonconst = [k for k, v in C.items() if any(e != Z4 for e in v)]
            Cc = {k: v.get(Z4, K0) for k, v in C.items()}
            Cc = {k: v for k, v in Cc.items() if not kiszero(v)}
            repro = (not nonconst) and Cc == tgt
            line = f'{name}_{i}{a}: rows {len(lam)}, equilibrium (1053 columns) {"OK" if not eq else "FAIL %d" % len(eq)}, ' \
                   f'reproduction (300 monomials) {"OK" if repro else "FAIL"}'
            if not repro:
                diff = set(Cc) ^ set(tgt) | {k for k in set(Cc) & set(tgt) if Cc[k] != tgt[k]}
                line += f' [nonconstant {len(nonconst)}, differing {sorted(diff)[:6]}]'
            if name == 'h':
                pts = {3*f+t for f in ANCH[i] for t in range(3)}
                bad = support_defect(lam, pts)
                rowsup = all((r < 1053 and r // 27 in pts) or (1053 <= r < 1092 and r - 1053 in pts) or
                             (1092 <= r < 1443 and (r - 1092) // 27 in ANCH[i]) or
                             (r >= 1443 and set(PAIR_OF[r]) <= pts) for r in lam)
                line += f', support on quartet {i} {"OK" if not bad else "FAIL " + str(bad[:4])}' \
                        f' (row support {"OK" if rowsup else "no"})'
                allok &= not bad
            allok &= (not eq) and repro
            print(line, f'[{time.time()-t0:.0f}s]', flush=True)
    print('ALL CHECKS PASSED' if allok and not Fp and kerok and imag_ok and velok else 'SOME CHECK FAILED', '-', mode)

if __name__ == '__main__':
    main()