#!/usr/bin/env python3
"""
test_331 -- A LOWER-DEGREE POLYNOMIAL FORM OF THE REALIZABILITY OPERATOR: sqrt(det M) F.M, QUADRATIC IN THE CURVATURE
            AND LINEAR IN ITS DERIVATIVE (THE BIDEGREE OF RIVERTZ'S SIX EQUALITIES)

For MA section 8.1 (prior-art sweep, 29 Oct 2026: Rivertz 1999 / Agaoka-Hashinaga 2022 characterize realizability by
six equalities quadratic in R and linear in nabla R; the Thomas form E of Proposition 8.1 has bidegree (4,1)).

With phi = sqrt(det M), K~ = phi M^{-1}, P = adj M (endomorphisms), K~ M = phi I. Differentiating,
(nabla_a K~) M = (d_a phi) I - K~ nabla_a M, so for F_ab^c = nabla_a K~_b^c - nabla_b K~_a^c,

    phi (F M)_ab^d = (1/2)(d_a det M delta_b^d - d_b det M delta_a^d) - (P nabla_a M)_b^d + (P nabla_b M)_a^d  =: R_ab^d,

and d_a det M = tr(P nabla_a M) (Jacobi). R is polynomial: degree 2 in the curvature, 1 in its derivative.

PRE-REGISTERED:
  P1 On a generic (non-realizable) metric with det M > 0 at the point, phi (F M) = R exactly, all nine components,
         and d det M = tr(P nabla M); F is not zero there (the check is not vacuous).
  P2 On a realizable metric (the round S^3, stereographic chart, at the origin) M = I and R = 0.
  P3 At the audit's rank-2 point (g = dr^2 + (1+r+r^3)^2 dx^2 + (1+r+r^2)^2 dy^2, r = 0): R_rx^x = -6 and
         R_ry^y = 6. So R, like E, is a genuine condition at a rank-2 point, but it is a different polynomial extension
         from E there (E has the single component E_ryy = 24).
  P4 On the realizable rank-one cone over a Clifford torus, R = 0 at the point (P = 0 and det M = 0 identically).
  P5 (added after the first run, which passed P1-P4: on the round S^3, M = I everywhere and P2 holds trivially.)
         On a non-umbilic realizable graph x_4 = f, f = (x0^2 + 2 x1^2 + 3 x2^2)/2 + x0 x1 x2/3 + x1^3/5, with
         g = delta + df df: at the origin M = diag(6, 3, 2), nabla M != 0, and F = 0 and R = 0 (all components).
  P6 (added with P5) The identities of R, at the generic metric of P1: the trace R_ab^b = 0 for each a (the
         contracted Bianchi identity, div M = 0, gives ((div M) P)_a); and the cyclic identity
         sum_cyc(a,b,c) R_ab^d P_dc = 0 (lowered with g), inherited from Cod of a symmetric tensor. So the nine components
         carry six after the three trace identities, and five after the cyclic one too.

METHOD. All quantities are exact jets at the origin of a polynomial metric, by truncated power series (the 3-jet of g
fixes M, nabla M and Gamma there). Kt = sqrt(det M) M^{-1} is differentiated directly along coordinate lines, not
through the identity under test. A first version that built the full symbolic curvature of rational metrics was
stopped after ten minutes with no check run; the checks themselves are unchanged.

BUNDLE EDIT (29 October, before the bundle's final run; claims unchanged): sympy's nsimplify removed from the
comparison helper. On a number with no free symbols nsimplify evaluates numerically and can return a different
algebraic number agreeing to about 15 digits; the compared values here are exact. The check now compares exact
expressions directly. The archive's [M] labels were removed from the claims.
"""
import sys
import sympy as sp

FAILURES = []


def check(name, ok, detail=""):
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}" + (f"  -- {detail}" if detail else ""))
    if not ok:
        FAILURES.append(name)


X = sp.symbols("u0:3")


def trunc(e, deg):
    """Keep the monomials of total degree <= deg in X (e polynomial in X)."""
    e = sp.expand(e)
    if e == 0:
        return sp.Integer(0)
    return sum((c * sp.prod([x**k for x, k in zip(X, m)]) for m, c in sp.Poly(e, *X).terms() if sum(m) <= deg),
               sp.Integer(0))


def tmat(Mx, deg):
    return Mx.applyfunc(lambda e: trunc(e, deg))


def geometry_jet(g):
    """Exact 1-jet of M and 0-jet of Gamma at the origin, from the 3-jet of a polynomial metric g (truncated power
    series; no rational functions are formed). Returns M0, [d_a M (0)], Gamma0[i][a][e] = Gamma^i_ae(0)."""
    g = tmat(g, 3)
    g0 = g.subs({x: 0 for x in X})
    g0i = g0.inv()
    A = g - g0
    gi, term = g0i, g0i
    for _ in range(2):
        term = tmat(-g0i * A * term, 2)
        gi = gi + term
    gi = tmat(gi, 2)
    dg = [[[sp.diff(g[b, c], X[a]) for c in range(3)] for b in range(3)] for a in range(3)]   # dg[a][b][c] = d_a g_bc
    Gam = [[[trunc(sum(gi[i, e] * (dg[c][e][b] + dg[b][e][c] - dg[e][b][c]) for e in range(3)) / 2, 2)
             for c in range(3)] for b in range(3)] for i in range(3)]
    Riem = lambda a, b, c, d: trunc(sp.diff(Gam[a][d][b], X[c]) - sp.diff(Gam[a][c][b], X[d])
                                    + sum(Gam[a][c][e] * Gam[e][d][b] - Gam[a][d][e] * Gam[e][c][b] for e in range(3)), 1)
    Ric = sp.Matrix(3, 3, lambda b, d: sum(Riem(a, b, a, d) for a in range(3)))
    Rs = trunc(sum(gi[b, d] * Ric[b, d] for b in range(3) for d in range(3)), 1)
    G = tmat(Ric - Rs * g / 2, 1)
    M = tmat(-gi * G, 1)                     # endomorphism M^i_j (row i, column j)
    o = {x: 0 for x in X}
    M0 = M.subs(o)
    M1 = [M.diff(X[a]).subs(o) for a in range(3)]
    G0 = [[[Gam[i][a][e].subs(o) for e in range(3)] for a in range(3)] for i in range(3)]
    return M0, M1, G0


def cov(T0, dT, G0, a):
    """(nabla_a T)^i_j at the origin from T(0) and d_a T(0)."""
    return sp.Matrix(3, 3, lambda i, j: dT[i, j] + sum(G0[i][a][e] * T0[e, j] - G0[e][a][j] * T0[i, e] for e in range(3)))


def R_form(M0, M1, G0):
    """R_ab^d = (1/2)(d_a det M delta_b^d - d_b det M delta_a^d) - ((nabla_a M) P)^d_b + ((nabla_b M) P)^d_a."""
    P = M0.adjugate()
    nM = [cov(M0, M1[a], G0, a) for a in range(3)]
    ddet = [(P * nM[a]).trace() for a in range(3)]
    R = {}
    for a in range(3):
        for b in range(3):
            for d in range(3):
                R[(a, b, d)] = sp.simplify(sp.Rational(1, 2) * (ddet[a] * int(b == d) - ddet[b] * int(a == d))
                                           - (nM[a] * P)[d, b] + (nM[b] * P)[d, a])
    return R, ddet, P


def F_contracted(M0, M1, G0):
    """phi (F M)_ab^d, F_ab^c = nabla_a Kt^c_b - nabla_b Kt^c_a, with Kt = sqrt(det M) M^{-1} differentiated directly
    (not through the identity under test): along each coordinate line M(t) = M0 + t d_a M(0) carries the exact first
    derivative, and Kt(t) = sqrt(det M(t)) M(t)^{-1} is differentiated in t."""
    tt = sp.symbols("tt")
    phi = sp.sqrt(M0.det())
    Kt0 = phi * M0.inv()
    nK = []
    for a in range(3):
        Mt = M0 + tt * M1[a]
        dK = (sp.sqrt(Mt.det()) * Mt.inv()).diff(tt).subs(tt, 0)
        nK.append(cov(Kt0, dK, G0, a))
    out, Fvals = {}, []
    for a in range(3):
        for b in range(3):
            Fab = sp.Matrix(3, 1, lambda c, _: nK[a][c, b] - nK[b][c, a])
            Fvals += list(Fab)
            FM = M0 * Fab                      # F_ab^c M_c^d  (M self-adjoint: M_c^d = M^d_c)
            for d in range(3):
                out[(a, b, d)] = sp.simplify(phi * FM[d])
    return out, Fvals


def agree(A, B):
    return all(sp.simplify(A[k] - B[k]) == 0 for k in A)


# ---------------------------------------------------------------- P1
x0, x1, x2 = X
rr = x0**2 + x1**2 + x2**2
round3 = sp.eye(3) * (1 - rr / 2)          # 3-jet of the round S^3 in stereographic coordinates at the origin
print("P1  phi (F M) = R on a generic metric with det M > 0")
g1 = round3.copy()
g1[0, 2] = g1[2, 0] = x0 * x1 / 10 + x2**3 / 6
g1[1, 1] = g1[1, 1] + x2**2 / 7 + x0 * x1 * x2 / 5
g1[0, 1] = g1[1, 0] = x0**2 * x1 / 9
M0, M1, G0 = geometry_jet(g1)
R1, ddet1, _ = R_form(M0, M1, G0)
FM1, Fv1 = F_contracted(M0, M1, G0)
tt = sp.symbols("tt")
jac = all(sp.simplify(ddet1[a] - (M0 + tt * M1[a]).det().diff(tt).subs(tt, 0)) == 0 for a in range(3))
Fnonzero = any(sp.simplify(v) != 0 for v in Fv1)
print(f"      det M at the point = {M0.det()}")
check("det M > 0; phi (F M)_ab^d = R_ab^d (nine components); d det M = tr(P nabla M); F != 0 there",
      M0.det() > 0 and agree(FM1, R1) and jac and Fnonzero,
      f"max |F| ~ {max(abs(sp.N(v)) for v in Fv1):.3e}")

# ---------------------------------------------------------------- P2
print("P2  R = 0 on the round S^3")
M0, M1, G0 = geometry_jet(round3)
R2, _, _ = R_form(M0, M1, G0)
check("M = I and R = 0 at the point", M0 == sp.eye(3) and all(v == 0 for v in R2.values()))

# ---------------------------------------------------------------- P3
print("P3  the audit's rank-2 point")
r = X[0]
g3 = sp.diag(1, (1 + r + r**3)**2, (1 + r + r**2)**2)
M0, M1, G0 = geometry_jet(g3)
R3, _, _ = R_form(M0, M1, G0)
nz = {k: v for k, v in R3.items() if v != 0}
check("M = diag(-1,-2,0); R_rx^x = -6, R_ry^y = 6 (and antisymmetric partners)",
      M0 == sp.diag(-1, -2, 0) and R3[(0, 1, 1)] == -6 and R3[(0, 2, 2)] == 6 and R3[(1, 0, 1)] == 6
      and R3[(2, 0, 2)] == -6, f"nonzero {nz}")

# ---------------------------------------------------------------- P4
print("P4  the realizable rank-one cone")
t = sp.Rational(3, 2) + X[0]
g4 = sp.diag(1, t**2 * sp.Rational(16, 25), t**2 * sp.Rational(9, 25))
M0, M1, G0 = geometry_jet(g4)
R4, _, _ = R_form(M0, M1, G0)
check("M = diag(-1/t^2, 0, 0) and R = 0 at the point",
      M0 == sp.diag(-sp.Rational(4, 9), 0, 0) and all(v == 0 for v in R4.values()))

# ---------------------------------------------------------------- P6 (data from P1)
print("P6  the identities of R")
M0, M1, G0 = geometry_jet(g1)
g0 = g1.subs({x: 0 for x in X})
P = M0.adjugate()
trace_ok = all(sp.simplify(sum(R1[(a, b, b)] for b in range(3))) == 0 for a in range(3))
RP = lambda a, b, c: sum(R1[(a, b, d)] * (g0 * P)[d, c] for d in range(3))     # R_ab^d P_dc, P lowered with g
cyc = sp.simplify(RP(0, 1, 2) + RP(1, 2, 0) + RP(2, 0, 1))
check("R_ab^b = 0 (three); cyclic sum of R_ab^d P_dc = 0", trace_ok and cyc == 0, f"cyclic = {cyc}")

# ---------------------------------------------------------------- P5
print("P5  a non-umbilic realizable graph")
f = (x0**2 + 2 * x1**2 + 3 * x2**2) / 2 + x0 * x1 * x2 / 3 + x1**3 / 5
df = sp.Matrix([sp.diff(f, x) for x in X])
g5 = sp.eye(3) + df * df.T
M0, M1, G0 = geometry_jet(g5)
R5, _, _ = R_form(M0, M1, G0)
FM5, Fv5 = F_contracted(M0, M1, G0)
nMnz = any(cov(M0, M1[a], G0, a) != sp.zeros(3, 3) for a in range(3))
check("M = diag(6,3,2), nabla M != 0, F = 0, R = 0",
      M0 == sp.diag(6, 3, 2) and nMnz and all(sp.simplify(v) == 0 for v in Fv5) and all(v == 0 for v in R5.values()))

print()
print(f"test_331: {'ALL PASS' if not FAILURES else 'FAILURES: ' + ', '.join(FAILURES)}")
sys.exit(1 if FAILURES else 0)