"""
verify_P210.py — Verification for Addendum P210
OP-3 Resolution Attempt: h∨(A₂) and the Lepton Face Correction

All assertions at mp.dps = 60.
© Léon Fernando Vlegels. MIT License.
"""

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

# ─── Constants ──────────────────────────────────────────────────────────────
ALPHA_INV = 4*pi**3 + pi**2 + pi       # monad μ
OMEGA_0   = pi**3 / 4                   # fundamental frequency ω
E_e       = pi                          # electron energy e
TARGET    = mpf('206.7682830')          # CODATA muon/electron mass ratio

# Lie-algebraic quantities for A₂
HV_A2     = mpf('3')                    # h∨(A₂) = dual Coxeter number
RANK_A2   = mpf('2')                    # rank(A₂)
DIM_A2    = mpf('8')                    # dim(A₂) = 8
WEYL_A2   = mpf('6')                    # |W(A₂)| = 3! = 6
ROOTS_A2  = mpf('6')                    # |Φ(A₂)| = 6
POS_ROOTS = mpf('3')                    # |Φ⁺(A₂)| = 3
COXETER_A2 = mpf('3')                   # h(A₂) = Coxeter number = 3 = h∨ for A₂

TOL = mpf('1e-50')   # tight tolerance for exact equalities
PPM = mpf('1e-6')

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

# ─── Section 1: Core constant relations ─────────────────────────────────────

# A1: μ = 4π³ + π² + π
check(1, "A1: μ definition", fabs(ALPHA_INV - (4*pi**3 + pi**2 + pi)) < TOL)

# A2: ω = π³/4 exactly
check(2, "A2: ω = π³/4", fabs(OMEGA_0 - pi**3 / 4) < TOL)

# A3: e = π
check(3, "A3: e = π", fabs(E_e - pi) < TOL)

# A4: e³ = 4ω (exact)
check(4, "A4: e³ = 4ω", fabs(E_e**3 - 4*OMEGA_0) < TOL)

# A5: μ = 4e³ + e² + e (monad polynomial in e)
check(5, "A5: μ = 4e³+e²+e", fabs(ALPHA_INV - (4*E_e**3 + E_e**2 + E_e)) < TOL)

# A6: ω/e³ = 1/4 exactly
check(6, "A6: ω/e³ = 1/4", fabs(OMEGA_0 / E_e**3 - mpf('1')/4) < TOL)

# ─── Section 2: R★ computation ──────────────────────────────────────────────

R_star = 207 * (1 + (HV_A2 - OMEGA_0) / (E_e**3 * ALPHA_INV))
gap_star_ppm = fabs(R_star - TARGET) / TARGET * 1e6

# A7: R* agrees with CODATA within 2 ppm
check(7, f"A7: R* agrees with CODATA within 2 ppm (gap {float(gap_star_ppm):.4f} ppm)", gap_star_ppm < 2)

# A8: R* agrees within 1.2 ppm (tight)
check(8, f"A8: R* agrees within 1.2 ppm (gap {float(gap_star_ppm):.6f} ppm)", gap_star_ppm < mpf('1.2'))

# A9: R* is close to but strictly greater than TARGET
check(9, "A9: R* > TARGET", R_star > TARGET)

# A10: Equivalent sign-flipped form gives same result
R_star2 = 207 * (1 - (OMEGA_0 - HV_A2) / (E_e**3 * ALPHA_INV))
check(10, "A10: sign-flip equivalence", fabs(R_star - R_star2) < TOL)

# A11: R* is in the ball (TARGET-1ppm, TARGET+5ppm)
check(11, "A11: R* within 5 ppm of TARGET", fabs(R_star - TARGET) / TARGET < 5*PPM)

# ─── Section 2: Geometric quantities ────────────────────────────────────────

diff = OMEGA_0 - HV_A2  # ω - h∨

# A12: ω - h∨ = π³/4 - 3
check(12, "A12: ω - h∨ = π³/4 - 3", fabs(diff - (pi**3/4 - 3)) < TOL)

# A13: ω - h∨ is positive (ω ≈ 7.75 > 3)
check(13, "A13: ω > h∨", diff > 0)

# A14: 4(ω - h∨) = π³ - 12
check(14, "A14: 4(ω-h∨) = π³-12", fabs(4*diff - (pi**3 - 12)) < TOL)

# A15: 4(ω - h∨) = e³ - 4h∨
check(15, "A15: 4(ω-h∨) = e³-4h∨", fabs(4*diff - (E_e**3 - 4*HV_A2)) < TOL)

# A16: Correction reformulation — (ω-h∨)/(e³μ) = (1-h∨/ω)/(4μ)
lhs = diff / (E_e**3 * ALPHA_INV)
rhs = (1 - HV_A2/OMEGA_0) / (4*ALPHA_INV)
check(16, "A16: correction reformulation", fabs(lhs - rhs) < TOL)

# A17: h∨/ω = 12/π³
check(17, "A17: h∨/ω = 12/π³", fabs(HV_A2/OMEGA_0 - 12/pi**3) < TOL)

# A18: h∨/ω < 1 (defect is positive)
check(18, "A18: h∨ < ω", HV_A2/OMEGA_0 < 1)

# A19: h∨/ω is approximately 0.387
check(19, "A19: h∨/ω ≈ 0.387", fabs(HV_A2/OMEGA_0 - mpf('0.387')) < mpf('0.001'))

# A20: defect δ = 1 - h∨/ω ≈ 0.613
delta = 1 - HV_A2/OMEGA_0
check(20, "A20: δ ≈ 0.613", fabs(delta - mpf('0.613')) < mpf('0.001'))

# A21: The correction (ω-h∨)/(e³μ) is of order 10⁻³
corr = diff / (E_e**3 * ALPHA_INV)
check(21, "A21: correction order 10⁻³", corr > mpf('1e-4') and corr < mpf('1e-2'))

# A22: Correction × 10⁶ ≈ 1118 (in ppm units, unnormalised)
check(22, "A22: correction ≈ 1118 × 10⁻⁶", fabs(corr * 1e6 - 1118) < 1)

# ─── Section 3: Gram matrix of lepton face ──────────────────────────────────

# Gram matrix: G_ii = 1, G_ij = -1/4
a = mpf('1')
b = mpf('-1') / 4

# A23: Eigenvalue λ_max = 5/4 (multiplicity 2)
ev_max = a - b   # = 1 + 1/4 = 5/4
check(23, "A23: λ_max = 5/4", fabs(ev_max - mpf('5')/4) < TOL)

# A24: Eigenvalue λ_min = 1/2 (multiplicity 1)
ev_min = a + 2*b   # = 1 - 1/2 = 1/2
check(24, "A24: λ_min = 1/2", fabs(ev_min - mpf('1')/2) < TOL)

# A25: Trace = 2*(5/4) + 1/2 = 3 = h∨(A₂)
trace = 2*ev_max + ev_min
check(25, "A25: tr(G) = h∨(A₂) = 3", fabs(trace - HV_A2) < TOL)

# A26: det(G) = (5/4)² * (1/2) = 25/32
det_G = ev_max**2 * ev_min
check(26, "A26: det(G) = 25/32", fabs(det_G - mpf('25')/32) < TOL)

# A27: All eigenvalues are positive (Gram matrix is positive definite)
check(27, "A27: all eigenvalues positive (λ_max > 0, λ_min > 0)", ev_max > 0 and ev_min > 0)

# A28: λ_max / λ_min = 5/2
check(28, "A28: λ_max/λ_min = 5/2", fabs(ev_max/ev_min - mpf('5')/2) < TOL)

# A29: off-diagonal = inner product of lepton face vertices = -1/4
check(29, "A29: G_ij = -1/4", fabs(b - mpf('-1')/4) < TOL)

# A30: Euler characteristic of boundary triangle: V - E = 0
V_tri, E_tri = 3, 3
chi_boundary = V_tri - E_tri
check(30, "A30: χ_boundary = V-E = 0", chi_boundary == 0)

# ─── Section 3: Integer scan ─────────────────────────────────────────────────

def R_formula(n):
    return 207 * (1 + (n - OMEGA_0) / (E_e**3 * ALPHA_INV))

def gap_ppm(n):
    return fabs(R_formula(n) - TARGET) / TARGET * 1e6

# A31: n=3 gives the minimum gap in n=1..10
gaps = {n: gap_ppm(n) for n in range(1, 11)}
check(31, "A31: n=3 gives the minimum gap in n=1..10",
      all(gaps[3] < gaps[n] for n in range(1, 11) if n != 3))

# A32: n=3 gap < 2 ppm
check(32, "A32: n=3 gap < 2 ppm", gaps[3] < 2)

# A33: n=3 gap < n=2 gap (h∨ better than rank)
check(33, "A33: h∨ better than rank", gaps[3] < gaps[2])

# A34: n=3 gap < n=6 gap (h∨ better than Weyl group order)
check(34, "A34: h∨ better than |W(A₂)|", gaps[3] < gaps[6])

# A35: n=3 gap < n=8 gap (h∨ better than dim(A₂))
check(35, "A35: h∨ better than dim(A₂)", gaps[3] < gaps[8])

# A36: n=3 gap is at least 100× better than n=4
check(36, "A36: n=3 at least 100× better than n=4", gaps[3] < gaps[4] / 100)

# A37: n=3 gap is at least 100× better than n=2
check(37, "A37: n=3 at least 100× better than n=2", gaps[3] < gaps[2] / 100)

# A38: Gaps increase monotonically away from n=3 for n=3..10
check(38, "A38: gaps increase monotonically away from n=3 for n=3..10",
      all(gaps[n] < gaps[n+1] for n in range(3, 10)))

# A39: Gaps decrease monotonically approaching n=3 from below (n=1,2,3)
check(39, "A39: gaps decrease monotonically approaching n=3 from below (gap(1) > gap(2) > gap(3))",
      gaps[1] > gaps[2] and gaps[2] > gaps[3])

# A40: R(n=3) is between 206.768 and 206.769
R3 = R_formula(3)
check(40, f"A40: R(3) = {float(R3):.6f} in (206.768, 206.769)",
      R3 > mpf('206.768') and R3 < mpf('206.769'))

# ─── Section 4: A₂ data consistency ─────────────────────────────────────────

# A41: h∨(A₂) = h(A₂) = |Φ⁺(A₂)| = 3 (triple coincidence for A₂)
check(41, "A41: h∨(A₂) = h(A₂) = |Φ⁺(A₂)| = 3 (triple coincidence)",
      fabs(HV_A2 - COXETER_A2) < TOL and fabs(HV_A2 - POS_ROOTS) < TOL)

# A42: dim(A₂) = rank + 2|Φ⁺| = 2 + 6 = 8
check(42, "A42: dim = rank + 2|Φ⁺|", fabs(DIM_A2 - (RANK_A2 + 2*POS_ROOTS)) < TOL)

# A43: |W(A₂)| = h∨(A₂)! = 6
import math
check(43, "A43: |W| = h∨! = 6", fabs(WEYL_A2 - math.factorial(int(HV_A2))) < TOL)

# A44: |Φ(A₂)| = 2|Φ⁺(A₂)| = 6
check(44, "A44: |Φ| = 2|Φ⁺|", fabs(ROOTS_A2 - 2*POS_ROOTS) < TOL)

# A45: The heat budget sums to μ
E_Fork    = 4*pi**3
E_Weld    = pi**2
E_Plateau = pi
heat_total = E_Fork + E_Weld + E_Plateau
check(45, "A45: E_F+E_W+E_P = μ", fabs(heat_total - ALPHA_INV) < TOL)

# A46: E_Fork = 4e³ = 16ω
check(46, "A46: E_F = 4e³ = 16ω",
      fabs(E_Fork - 4*E_e**3) < TOL and fabs(E_Fork - 16*OMEGA_0) < TOL)

# A47: E_Weld = e² = 4·(e²/4) — just checks the value
check(47, "A47: E_W = e²", fabs(E_Weld - E_e**2) < TOL)

# A48: E_Plateau = e
check(48, "A48: E_P = e", fabs(E_Plateau - E_e) < TOL)

# A49: Fork dominates: E_Fork > 90% of total
check(49, "A49: E_Fork > 90% of μ", E_Fork / heat_total > mpf('0.9'))

# A50: The full correction is of order 1/(4μ)
correction_full = (OMEGA_0 - HV_A2) / (E_e**3 * ALPHA_INV)
one_over_4mu = 1 / (4*ALPHA_INV)
# correction_full = one_over_4mu * delta
check(50, "A50: correction = δ/(4μ)",
      fabs(correction_full - one_over_4mu * delta) < TOL)

# ─── Final report ──────────────────────────────────────────────────────────

print(f"R★  = {nstr(R_star, 20)}")
print(f"Gap = {nstr(gap_star_ppm, 8)} ppm")
print(f"ω   = {nstr(OMEGA_0, 20)}")
print(f"h∨(A₂) = 3")
print(f"ω − h∨ = {nstr(diff, 20)}")
print(f"Defect δ = 1 − h∨/ω = {nstr(delta, 15)}")
print(f"Correction = {nstr(correction_full, 15)}")
print(f"(mp.dps={mp.dps})")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
import sys
sys.exit(0 if FAIL == 0 else 1)