"""test_627 -- zeta_TT(0) = +10 as a theorem: the regularised counts of the trace-free sector of the unit S^3.

The proof (derivation §89.3 and MA Remark 6.1): on a closed odd-dimensional manifold a self-adjoint elliptic
differential operator has no constant term in its heat expansion (Gilkey), so the regularised count of all its modes
is zero and zeta(0) = -dim ker. L_2 = -Delta + 2 on trace-free symmetric tensors has no kernel; the split
Sym^2_0 = TT + im K (K the conformal Killing operator) is preserved by every SO(4)-equivariant operator, since the two
lie in different SO(4) types; on the Einstein S^3 the Lichnerowicz Laplacian intertwines K, Delta_L K = K Delta_H, and
Delta_L = nabla*nabla + 6 on trace-free tensors, so L_2 K = K (Delta_H - 4): L_2 on im K is a differential operator on
vector fields with the ten conformal Killing modes removed, zeta = -10, and zeta_TT(0) = +10.

This script checks the bookkeeping with the explicit spectra (regularised sums by Hurwitz zeta, exact rationals):
  P1  TT: levels n >= 3, multiplicity 2(n^2 - 4), eigenvalue n^2 - 1: zeta(0) = +10.
  P2  im K, exact part: K(grad f), f at level n >= 3 (the constants and the four gradient conformal Killing fields
      removed), multiplicity n^2, eigenvalue n^2 - 1 - 4 = n^2 - 5: zeta(0) = -5.
  P3  im K, coexact part: coexact fields at level m = k + 1 >= 3 (the six Killing fields, k = 1, removed), multiplicity
      2(m^2 - 1), eigenvalue m^2 - 4: zeta(0) = -5; so im K counts -10 and Sym^2_0 counts 0, as the lemma requires.
  P4  scalar shape modes: levels n >= 3, multiplicity n^2: zeta(0) = -5 = -dim ker of Obata's trace-free Hessian
      (the constants and first harmonics). Controls of the lemma: all scalars (n >= 1, n^2) and all 1-forms (exact
      n >= 2, n^2; coexact m >= 2, 2(m^2 - 1)) count 0.
  P5  The counts do not depend on the eigenvalues within the class c n^p (1 + a_1 n^-2 + ...): P1 continued in s with
      lambda = n^2 - 1, n(n^2 - 1) and n^2 (n^2 - 1)^2 / 2 (MA Remark 6.1's multipliers), by subtracting the large-n
      expansion through n^-3 and summing the remainder. Control: lambda = n(n + 1), with an odd power, gives another value.
For lambda_n = n^2 (1 + even powers), zeta(0) of sum_n m(n) lambda_n^{-s} is the regularised sum of the polynomial m(n),
since no binomial term reaches a pole; this is MA Remark 6.1's continuation, applied to each sector.
History: Gemini (12 November 2026) proposed the result with an isospectrality step (K K* against K* K); the proof here
uses the Lichnerowicz intertwining instead, which also shows that L_2 preserves the split.
"""
import sys
import mpmath as mp
import sympy as sp
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)

def reg_sum(poly_coeffs, n0):
    """regularised sum over n >= n0 of sum_k c_k n^k = sum_k c_k zeta_H(-k, n0)"""
    return mp.fsum(c*mp.zeta(-k, n0) for k, c in poly_coeffs.items())

def zeta_cont(mcoef, C, p, rpoly, n0, s, J=8):
    """zeta(s) = sum_{n>=n0} m(n) lam(n)^{-s}, lam = C n^p (1 + r(1/n)), continued to small s: the large-n expansion of
    m(n) lam^{-s} through n^{-3} is subtracted (its sum is a Hurwitz zeta) and the convergent remainder summed."""
    uu, ss = sp.symbols('u s')
    G = sp.series((1 + rpoly(uu))**(-ss), uu, 0, J).removeO()
    terms = []
    for d, c in mcoef.items():
        for j in range(J):
            g = G.coeff(uu, j)
            if g != 0 and d - j >= -3:
                terms.append((d - j, c*g))
    sv = mp.mpf(s)
    tf = [(e, mp.mpf(sp.N(cg.subs(ss, sp.Float(str(sv), 40)), 40))) for e, cg in terms]
    lam = lambda n: mp.mpf(C)*mp.mpf(n)**p*(1 + rnum(rpoly, n))
    m = lambda n: mp.fsum(c*mp.mpf(n)**d for d, c in mcoef.items())
    rem = mp.nsum(lambda n: m(n)*lam(n)**(-sv) - mp.mpf(C)**(-sv)*mp.fsum(cg*mp.mpf(n)**(e - p*sv) for e, cg in tf), [n0, mp.inf])
    hur = mp.mpf(C)**(-sv)*mp.fsum(cg*mp.zeta(p*sv - e, n0) for e, cg in tf)
    return rem + hur

def rnum(rpoly, n):
    return mp.mpf(rpoly(1/mp.mpf(n)))

class Poly:
    def __init__(self, coeffs): self.coeffs = coeffs
    def __call__(self, n): return mp.fsum(c*mp.mpf(n)**k for k, c in self.coeffs.items())

TT = Poly({2: 2, 0: -8}); EX = Poly({2: 1}); CO = Poly({2: 2, 0: -2}); SC = Poly({2: 1})
z_tt = reg_sum(TT.coeffs, 3); z_ex = reg_sum(EX.coeffs, 3); z_co = reg_sum(CO.coeffs, 3); z_sc = reg_sum(SC.coeffs, 3)
check("P1 TT: sum_{n>=3} 2(n^2 - 4) regularised = +10", abs(z_tt - 10) < mp.mpf('1e-25'), mp.nstr(z_tt, 10))
check("P2 im K exact part: sum_{n>=3} n^2 regularised = -5", abs(z_ex + 5) < mp.mpf('1e-25'), mp.nstr(z_ex, 10))
check("P3 im K coexact part: sum_{m>=3} 2(m^2 - 1) regularised = -5; im K = -10; Sym^2_0 = TT + im K = 0",
      abs(z_co + 5) < mp.mpf('1e-25') and abs(z_tt + z_ex + z_co) < mp.mpf('1e-25'), mp.nstr(z_co, 10))
z_scal_all = reg_sum({2: 1}, 1); z_1f_all = reg_sum({2: 1}, 2) + reg_sum({2: 2, 0: -2}, 2)
check("P4 scalar shape modes n >= 3: -5 = -dim ker(Obata's trace-free Hessian); all scalars and all 1-forms count 0",
      abs(z_sc + 5) < mp.mpf('1e-25') and abs(z_scal_all) < mp.mpf('1e-25') and abs(z_1f_all) < mp.mpf('1e-25'),
      f"scalars {mp.nstr(z_scal_all, 5)}, 1-forms {mp.nstr(z_1f_all, 5)}")
ok5 = True; vals = []
for C, p, r in [(1, 2, lambda u: -u**2), (1, 3, lambda u: -u**2), (sp.Rational(1, 2), 6, lambda u: (1 - u**2)**2 - 1)]:
    z = zeta_cont(TT.coeffs, C, p, r, 3, '1e-20')
    vals.append(z); ok5 &= abs(z - 10) < mp.mpf('1e-12')
zc = zeta_cont(TT.coeffs, 1, 2, lambda u: u, 3, '1e-20')   # control: lam = n^2 + n has an odd power
check("P5 the TT count is +10 for lambda = n^2 - 1, n(n^2 - 1), n^2 (n^2 - 1)^2/2 (continuation in s); "
      "control: an odd power, lambda = n(n + 1), shifts it", ok5 and abs(zc - 10) > mp.mpf('0.1'),
      ", ".join(mp.nstr(v, 12) for v in vals) + f"; control {mp.nstr(zc, 12)}")
print("\nALL PASS" if not FAIL else f"\nFAILURES: {FAIL}")
sys.exit(1 if FAIL else 0)