"""test_625 -- at the round S^3 the linearised realizability operator is the linearised Cotton operator. F(g) = Cod_g(K~), K~ = sqrt(det M) M^{-1} lowered with g, M = -g^{-1} G (MA (1.3)); C(g) = Cod_g(P), P = Ric - R g/4 (the Schouten tensor in dimension three, so C is the Cotton tensor). At the unit round metric gamma, K~ = gamma and P = gamma/2, both Codazzi, so F(gamma) = C(gamma) = 0. P1 [exact] DF_gamma(H) = DC_gamma(H) componentwise, for general polynomial H (integer coefficients, degree <= 3 in each component) at two rational points of S^3 in stereographic coordinates, gamma = 4|dy|^2/(1 + |y|^2)^2. Both sides are computed from scratch (Christoffels, Ricci, M, the root, Codazzi; and Ricci, Schouten, Codazzi) in exact dual-number arithmetic, a + eps b, to first order in eps. P2 [exact] Controls: on H = 2u gamma (u a polynomial) and on H = L_X gamma (X a polynomial field) both sides vanish, as they must (C vanishes on conformally flat metrics; F's kernel). So MA Remark 6.1's coincidence of multipliers (DF_gamma on TT has singular values n(n^2 - 1)/sqrt 2, Branson's conformal P_3) is an identity of operators: on TT, DF_gamma is the linearised Cotton operator, which is conformally covariant. Proposed by Gemini (12 November 2026); checked here. A first version of this script passed values through sympy's nsimplify, which replaced exact rationals by nearby radicals and made the Lie control fail spuriously; all arithmetic is now exact rational, with no numerical identification. """ import sys import itertools import random 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) class Du: """a + eps b, exact.""" __slots__ = ("a", "b") def __init__(self, a, b=0): self.a = sp.sympify(a); self.b = sp.sympify(b) def __add__(s, o): o = o if isinstance(o, Du) else Du(o); return Du(s.a + o.a, s.b + o.b) __radd__ = __add__ def __sub__(s, o): o = o if isinstance(o, Du) else Du(o); return Du(s.a - o.a, s.b - o.b) def __rsub__(s, o): return Du(o) - s def __neg__(s): return Du(-s.a, -s.b) def __mul__(s, o): o = o if isinstance(o, Du) else Du(o); return Du(s.a*o.a, s.a*o.b + s.b*o.a) __rmul__ = __mul__ def inv(s): return Du(1/s.a, -s.b/s.a**2) def __truediv__(s, o): o = o if isinstance(o, Du) else Du(o); return s*o.inv() def sqrt(s): r = sp.sqrt(s.a); return Du(r, s.b/(2*r)) y = sp.symbols('y0:3', real=True) eps = sp.Symbol('eps') n = 3 def gamma_stereo(): c = 4/(1 + sum(v**2 for v in y))**2 return sp.Matrix(3, 3, lambda i, j: c if i == j else 0) def run(Hexpr, pt): """Return (F_eps, C_eps, formula_eps): dicts of the eps-coefficients of the Codazzi components.""" g = gamma_stereo() + eps*Hexpr ev = lambda e: Du(sp.cancel(e.subs(eps, 0).subs(pt)), sp.cancel(sp.diff(e, eps).subs(eps, 0).subs(pt))) cache = {} def Dg(i, j, *ks): key = (min(i, j), max(i, j)) + tuple(sorted(ks)) if key not in cache: cache[key] = ev(sp.diff(g[i, j], *[y[k] for k in ks]) if ks else g[i, j]) return cache[key] G0 = [[Dg(i, j) for j in range(n)] for i in range(n)] # inverse of a 3x3 Du matrix by adjugate def det3(A): return (A[0][0]*(A[1][1]*A[2][2] - A[1][2]*A[2][1]) - A[0][1]*(A[1][0]*A[2][2] - A[1][2]*A[2][0]) + A[0][2]*(A[1][0]*A[2][1] - A[1][1]*A[2][0])) def inv3(A): d = det3(A) cof = [[None]*3 for _ in range(3)] for i in range(3): for j in range(3): m = [[A[r][c] for c in range(3) if c != j] for r in range(3) if r != i] cof[i][j] = (m[0][0]*m[1][1] - m[0][1]*m[1][0])*(1 if (i + j) % 2 == 0 else -1) return [[cof[j][i]/d for j in range(3)] for i in range(3)], d gi, _ = inv3(G0) dgi = [[[sum((-(gi[i][p]*Dg(p, q, m)*gi[q][l]) for p in range(n) for q in range(n)), Du(0)) for m in range(n)] for l in range(n)] for i in range(n)] def ddgi(i, l, m, k): return sum((-(dgi[i][p][k]*Dg(p, q, m)*gi[q][l] + gi[i][p]*Dg(p, q, m, k)*gi[q][l] + gi[i][p]*Dg(p, q, m)*dgi[q][l][k]) for p in range(n) for q in range(n)), Du(0)) def C3(l, j, k, *ms): return (Dg(l, j, k, *ms) + Dg(l, k, j, *ms) - Dg(j, k, l, *ms))*Q(1, 2) Gm = [[[sum((gi[i][l]*C3(l, j, k) for l in range(n)), Du(0)) for k in range(n)] for j in range(n)] for i in range(n)] dG = [[[[sum((dgi[i][l][m]*C3(l, j, k) + gi[i][l]*C3(l, j, k, m) for l in range(n)), Du(0)) for m in range(n)] for k in range(n)] for j in range(n)] for i in range(n)] def ddG(i, j, k, m, q): return sum((ddgi(i, l, m, q)*C3(l, j, k) + dgi[i][l][m]*C3(l, j, k, q) + dgi[i][l][q]*C3(l, j, k, m) + gi[i][l]*C3(l, j, k, m, q) for l in range(n)), Du(0)) def Ric(j, k): return sum((dG[i][j][k][i] - dG[i][j][i][k] + sum((Gm[i][i][l]*Gm[l][j][k] - Gm[i][k][l]*Gm[l][j][i] for l in range(n)), Du(0)) for i in range(n)), Du(0)) def dRic(j, k, q): return sum((ddG(i, j, k, i, q) - ddG(i, j, i, k, q) + sum((dG[i][i][l][q]*Gm[l][j][k] + Gm[i][i][l]*dG[l][j][k][q] - dG[i][k][l][q]*Gm[l][j][i] - Gm[i][k][l]*dG[l][j][i][q] for l in range(n)), Du(0)) for i in range(n)), Du(0)) Rc = [[Ric(j, k) for k in range(n)] for j in range(n)] dRc = [[[dRic(j, k, q) for q in range(n)] for k in range(n)] for j in range(n)] Rs = sum((gi[j][k]*Rc[j][k] for j in range(n) for k in range(n)), Du(0)) dRs = [sum((dgi[j][k][q]*Rc[j][k] + gi[j][k]*dRc[j][k][q] for j in range(n) for k in range(n)), Du(0)) for q in range(n)] def cod(T, dT): # T_jk lowered, dT[j][k][a] = d_a T_jk out = {} for a in range(n): for j in range(n): if not a < j: continue for k in range(n): t = dT[j][k][a] - dT[a][k][j] for m in range(n): t = t - Gm[m][a][j]*T[m][k] - Gm[m][a][k]*T[j][m] + Gm[m][j][a]*T[m][k] + Gm[m][j][k]*T[a][m] out[(a, j, k)] = t return out # Cotton P = [[Rc[i][j] - Rs*G0[i][j]*Q(1, 4) for j in range(n)] for i in range(n)] dP = [[[dRc[i][j][q] - (dRs[q]*G0[i][j] + Rs*Dg(i, j, q))*Q(1, 4) for q in range(n)] for j in range(n)] for i in range(n)] Cc = cod(P, dP) # realizability: Glow = Ric - R g/2, M^a_b = -g^{ac} G_cb Gl = [[Rc[i][j] - Rs*G0[i][j]*Q(1, 2) for j in range(n)] for i in range(n)] dGl = [[[dRc[i][j][q] - (dRs[q]*G0[i][j] + Rs*Dg(i, j, q))*Q(1, 2) for q in range(n)] for j in range(n)] for i in range(n)] Mm = [[-sum((gi[a][c]*Gl[c][b] for c in range(n)), Du(0)) for b in range(n)] for a in range(n)] dMm = [[[-sum((dgi[a][c][q]*Gl[c][b] + gi[a][c]*dGl[c][b][q] for c in range(n)), Du(0)) for q in range(n)] for b in range(n)] for a in range(n)] Mi, dM = inv3(Mm) sq = dM.sqrt() trMidM = [sum((Mi[i][j]*dMm[j][i][q] for i in range(n) for j in range(n)), Du(0)) for q in range(n)] dsq = [sq*trMidM[q]*Q(1, 2) for q in range(n)] dMi = [[[-sum((Mi[i][p]*dMm[p][r][q]*Mi[r][j] for p in range(n) for r in range(n)), Du(0)) for q in range(n)] for j in range(n)] for i in range(n)] Kmix = [[sq*Mi[a][b] for b in range(n)] for a in range(n)] dKmix = [[[dsq[q]*Mi[a][b] + sq*dMi[a][b][q] for q in range(n)] for b in range(n)] for a in range(n)] Kl = [[sum((G0[j][c]*Kmix[c][k] for c in range(n)), Du(0)) for k in range(n)] for j in range(n)] dKl = [[[sum((Dg(j, c, q)*Kmix[c][k] + G0[j][c]*dKmix[c][k][q] for c in range(n)), Du(0)) for q in range(n)] for k in range(n)] for j in range(n)] Fc = cod(Kl, dKl) # Gemini's formula: Cod_gamma(dRic - dR gamma/4 - 2H): the eps-part of the tensor T = Ric - R g/4 - 2 eps H ... built # as T = P - (3/2)(g - gamma) (so dT = dRic - dR gamma/4 - (1/4) R0 H - (3/2) H = ... - 2H since R0 = 6), at eps order 1 T = [[P[i][j] - (G0[i][j] - Du(G0[i][j].a))*Q(3, 2) for j in range(n)] for i in range(n)] dT = [[[dP[i][j][q] - (Dg(i, j, q) - Du(Dg(i, j, q).a))*Q(3, 2) for q in range(n)] for j in range(n)] for i in range(n)] # Codazzi of the eps-part with the round connection: Cod_gamma(T1) = eps-part of Cod_g(T) minus Cod-variation acting on T0 = gamma/2, # which is -(1/2) Cod_gamma(H); so Cod_gamma(T1) = [Cod_g(T)]_eps + (1/2) Cod_gamma(H) ... computed via the identity below Tc = cod(T, dT) return ({k: sp.simplify(v.b) for k, v in Fc.items()}, {k: sp.simplify(v.b) for k, v in Cc.items()}, {k: sp.simplify(v.a) for k, v in Fc.items()}, {k: sp.simplify(v.a) for k, v in Cc.items()}, Tc) random.seed(625) def rand_poly(deg=3): terms = [sp.Integer(1)] + list(y) + [y[i]*y[j] for i in range(3) for j in range(i, 3)] + [y[i]*y[j]*y[k] for i in range(3) for j in range(i, 3) for k in range(j, 3)] return sum(random.randint(-3, 3)*t for t in terms) 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)}] ok1 = True; ok0 = True; detail = "" for pt in pts: 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() Fe, Ce, F0, C0, _ = run(H, pt) ok0 &= all(v == 0 for v in F0.values()) and all(v == 0 for v in C0.values()) ok1 &= all(sp.simplify(Fe[k] - Ce[k]) == 0 for k in Fe) and any(Fe[k] != 0 for k in Fe) detail = f"e.g. (0,1,2): DF = {Fe[(0, 1, 2)]}, DC = {Ce[(0, 1, 2)]}" check("P1 F(gamma) = C(gamma) = 0, and DF_gamma(H) = DC_gamma(H) for general H at two rational points", ok0 and ok1, detail) # P2 controls pt = pts[0] gam = gamma_stereo() u = rand_poly(2) H_conf = 2*u*gam X = [rand_poly(2) for _ in range(3)] Xl = [sum(gam[i, j]*X[j] for j in range(3)) for i in range(3)] Gs = [[[sum(gam.inv()[l, m]*(sp.diff(gam[m, i], y[j]) + sp.diff(gam[m, j], y[i]) - sp.diff(gam[i, j], y[m])) for m in range(3))/2 for j in range(3)] for i in range(3)] for l in range(3)] H_lie = sp.Matrix(3, 3, lambda i, j: sp.diff(Xl[j], y[i]) + sp.diff(Xl[i], y[j]) - 2*sum(Gs[l][i][j]*Xl[l] for l in range(3))) ok2 = True for Hc in (H_conf, H_lie): Fe, Ce, _, _, _ = run(Hc, pt) ok2 &= all(v == 0 for v in Fe.values()) and all(v == 0 for v in Ce.values()) check("P2 controls: both sides vanish on 2u gamma and on L_X gamma", ok2) print("\nALL PASS" if not FAIL else f"\nFAILURES: {FAIL}") sys.exit(1 if FAIL else 0)