#!/usr/bin/env python3
"""
test_330 -- THE SMALL DISPLAYED IDENTITIES OF THE LOCAL MATHS PAPER (MA) THAT HAD NO TEST OF THEIR OWN

For MA's verification appendix. Each is exact (symbolic or rational).

PRE-REGISTERED:
  P1 (§2.4) det of D adj_K on Sym^2 is -2 (det K)^2 for every symmetric 3x3 K (symbolic, general K).
  P2 (Theorem A(2), indefinite) for M of signature (+,-,-) with det M > 0, K~ = sqrt(det M) M^{-1} has
         adj K~ = M exactly (random rational eigenvalues, random rational rotation).
  P3 (Example 3.3, the Flamm hypersurface) with profile z^2 = 8m(rho - 2m): kappa_mer = -sqrt(2m)/(2 rho^{3/2}),
         kappa_sph = sqrt(2m)/rho^{3/2}; M has eigenvalues (2m/rho^3, -m/rho^3, -m/rho^3); det M = 2m^3/rho^9, which
         is 1/(256 m^6) at the throat rho = 2m; the radial distance from rho = 2m to 3m is finite; the induced metric
         is (1 - 2m/rho)^{-1} d rho^2 + rho^2 dOmega^2; T_1 = diag(2k, k/2, k/2) with k = kappa_sph.
  P4 (§10) n^2(n^2-1)/12 - n(n+1)/2 = n(n-3)(n+1)(n+2)/12, zero for n >= 1 only at n = 3.
  P5 (Lemma 7.3) E'/E = -W.n exactly, W = (P', Q', S')/S, n the direction cosines, for general P, Q, S.
  P6 (Proposition 7.6) for V = (1, -P, -Q, P^2+Q^2+S^2)/(2S) and Q(V) = 4(V0 V3 - V1^2 - V2^2): Q(V) = 1;
         -Q(V') = (P'^2 + Q'^2 + S'^2)/S^2; E = B(V, N)/2 with N = (1, -p, -q, p^2+q^2), Q(N) = 0;
         Q((Phi V)') = Phi'^2 - Phi^2 |W|^2.

FIRST RUN: P1, P2, P4-P6 passed. P3 failed in sympy's simplification only: without the assumption rho > 2m it could
not cancel (rho - 2m)^{3/2} against (-1/(2m - rho))^{3/2}. Re-posed with rho = 2m + u, u > 0; the claims are unchanged.

BUNDLE EDIT (29 October): the archive's [M] labels were removed from the claims; nothing else changed.
"""
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)


SYM = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]

# P1
print("P1  det D adj_K = -2 (det K)^2")
q = sp.symbols("q0:6")
K = sp.zeros(3, 3)
for v, (i, j) in zip(q, SYM):
    K[i, j] = K[j, i] = v
t = sp.symbols("t")
cols = []
for (i, j) in SYM:
    E = sp.zeros(3, 3)
    E[i, j] = E[j, i] = 1
    D = sp.diff((K + t * E).adjugate(), t).subs(t, 0)
    cols.append([D[a, b] for (a, b) in SYM])
detD = sp.Matrix(cols).T.det()
check("det D adj_K on Sym^2 = -2 (det K)^2, general symmetric K", sp.expand(detD + 2 * K.det()**2) == 0)

# P2
print("P2  Theorem A(2) on the indefinite branch")
A = sp.Matrix([[0, -sp.Rational(1, 3), sp.Rational(2, 5)], [sp.Rational(1, 3), 0, -sp.Rational(1, 7)],
               [-sp.Rational(2, 5), sp.Rational(1, 7), 0]])
R = (sp.eye(3) - A) * (sp.eye(3) + A).inv()
ok2 = True
for ev in ((4, -1, -9), (sp.Rational(9, 4), -sp.Rational(1, 4), -1)):
    M = R * sp.diag(*ev) * R.T
    Kt = sp.sqrt(M.det()) * M.inv()
    ok2 &= sp.simplify(Kt.adjugate() - M) == sp.zeros(3, 3) and M.det() > 0
check("adj(sqrt(det M) M^{-1}) = M for M of signature (+,-,-), det M > 0", ok2)

# P3
print("P3  the Flamm hypersurface")
m, u = sp.symbols("m u", positive=True)
rho = 2 * m + u                               # rho > 2m built in
rs = sp.symbols("rs", positive=True)
z = sp.sqrt(8 * m * u)
z = z.subs(u, rs - 2 * m)
dz = sp.diff(z, rs).subs(rs, rho)
ds = sp.sqrt(1 + dz**2)                       # arc length per d rho of the profile (rho, z)
grr = sp.simplify(1 + dz**2)
ok3 = sp.simplify(grr - 1 / (1 - 2 * m / rho)) == 0
# principal curvatures of the revolution hypersurface (rho w, z(rho)), normal chosen (-z', 1)/ds
kmer = sp.simplify(sp.diff(sp.diff(z, rs), rs).subs(rs, rho) / ds**3)   # profile curvature
ksph = sp.simplify(dz / (rho * ds))                   # rotational curvature
ok3 &= sp.simplify(kmer + sp.sqrt(2 * m) / (2 * rho**sp.Rational(3, 2))) == 0
ok3 &= sp.simplify(ksph - sp.sqrt(2 * m) / rho**sp.Rational(3, 2)) == 0
Kf = sp.diag(kmer, ksph, ksph)
Mf = Kf.adjugate()
ok3 &= all(sp.simplify(a - b) == 0 for a, b in zip([Mf[0, 0], Mf[1, 1], Mf[2, 2]], [2 * m / rho**3, -m / rho**3, -m / rho**3]))
dM = sp.simplify(Mf.det())
ok3 &= sp.simplify(dM - 2 * m**3 / rho**9) == 0 and sp.simplify((2 * m**3 / rs**9).subs(rs, 2 * m) - 1 / (256 * m**6)) == 0
dist = sp.integrate(1 / sp.sqrt(1 - 2 * m / rs), (rs, 2 * m, 3 * m))
ok3 &= dist.is_finite
T1 = (Kf.trace() * sp.eye(3) - Kf).subs({})
ok3 &= sp.simplify(T1 - sp.diag(2 * ksph, ksph / 2, ksph / 2)) == sp.zeros(3, 3)
check("metric, principal curvatures, M's eigenvalues, det M (1/256m^6 at the throat), finite distance, T_1", ok3,
      f"distance from 2m to 3m = {sp.simplify(dist)}")

# P4
print("P4  the count")
n = sp.symbols("n", integer=True, positive=True)
lhs = n**2 * (n**2 - 1) / 12 - n * (n + 1) / 2
ok4 = sp.simplify(lhs - n * (n - 3) * (n + 1) * (n + 2) / 12) == 0 and [k for k in range(1, 50) if lhs.subs(n, k) == 0] == [3]
check("n^2(n^2-1)/12 - n(n+1)/2 = n(n-3)(n+1)(n+2)/12, zero only at n = 3", ok4)

# P5, P6
print("P5  Lemma 7.3")
r, p, qq = sp.symbols("r p q", real=True)
Pf, Qf, Sf, Phi = [sp.Function(nm)(r) for nm in ("P", "Q", "S", "Phi")]
Eq = ((p - Pf)**2 + (qq - Qf)**2 + Sf**2) / (2 * Sf)
nvec = sp.Matrix([(p - Pf) / Eq, (qq - Qf) / Eq, 1 - Sf / Eq])
W = sp.Matrix([sp.diff(Pf, r), sp.diff(Qf, r), sp.diff(Sf, r)]) / Sf
ok5 = sp.simplify(sp.diff(Eq, r) / Eq + (W.T * nvec)[0]) == 0 and sp.simplify((nvec.T * nvec)[0] - 1) == 0
check("E'/E = -W.n, |n| = 1", ok5)
print("P6  Proposition 7.6")
Qf4 = lambda Vv: 4 * (Vv[0] * Vv[3] - Vv[1]**2 - Vv[2]**2)
Bf4 = lambda U, Vv: (Qf4([u + v for u, v in zip(U, Vv)]) - Qf4(U) - Qf4(Vv)) / 2
V = [1 / (2 * Sf), -Pf / (2 * Sf), -Qf / (2 * Sf), (Pf**2 + Qf**2 + Sf**2) / (2 * Sf)]
Nn = [1, -p, -qq, p**2 + qq**2]
dV = [sp.diff(v, r) for v in V]
Y = [Phi * v for v in V]
dY = [sp.diff(v, r) for v in Y]
ok6 = (sp.simplify(Qf4(V) - 1) == 0
       and sp.simplify(-Qf4(dV) - (sp.diff(Pf, r)**2 + sp.diff(Qf, r)**2 + sp.diff(Sf, r)**2) / Sf**2) == 0
       and sp.simplify(Bf4(V, Nn) / 2 - Eq) == 0 and sp.simplify(Qf4(Nn)) == 0
       and sp.simplify(Qf4(dY) - (sp.diff(Phi, r)**2 - Phi**2 * (W.T * W)[0])) == 0)
check("Q(V) = 1, -Q(V') = |W|^2, E = B(V,N)/2 with N null, Q((Phi V)') = Phi'^2 - Phi^2 |W|^2", ok6)

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