"""test_629 -- 2 *DF_gamma on TT is Branson's third-order intertwiner, with eigenvalues +-n(n^2 - 1).

MA Remark 6.1: at the round S^3, DF_gamma = DC_gamma, the linearised Cotton operator, and on TT
DF_gamma = Cod o (1/2)(-Delta + 2) (MA Theorem E(2)). Dualising the antisymmetric pair, *T_cd = (1/2) eps_c^ab T_abd, turns
the Cotton tensor into the Cotton-York tensor Y, which is symmetric, trace-free and divergence-free for every metric, with
Y(gamma) = 0; so *DF_gamma = DY_gamma maps all of Sym^2 into TT, and on TT *Cod is the curl, (curl h)_cd = eps_c^ab nabla_a h_bd.
  P1 [floating point, exact representation theory] On S^3 = SU(2) with a left-invariant orthonormal frame
     (nabla_a e_b = eps_abk e_k, unit curvature), functions decompose as V_l (x) V_l^*; the frame fields act on V_l as
     -2i J_a. For 2l = 0, ..., 8 the operators nabla_a, div, curl and nabla^* nabla are assembled on V_l (x) Sym^2_0 and
     restricted to TT = ker div. On TT: curl preserves TT; curl^2 = nabla^* nabla + 3 as operators; every eigenvalue of
     curl is +-n with n an integer >= 3; and for each n <= 7 the multiplicities, weighted by dim V_l^* = 2l + 1, are
     n^2 - 4 for +n and n^2 - 4 for -n (the two chiral summands; total 2(n^2 - 4), test_627's multiplicity). Since
     -Delta + 2 = nabla^* nabla + 2 = n^2 - 1 on level n, 2 *DF_gamma = curl o (-Delta + 2) has eigenvalues +-n(n^2 - 1), and
     |DF_gamma h|^2 = (1/2) n^2 (n^2 - 1)^2 |h|^2 because |C|^2 = 2 |*C|^2: the 1/sqrt 2 of MA Theorem E(4).
  P2 [exact] Conformal covariance. For the Moebius dilation Phi(y) = 2y of the stereographic chart,
     Phi^* gamma = Omega^2 gamma with Omega = 2(1 + |y|^2)/(1 + 4|y|^2), and DC_gamma(Omega^-2 Phi^* H) = Phi^* DC_gamma(H) at two
     rational points for a general polynomial H (test_625's exact dual-number engine). With eps_{Omega^2 g} = Omega^3 eps_g
     this is DY_gamma(Omega^-2 Phi^* H) = Omega Phi^* DY_gamma(H): weights 0 and 3, Branson's third-order intertwiner.
So 0 -> so(4,1) -> X(S^3) --K--> Sym^2_0 --2*DF_gamma--> TT -> 0 is exact: at Sym^2_0 by MA Theorem E(1)-(3), and onto TT
because 2*DF_gamma is invertible on TT (P1).
History: proposed by Gemini (12 November 2026, extensions note section 1.1).
"""
import sys, os
import numpy as np
import sympy as sp
from sympy import Rational as Q
FAIL = []
def check(name, ok, detail=""):
    ok = bool(ok)
    if not ok: FAIL.append(name)
    print(f"[{'PASS' if ok else 'FAIL'}] {name}" + (f"   ({detail})" if detail else ""), flush=True)

# ---------- P1
def spin(l2):
    l = l2/2; m = np.arange(l, -l - 1, -1); d = len(m)
    Jp = np.zeros((d, d), complex)
    for k in range(1, d):
        Jp[k - 1, k] = np.sqrt(l*(l + 1) - m[k]*(m[k] + 1))
    Jx = (Jp + Jp.conj().T)/2; Jy = (Jp - Jp.conj().T)/(2j); Jz = np.diag(m).astype(complex)
    return [Jx, Jy, Jz]
eps3 = np.zeros((3, 3, 3))
for a, b, c in [(0, 1, 2), (1, 2, 0), (2, 0, 1)]:
    eps3[a, b, c] = 1; eps3[b, a, c] = -1
I3 = np.eye(3)
# orthonormal basis of Sym^2_0
S = []
for (i, j) in [(0, 1), (0, 2), (1, 2)]:
    M = np.zeros((3, 3)); M[i, j] = M[j, i] = 1/np.sqrt(2); S.append(M)
S.append(np.diag([1, -1, 0])/np.sqrt(2)); S.append(np.diag([1, 1, -2])/np.sqrt(6))
Bsym = np.array([M.reshape(9) for M in S]).T          # 9 x 5
counts = {}
ok_pres = ok_weitz = ok_int = True
worst = 0.0
for l2 in range(0, 9):
    J = spin(l2); dV = l2 + 1; IV = np.eye(dV)
    rho = [-2j*Ja for Ja in J]
    A = [eps3[a] for a in range(3)]                    # (A_a)_{bk} = eps_abk
    nab = [np.kron(rho[a], np.kron(I3, I3)) - np.kron(IV, np.kron(A[a], I3)) - np.kron(IV, np.kron(I3, A[a])) for a in range(3)]
    sel = [np.kron(IV, np.kron(I3[b:b + 1, :], I3)) for b in range(3)]   # pick first index b
    Div = sum(sel[b] @ nab[b] for b in range(3))
    E = [np.array([[eps3[c, a, b] for b in range(3)] for c in range(3)]) for a in range(3)]
    Curl = sum(np.kron(IV, np.kron(E[a], I3)) @ nab[a] for a in range(3))
    Rough = -sum(nab[a] @ nab[a] for a in range(3))
    B = np.kron(IV, Bsym)                               # 9 dV x 5 dV
    DB = Div @ B
    u, s, vh = np.linalg.svd(DB)
    null = vh[np.sum(s > 1e-9):].conj().T                 # 5 dV x k
    N = B @ null
    N, _ = np.linalg.qr(N)
    CN = Curl @ N
    resid = np.linalg.norm(CN - N @ (N.conj().T @ CN))
    ok_pres &= resid < 1e-9
    Ct = N.conj().T @ CN; Rt = N.conj().T @ Rough @ N
    ok_weitz &= np.linalg.norm(Ct @ Ct - Rt - 3*np.eye(N.shape[1])) < 1e-9
    ev = np.linalg.eigvals(Ct)
    for e in ev:
        r = round(e.real)
        worst = max(worst, abs(e - r))
        ok_int &= abs(r) >= 3
        counts[r] = counts.get(r, 0) + (l2 + 1)
ok_int &= worst < 1e-9
ok_mult = all(counts.get(nn, 0) == nn*nn - 4 and counts.get(-nn, 0) == nn*nn - 4 for nn in range(3, 8))
check("P1 on TT: curl preserves TT, curl^2 = nabla*nabla + 3, eigenvalues +-n (n >= 3 integer), multiplicities n^2 - 4 each sign (n <= 7)",
      ok_pres and ok_weitz and ok_int and ok_mult,
      f"max |eig - integer| = {worst:.1e}; weighted counts " + ", ".join(f"{nn}: {counts.get(nn, 0)}/{counts.get(-nn, 0)}" for nn in range(3, 8)))

# ---------- P2: conformal covariance with test_625's engine
here = os.path.dirname(os.path.abspath(__file__))
src = open(os.path.join(here, "test_625_dF_equals_dC.py")).read()
ns = {}
exec(src[:src.index("random.seed(625)")], ns)
run, y = ns["run"], ns["y"]
import random
random.seed(629)
def rand_poly(deg=2):
    terms = [sp.Integer(1)] + list(y) + [y[i]*y[j] for i in range(3) for j in range(i, 3)]
    return sum(random.randint(-3, 3)*t for t in terms)
lam = 2
r2 = sum(v**2 for v in y)
Om = lam*(1 + r2)/(1 + lam**2*r2)
H = sp.Matrix(3, 3, lambda i, j: 0)
for i in range(3):
    for j in range(i, 3):
        H[i, j] = H[j, i] = rand_poly()
scaled = {y[k]: lam*y[k] for k in range(3)}
Hpull = sp.Matrix(3, 3, lambda i, j: sp.cancel(Om**-2*lam**2*H[i, j].subs(scaled, simultaneous=True)))
# sanity: Phi^* gamma = Omega^2 gamma
gam = ns["gamma_stereo"]()
pull_gam = sp.Matrix(3, 3, lambda i, j: lam**2*gam[i, j].subs(scaled, simultaneous=True))
ok_om = sp.simplify(pull_gam - Om**2*gam) == sp.zeros(3, 3)
ok2 = ok_om
pts = [{y[0]: Q(1, 2), y[1]: Q(1, 3), y[2]: Q(1, 5)}, {y[0]: Q(-2, 3), y[1]: Q(1, 4), y[2]: Q(3, 7)}]
nonzero = False
for pt in pts:
    _, CL, _, _, _ = run(Hpull, pt)
    ptl = {k: lam*v for k, v in pt.items()}
    _, CR, _, _, _ = run(H, ptl)
    ok2 &= all(sp.simplify(CL[k] - lam**3*CR[k]) == 0 for k in CL)
    nonzero |= any(CL[k] != 0 for k in CL)
check("P2 DC_gamma(Omega^-2 Phi^* H) = Phi^* DC_gamma(H) for the Moebius dilation y -> 2y (exact, two points)", ok2 and nonzero)
print("\nALL PASS" if not FAIL else f"\nFAILURES: {FAIL}")
sys.exit(1 if FAIL else 0)