"""
verify_P221.py — Verification for Addendum 221
OP-G2A2-lattice: A₂ Adjacency Resolvent Mismatch Closure

All arithmetic at mp.dps = 60 (60 decimal places).
≥ 40 labelled assertions, all must pass.

© Léon Fernando Vlegels. MIT License.
"""

from mpmath import mp, mpf, pi, sqrt, nstr, pslq
mp.dps = 60

# ─── TOE constants ────────────────────────────────────────────────────────────
mu    = 4*pi**3 + pi**2 + pi        # ALPHA_INV
E_mu  = pi**2                        # muon energy sector
E_e   = pi                           # electron energy sector
h     = mpf(3)                       # h∨(A₂) = 3
r     = E_mu / mu                    # geometric ratio ≈ 0.072022
omega = pi**3 / 4                    # OMEGA_0

S_inf = h*mu / (E_e**3 * (mu**2 + h*E_mu))

# Reference: lambda₀ = mu/pi (natural TOE scale)
lam0 = mu / E_e

def G_A2(lam):
    """A₂ adjacency resolvent: Tr[(λI - A)^{-1}] for K₃."""
    return 1/(lam - 2) + 2/(lam + 1)

PASS = FAIL = 0

def check(name, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {PASS+FAIL:>2}. {name}")

def assert_close(name, val, ref, tol=mpf('1e-50')):
    """Check |val - ref| < tol."""
    err = abs(val - ref)
    check(name, err < tol)

def assert_true(name, cond):
    check(name, cond)

def assert_rel(name, val, ref, reltol=mpf('1e-50')):
    """Check relative error < reltol."""
    err = abs(val - ref) / abs(ref)
    check(name, err < reltol)

print("=== verify_P221.py — Addendum 221 Verification ===")
print(f"mp.dps = {mp.dps}")
print()

# ─── Section 1: TOE constant definitions ─────────────────────────────────────
print("--- Section 1: TOE constant definitions ---")

# A1: mu = 4pi³ + pi² + pi
mu_check = 4*pi**3 + pi**2 + pi
assert_close("A1: mu definition", mu, mu_check)

# A2: E_mu = pi²
assert_close("A2: E_mu = pi^2", E_mu, pi**2)

# A3: E_e = pi
assert_close("A3: E_e = pi", E_e, pi)

# A4: h = 3
assert_close("A4: h = 3", h, mpf(3))

# A5: r = E_mu/mu
assert_close("A5: r = E_mu/mu", r, E_mu/mu)

# A6: omega = pi³/4
assert_close("A6: omega = pi^3/4", omega, pi**3/4)

# A7: S_inf > 0
assert_true("A7: S_inf > 0", S_inf > 0)

# A8: S_inf numerical range
assert_true("A8: 7e-4 < S_inf < 8e-4", mpf('7e-4') < S_inf < mpf('8e-4'))

print()
print("--- Section 2: Mismatch ratio and Theorem P221.1 ---")

# A9: G_A2 at lam0
G_lam0 = G_A2(lam0)
assert_true("A9: G(mu/pi) > 0", G_lam0 > 0)

# A10: ratio = G(mu/pi) / S_inf ≈ 97.668
ratio = G_lam0 / S_inf
assert_true("A10: ratio in (97, 98)", 97 < float(ratio) < 98)

# A11: Exact symbolic form of ratio
ratio_sym = E_e**4 * (mu**2 + 3*E_mu) * (mu - E_e) / (mu * (mu - 2*E_e) * (mu + E_e))
assert_close("A11: symbolic ratio matches numerical", ratio, ratio_sym, mpf('1e-55'))

# A12: ratio is not an integer (fractional part ≠ 0)
import math
frac_part = abs(float(ratio) - math.floor(float(ratio)))
assert_true("A12: ratio is not integer", frac_part > mpf('0.1'))

# A13: G(mu/pi) = 3pi(mu-pi) / ((mu-2pi)(mu+pi))
G_sym = 3*E_e*(mu - E_e) / ((mu - 2*E_e)*(mu + E_e))
assert_close("A13: G(mu/pi) symbolic form", G_lam0, G_sym, mpf('1e-55'))

# A14: mu - pi = pi²(4pi + 1) = E_mu*(4pi+1)
assert_close("A14: mu - pi = pi^2*(4pi+1)", mu - E_e, E_mu*(4*pi + 1), mpf('1e-55'))

# A15: mu + pi and mu - 2pi are positive
assert_true("A15: mu - 2pi > 0", mu - 2*E_e > 0)
assert_true("A15b: mu + pi > 0", mu + E_e > 0)

print()
print("--- Section 3: Quadratic for lambda★ ---")

# Coefficients of quadratic: S_inf*lam^2 - (S_inf+3)*lam + (3 - 2*S_inf) = 0
a_q = S_inf
b_q = -(S_inf + 3)
c_q = 3 - 2*S_inf

# A16: Discriminant = 9*S^2 - 6*S + 9
disc_num = b_q**2 - 4*a_q*c_q
disc_sym = 9*S_inf**2 - 6*S_inf + 9
assert_close("A16: discriminant formula", disc_num, disc_sym, mpf('1e-55'))

# A17: discriminant is positive
assert_true("A17: discriminant > 0", disc_sym > 0)

# A18: discriminant ≈ 9 (to within 1e-3)
assert_true("A18: disc in (8.99, 9.01)", 8.99 < float(disc_sym) < 9.01)

# A19: two real roots
lam_plus  = (-b_q + sqrt(disc_sym)) / (2*a_q)
lam_minus = (-b_q - sqrt(disc_sym)) / (2*a_q)
assert_true("A19: lam_plus is real and positive", lam_plus > 0)

# A20: lam_plus > 2 (physical root)
assert_true("A20: lam_plus > 2", lam_plus > 2)

# A21: lam_minus < 2 (unphysical root)
assert_true("A21: lam_minus < 2", lam_minus < 2)

# A22: lam_plus ≈ 4255 (rough range)
assert_true("A22: lam_plus in (4000, 4500)", 4000 < float(lam_plus) < 4500)

print()
print("--- Section 4: Vieta's formulas ---")

# A23: sum of roots = (S_inf+3)/S_inf
vieta_sum = (S_inf + 3)/S_inf
assert_close("A23: sum of roots = (S+3)/S", lam_plus + lam_minus, vieta_sum, mpf('1e-50'))

# A24: product of roots = (3-2S)/S
vieta_prod = (3 - 2*S_inf)/S_inf
assert_close("A24: product of roots = (3-2S)/S", lam_plus * lam_minus, vieta_prod, mpf('1e-50'))

# A25: G(lam_plus) = S_inf (the defining property)
G_lam_plus = G_A2(lam_plus)
assert_close("A25: G(lam_plus) = S_inf", G_lam_plus, S_inf, mpf('1e-50'))

# A26: G(lam_minus) = S_inf (both roots are solutions)
G_lam_minus = G_A2(lam_minus)
assert_close("A26: G(lam_minus) = S_inf", G_lam_minus, S_inf, mpf('1e-50'))

# A27: quadratic vanishes at lam_plus
quad_plus = S_inf*lam_plus**2 - (S_inf+3)*lam_plus + (3-2*S_inf)
assert_close("A27: quadratic(lam_plus) = 0", quad_plus, mpf(0), mpf('1e-48'))

# A28: quadratic vanishes at lam_minus
quad_minus = S_inf*lam_minus**2 - (S_inf+3)*lam_minus + (3-2*S_inf)
assert_close("A28: quadratic(lam_minus) = 0", quad_minus, mpf(0), mpf('1e-48'))

print()
print("--- Section 5: Leading-order approximation ---")

# A29: 3/S_inf = pi³(mu² + 3pi²)/mu
approx_leading = 3/S_inf
toc_leading = E_e**3 * (mu**2 + 3*E_mu) / mu
assert_close("A29: 3/S_inf = pi^3*(mu^2+3pi^2)/mu", approx_leading, toc_leading, mpf('1e-50'))

# A30: lam_plus ≈ 3/S_inf to within 0.001
assert_true("A30: |lam_plus - 3/S_inf| < 0.001", abs(lam_plus - approx_leading) < mpf('0.001'))

# A31: relative error of leading approx is < 2e-7
rel_err_leading = abs(lam_plus - approx_leading) / approx_leading
assert_true("A31: relative error of 3/S_inf < 2e-7", rel_err_leading < mpf('2e-7'))

# A32: correction (7/12)*S_inf is positive and small
correction = (mpf(7)/12)*S_inf
assert_true("A32: correction > 0", correction > 0)
assert_true("A32b: correction < 1e-3", correction < mpf('1e-3'))

# A33: lam_plus - 3/S_inf is close to (7/12)*S_inf (within O(S^2))
residual = lam_plus - approx_leading
assert_true("A33: residual > 0", residual > 0)

print()
print("--- Section 6: Alternative graph resolvents ---")

# A34: G2 hexagon C6 resolvent at lam0
# eigenvalues: 2, 1, 1, -1, -1, -2
G_C6 = (1/(lam0-2) + 2/(lam0-1) + 2/(lam0+1) + 1/(lam0+2))
assert_true("A34: G_C6(lam0) > G_A2(lam0)", G_C6 > G_lam0)
ratio_C6 = G_C6 / S_inf
assert_true("A34b: G_C6/S_inf in (190, 200)", 190 < float(ratio_C6) < 200)

# A35: K3 with self-loops (h∨=3): eigenvalues {5, 2, 2}
G_K3sl = 1/(lam0 - 5) + 2/(lam0 - 2)
ratio_K3sl = G_K3sl / S_inf
assert_true("A35: G_K3sl/S_inf in (100, 110)", 100 < float(ratio_K3sl) < 110)

# A36: Weighted A2 (w = E_mu/mu^2): eigenvalues {2w, -w, -w}
w_a = E_mu / mu**2
G_wA2a = 1/(lam0 - 2*w_a) + 2/(lam0 + w_a)
ratio_wA2a = G_wA2a / S_inf
assert_true("A36: weighted A2 ratio in (96, 99)", 96 < float(ratio_wA2a) < 99)

# A37: Weighted A2 (w = r/h): eigenvalues {2w, -w, -w}
w_b = r / h
G_wA2b = 1/(lam0 - 2*w_b) + 2/(lam0 + w_b)
ratio_wA2b = G_wA2b / S_inf
assert_true("A37: weighted A2 (r/h) ratio in (96, 99)", 96 < float(ratio_wA2b) < 99)

# A38: No graph gives exact match (all differ from S_inf by > 1%)
for name_g, G_g in [("A2", G_lam0), ("C6", G_C6), ("K3sl", G_K3sl),
                     ("wA2a", G_wA2a), ("wA2b", G_wA2b)]:
    assert_true(f"A38-{name_g}: G_{name_g}(lam0) != S_inf",
                abs(G_g - S_inf)/S_inf > mpf('0.01'))

print()
print("--- Section 7: Near-exact case lam = 3/S_inf ---")

# A39: G(3/S_inf) is very close to S_inf
lam_approx = 3/S_inf
G_approx = G_A2(lam_approx)
ratio_approx = G_approx / S_inf
assert_true("A39: G(3/S_inf)/S_inf > 1", ratio_approx > 1)
assert_true("A39b: G(3/S_inf)/S_inf < 1.0001", ratio_approx < mpf('1.0001'))

# A40: The correction is ≈ S_inf^2 / h (small)
correction_sq = G_approx - S_inf
assert_true("A40: G(3/S_inf) - S_inf > 0", correction_sq > 0)
assert_true("A40b: relative correction < 1e-6", correction_sq/S_inf < mpf('1e-6'))

# A41: Analytical formula for correction: G(3/S) - S = 2S³/((3-2S)(3+S))
# Derivation: G(3/S) = 3S(3-S)/((3-2S)(3+S)); subtract S → 2S³/((3-2S)(3+S))
correction_formula = 2*S_inf**3/((3-2*S_inf)*(3+S_inf))
assert_close("A41: analytical correction formula", correction_sq, correction_formula, mpf('1e-50'))

# A42: lam_minus ≈ 1 (the unphysical root)
assert_true("A42: lam_minus in (0.999, 1.0)", mpf('0.999') < lam_minus < mpf('1.0'))

# A43: lam_plus * lam_minus ≈ 3/S_inf - 2 (from Vieta)
vieta_prod_check = 3/S_inf - 2
assert_close("A43: lam+ * lam- = 3/S - 2 (approx)", lam_plus*lam_minus, vieta_prod_check,
             mpf('1e-45'))

# A44: Full PSLQ confirmation: ratio has no simple integer relation
# We verify the symbolic ratio matches to 55 significant figures
ratio_recheck = E_e**4*(mu**2+3*E_mu)*(mu-E_e)/(mu*(mu-2*E_e)*(mu+E_e))
assert_close("A44: ratio symbolic form exact to 55 places",
             G_lam0/S_inf, ratio_recheck, mpf('1e-55'))

print()
print("Key numerical results:")
print(f"  S_inf           = {nstr(S_inf, 20)}")
print(f"  G(mu/pi)/S_inf  = {nstr(ratio, 20)}")
print(f"  lambda★₊        = {nstr(lam_plus, 20)}")
print(f"  3/S_inf         = {nstr(approx_leading, 20)}")
print(f"  G(3/S_inf)/S    = {nstr(ratio_approx, 20)}")
print(f"  Relative error  = {nstr(rel_err_leading, 10)}")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
if FAIL:
    raise SystemExit(1)