"""test_628 -- MB open problem 7 in the O(2)-invariant sector: (df (x) df)^TT = 0 iff the graph of F is developable.

Let f on the unit S^3 be invariant under the O(2) acting on (x2, x3), so f = F(x0, x1)|_{S^3} with F smooth on a
neighbourhood of the closed unit disk D. Then (df (x) df)^TT = 0 if and only if det D^2 F = F_00 F_11 - F_01^2 = 0 on D.
Proof. If df (x) df = 2 w gamma + L_X gamma, average (w, X) over O(2) (the map is equivariant and f, gamma invariant), so
X = Y^T with Y = A d_0 + B d_1, A and B functions of (x0, x1); the reflection x3 -> -x3 removes a d_phi part. Since
iota^* L_Y delta = L_{Y^T} gamma + 2 (Y . x) gamma and L_Y delta has only dx0, dx1 components, the phi-phi component forces
w = x0 A + x1 B, and the equation becomes F_0^2 = 2 A_0, F_0 F_1 = A_1 + B_0, F_1^2 = 2 B_1 on the open disk (dx0, dx1 are
chart forms there). Two-dimensional Saint-Venant: d_1^2 (F_0^2) - 2 d_0 d_1 (F_0 F_1) + d_0^2 (F_1^2) = -2 det D^2 F, so the
determinant vanishes. Conversely, if det D^2 F = 0 on the open disk, the same identity and Saint-Venant on the
(simply connected) disk give A, B, hence the decomposition over the disk, so DF_gamma(df (x) df) = 0 there; DF_gamma is a
differential operator, so it vanishes on all of S^3 by continuity, and MA Theorem E(1) gives no TT part.
So in this sector the solutions are exactly the developable graphs z = F(x0, x1) over D: cylinders, cones and tangent
developables (and their smooth gluings along rulings).

  P1 [symbolic] d_1^2 (F_0^2) - 2 d_0 d_1 (F_0 F_1) + d_0^2 (F_1^2) = -2 det D^2 F for a general F(x0, x1).
  P2 [symbolic] In the chart (x0, x1, phi) of S^3, for Y = A(x0, x1) d_0 + B(x0, x1) d_1: L_{Y^T} gamma =
     iota^* L_Y delta - 2 (x0 A + x1 B) gamma; iota^* L_Y delta = 2 A_0 dx0^2 + 2 (A_1 + B_0) dx0 dx1 + 2 B_1 dx1^2 has no phi
     components; so the phi-phi component of df (x) df = 2 w gamma + L_{Y^T} gamma forces w = x0 A + x1 B.
  P3 [numeric] The converse direction at work: O(2)-invariant F with det D^2 F != 0 give Q < 0 -- x0 x1 (det -1, the
     control -1/24), x0^2 + x1^2 (det 4), x0^2 - x1^2 + x0 x1 (det -5).
  P4 [numeric] A tangent developable: the plane curve gamma(t) = (-3 + cos t, sin t) lifted to the space curve
     (gamma(t), t^2); F(gamma(t) + s gamma'(t)) = t^2 + 2 s t on the branch s > 0, smooth on the disk of radius 3/2 (the
     convex set {x0^2 + x1^2 < 9/4} x R^2 contains S^3). det D^2 F = 0 at 25 random points to 1e-25 (30 digits);
     the rulings are not concurrent and not parallel (a tangent developable, neither cone nor cylinder); Q = 0 to
     1e-12; not zonal about any axis in the (x0, x1) plane.
History: prompted by Gemini's note (12 November 2026, section 3.2), which sketched the equivalence. P1-P3 passed on the
first run; P4 stopped on a numpy cross-product call for 2-vectors, replaced by the explicit 2x2 determinant, claim unchanged.
"""
import sys, os
import numpy as np
import sympy as sp
import mpmath as mp
mp.mp.dps = 30
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)

x = sp.symbols('x0:4', real=True)
# P1
a0, a1 = sp.symbols('a0 a1', real=True)
Fg = sp.Function('F')(a0, a1)
F0, F1 = sp.diff(Fg, a0), sp.diff(Fg, a1)
lhs = sp.diff(F0**2, a1, 2) - 2*sp.diff(F0*F1, a0, a1) + sp.diff(F1**2, a0, 2)
det = sp.diff(Fg, a0, 2)*sp.diff(Fg, a1, 2) - sp.diff(Fg, a0, a1)**2
check("P1 the 2-d Saint-Venant operator on dF (x) dF is -2 det D^2 F", sp.simplify(lhs + 2*det) == 0)

# P2
u0, u1, ph = sp.symbols('u0 u1 phi', real=True)
rho = sp.sqrt(1 - u0**2 - u1**2)
emb = sp.Matrix([u0, u1, rho*sp.cos(ph), rho*sp.sin(ph)])
U = [u0, u1, ph]
E = emb.jacobian(U)                      # columns e_a
gam = sp.simplify(E.T*E)
A = sp.Function('A')(u0, u1); B = sp.Function('B')(u0, u1)
Yamb = sp.Matrix([A, B, 0, 0])
cov = sp.simplify(E.T*Yamb)              # Y . e_a
Xc = sp.simplify(gam.inv()*cov)          # components of Y^T in the chart
def lie(Xv, g):
    return sp.Matrix(3, 3, lambda a, b: sum(Xv[c]*sp.diff(g[a, b], U[c]) + g[c, b]*sp.diff(Xv[c], U[a])
                                            + g[a, c]*sp.diff(Xv[c], U[b]) for c in range(3)))
LX = lie(Xc, gam)
# iota^* L_Y delta: (dY_i/dx_j + dY_j/dx_i) e_a^i e_b^j; Y depends on ambient x0, x1 = u0, u1
J = sp.Matrix(4, 4, lambda i, j: 0)
Ycomp = [A, B, 0, 0]
amb_d = {0: u0, 1: u1}
for i in range(4):
    for j in range(4):
        dij = sp.diff(Ycomp[i], amb_d[j]) if j in amb_d else 0
        dji = sp.diff(Ycomp[j], amb_d[i]) if i in amb_d else 0
        J[i, j] = dij + dji
pull = sp.simplify(E.T*J*E)
ydotx = u0*A + u1*B
ok_id = sp.simplify(LX - (pull - 2*ydotx*gam)) == sp.zeros(3, 3)
ok_phi = all(sp.simplify(pull[a, 2]) == 0 for a in range(3))
target = sp.Matrix([[2*sp.diff(A, u0), sp.diff(A, u1) + sp.diff(B, u0), 0],
                    [sp.diff(A, u1) + sp.diff(B, u0), 2*sp.diff(B, u1), 0], [0, 0, 0]])
ok_form = sp.simplify(pull - target) == sp.zeros(3, 3)
check("P2 L_{Y^T} gamma = iota^* L_Y delta - 2 (x0 A + x1 B) gamma; iota^* L_Y delta has no phi components", ok_id and ok_phi and ok_form)

# Q by quadrature (test_624's routine, the same integrand as test_353)
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
src = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "test_624_nonzonal_tt_free.py")).read()
ns = {"np": np, "sp": sp, "x": x}
exec(src[src.index("def Q_of"):src.index("ctrl1 = Q_of")], ns)
Q_of = ns["Q_of"]

# P3
ctrl = {"x0 x1": x[0]*x[1], "x0^2 + x1^2": x[0]**2 + x[1]**2, "x0^2 - x1^2 + x0 x1": x[0]**2 - x[1]**2 + x[0]*x[1]}
rq = {k: Q_of(F_, n=120)[3] for k, F_ in ctrl.items()}
check("P3 O(2)-invariant F with det D^2 F != 0 have Q < 0", abs(rq["x0 x1"] + 1/24) < 1e-12 and all(r < -1e-6 for r in rq.values()),
      "; ".join(f"{k}: {r:.4e}" for k, r in rq.items()))

# P4: the tangent developable
X0, X1 = x[0], x[1]
d = sp.sqrt((X0 + 3)**2 + X1**2)
phi_ = sp.atan(X1/(X0 + 3))
t_ = phi_ - sp.acos(1/d)
s_ = sp.sqrt(d**2 - 1)
Ftd = t_**2 + 2*s_*t_
# det at random points, 30 digits
H2 = [[sp.diff(Ftd, X0, X0), sp.diff(Ftd, X0, X1)], [sp.diff(Ftd, X0, X1), sp.diff(Ftd, X1, X1)]]
detf = sp.lambdify((X0, X1), H2[0][0]*H2[1][1] - H2[0][1]**2, 'mpmath')
scale = sp.lambdify((X0, X1), abs(H2[0][0]) + abs(H2[1][1]) + abs(H2[0][1]), 'mpmath')
rng = np.random.default_rng(628)
worst = 0
for _ in range(25):
    r_, th = 1.2*np.sqrt(rng.random()), 2*np.pi*rng.random()
    p0, p1 = mp.mpf(r_*np.cos(th)), mp.mpf(r_*np.sin(th))
    worst = max(worst, abs(detf(p0, p1))/scale(p0, p1)**2)
# rulings: ruling through gamma(t) has direction gamma'(t); concurrency / parallelism fails for three rulings
ts = [0.3, 1.0, 1.7]
lines = [(np.array([-3 + np.cos(t), np.sin(t)]), np.array([-np.sin(t), np.cos(t)])) for t in ts]
def inter(L1, L2):
    (p, u), (q, w) = L1, L2
    M = np.column_stack([u, -w]); a_ = np.linalg.solve(M, q - p); return p + a_[0]*u
P12, P13 = inter(lines[0], lines[1]), inter(lines[0], lines[2])
not_cone = np.linalg.norm(P12 - P13) > 1e-3
u_, w_ = lines[0][1], lines[1][1]; not_cyl = abs(u_[0]*w_[1] - u_[1]*w_[0]) > 1e-3
qtd = Q_of(Ftd, n=120)
ftd = sp.lambdify((X0, X1), Ftd, 'numpy')
gapmin = np.inf
for beta in np.linspace(0, np.pi, 180, endpoint=False):
    a = np.array([np.cos(beta), np.sin(beta), 0, 0]); e = np.array([-np.sin(beta), np.cos(beta), 0, 0]); e2 = np.array([0, 0, 1.0, 0])
    gap = max(abs(ftd(*(uu*a + np.sqrt(1 - uu*uu)*e)[:2]) - ftd(*(uu*a + np.sqrt(1 - uu*uu)*e2)[:2])) for uu in np.linspace(-0.95, 0.95, 39))
    gapmin = min(gapmin, gap)
check("P4 a tangent developable: det D^2 F = 0, neither cone nor cylinder, Q = 0, not zonal",
      worst < mp.mpf('1e-25') and not_cone and not_cyl and np.max(np.abs(qtd)) < 1e-12 and gapmin > 1e-3,
      f"max |det|/scale^2 = {mp.nstr(worst, 3)}; Q/2pi^2 = {qtd[3]:.1e}; gap {gapmin:.3f}")
print("\nALL PASS" if not FAIL else f"\nFAILURES: {FAIL}")
sys.exit(1 if FAIL else 0)