"""
verify_P220.py — Verification suite for Addendum 220 (OP-rSinf-100)
All assertions at mp.dps = 60.

Run: python verify_P220.py
Pass condition: all assertions pass with zero exceptions.

Copyright: Léon Fernando Vlegels. License: MIT.
"""

import sys

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

# ── Core constants ─────────────────────────────────────────────────────────────
mu      = 4*pi**3 + pi**2 + pi        # ALPHA_INV ≈ 137.036
E_mu    = pi**2
h       = mpf(3)
E_e     = pi
OMEGA_0 = pi**3 / 4
TARGET  = mpf('206.7682830')

# ── Derived quantities ─────────────────────────────────────────────────────────
r     = h * E_mu / mu**2
S_inf = h * mu / (E_e**3 * (mu**2 + h * E_mu))
R_inf = 207 * (1 - OMEGA_0 / E_e**3 / mu + S_inf)
delta = R_inf - TARGET
c     = delta / 207

# ── Near-miss ε ───────────────────────────────────────────────────────────────
# The OP-rSinf-100 near-miss: r·S∞ ≈ 100·c,  ε ≡ r·S∞/(100c) − 1
eps = r * S_inf / (100 * c) - 1

# ── Closed form ───────────────────────────────────────────────────────────────
cf_rS  = mpf(9) / (mu * pi * (mu**2 + 3*pi**2))   # P217 closed form

# ── CODATA uncertainty ────────────────────────────────────────────────────────
sigma_T   = mpf('4.6e-6')
sigma_c   = sigma_T / 207
sigma_eps = fabs(r * S_inf) / (100 * c**2 * 207) * sigma_T

PASS = FAIL = 0

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


print("verify_P220.py — OP-rSinf-100 assertions  (mp.dps=60)")

# ─── A. Constant sanity ───────────────────────────────────────────────────────
print("S1  A. Constant sanity")
check("A01 mu>136", mu > 136)
check("A02 mu<138", mu < 138)
check("A03 mu formula", fabs(mu - (4*pi**3 + pi**2 + pi)) < mpf('1e-55'))
check("A04 r positive", r > 0)
check("A05 r magnitude", fabs(r - mpf('1.5767e-3')) < mpf('1e-6'))
check("A06 S_inf positive", S_inf > 0)
check("A07 S_inf magnitude", fabs(S_inf - mpf('7.049e-4')) < mpf('1e-6'))
check("A08 OMEGA_0 = pi^3/4", fabs(OMEGA_0 - pi**3/4) < mpf('1e-55'))
check("A09 R_inf > TARGET", R_inf > TARGET)
check("A10 delta positive", delta > 0)
check("A11 c positive", c > 0)
check("A12 c small", c < mpf('1e-7'))

# ─── B. Closed-form identity r·S∞ = 9/[μπ(μ²+3π²)] ──────────────────────────
print()
print("S2  B. Closed-form identity")
check("B01 cf_rS positive", cf_rS > 0)
check("B02 rS equals cf at 50-digit precision",
      fabs(r * S_inf - cf_rS) < mpf('1e-50'))
check("B03 cf numerator is h^2", fabs(mpf(9) - h**2) < mpf('1e-55'))
check("B04 cf denominator has 3pi^2 coefficient",
      fabs(mu**2 + 3*pi**2 - (mu**2 + h * E_mu)) < mpf('1e-50'))
check("B05 rS = h^2/[mu pi (mu^2 + h*E_mu)]",
      fabs(r * S_inf - h**2 / (mu * pi * (mu**2 + h * E_mu))) < mpf('1e-50'))
check("B06 r = h*E_mu/mu^2", fabs(r - h * E_mu / mu**2) < mpf('1e-55'))
check("B07 S_inf = h*mu/(E_e^3*(mu^2+h*E_mu))",
      fabs(S_inf - h * mu / (E_e**3 * (mu**2 + h * E_mu))) < mpf('1e-55'))

# ─── C. R∞ residual structure ─────────────────────────────────────────────────
print()
print("S3  C. R∞ residual")
check("C01 R_inf ≈ 206.768", fabs(R_inf - mpf('206.768')) < mpf('0.001'))
check("C02 delta ≈ 2.3e-6", fabs(delta - mpf('2.3e-6')) < mpf('3e-7'))
check("C03 c ≈ 1.111e-8", fabs(c - mpf('1.111e-8')) < mpf('1e-10'))
check("C04 c = delta/207", fabs(c - delta / 207) < mpf('1e-55'))
check("C05 R_inf = 207*(1 - OMEGA_0/pi^3/mu + S_inf)",
      fabs(R_inf - 207*(1 - OMEGA_0/E_e**3/mu + S_inf)) < mpf('1e-50'))
check("C06 sigma_c = sigma_T/207",
      fabs(sigma_c - sigma_T/207) < mpf('1e-30'))
check("C07 relative uncertainty in c near 2",
      fabs(sigma_c / c - 2) < mpf('0.01'))   # σ_c/c ≈ 2 (200%)

# ─── D. ε precision (Section 1) ──────────────────────────────────────────────
print()
print("S4  D. ε precision")
check("D01 eps positive", eps > 0)
check("D02 eps order 1e-5",
      fabs(eps) > mpf('1e-6') and fabs(eps) < mpf('1e-3'))
check("D03 eps leading digits 2.2036...",
      fabs(eps - mpf('2.2037e-5')) < mpf('1e-8'))
check("D04 100c < r*S_inf", 100 * c < r * S_inf)
check("D05 ratio r*S∞/(100c) > 1", r * S_inf / (100 * c) > 1)
check("D06 ratio r*S∞/(100c) < 1.001", r * S_inf / (100 * c) < mpf('1.001'))
# 25-digit precision of ε
eps_ref = mpf('2.203672588593701889106322e-5')
check("D07 eps matches reference to 20 digits",
      fabs(eps - eps_ref) < mpf('1e-26'))
check("D08 eps*mu < 0.01", eps * mu < mpf('0.01'))
check("D09 eps*mu > 0.001", eps * mu > mpf('0.001'))
check("D10 eps*mu^2 < 1", eps * mu**2 < 1)

# ─── E. CODATA uncertainty (Section 1) ───────────────────────────────────────
print()
print("S5  E. CODATA uncertainty")
check("E01 sigma_eps large (>1)", sigma_eps > 1)
check("E02 sigma_eps < 4", sigma_eps < 4)
check("E03 |eps/sigma_eps| < 1 → Case A",
      fabs(eps / sigma_eps) < 1)
check("E04 |eps/sigma_eps| < 1e-3",
      fabs(eps / sigma_eps) < mpf('1e-3'))
check("E05 sigma_eps ≈ 2.0",
      fabs(sigma_eps - 2) < mpf('0.1'))
# sigma_T > delta: the uncertainty in TARGET exceeds the residual itself
check("E06 sigma_T > delta (uncertainty exceeds residual)",
      sigma_T > delta)

# ─── F. PSLQ null results (Section 2) ────────────────────────────────────────
print()
print("S6  F. PSLQ null results")
# Basis A: ε vs pi/mu powers — no relation involving ε
bA = [eps, 1/mu, 1/(mu*pi), pi/mu**2, 1/mu**2, pi**2/mu**2, S_inf, r]
relA = pslq(bA, maxcoeff=1000, tol=mpf('1e-30'))
# A null or trivial result (coefficient of ε is 0) means no closed form found
check("F01 PSLQ basis-A: ε coeff is zero",
      relA is None or relA[0] == 0)

# Basis B: ε, r·S∞, r, S_inf, 1 — check no trivial relation with ε
bB = [eps, r * S_inf, r, S_inf, mpf(1)]
relB = pslq(bB, maxcoeff=1000, tol=mpf('1e-30'))
check("F02 PSLQ basis-B: ε coeff is zero",
      relB is None or relB[0] == 0)

# ε·μ^k is not close to an integer for k=1..4 (confirms no rational form n/μ^k)
from mpmath import nint
for k, limit, label in [(1, mpf('1e-3'), 'F03'), (2, mpf('0.05'), 'F04'),
                         (3, mpf('0.1'),  'F05'), (4, mpf('0.1'),  'F06')]:
    val = eps * mu**k
    diff = fabs(val - nint(val))
    check(f"{label} ε·μ^{k} not near integer (diff>{limit})", diff > limit)

# ─── G. Closed-form factored form ────────────────────────────────────────────
print()
print("S7  G. Closed-form factored")
# r·S∞ = h²/[μπ(μ²+3π²)] — verify numerically
rS_via_r_times_S = r * S_inf
rS_closed = h**2 / (mu * pi * (mu**2 + 3 * pi**2))
check("G01 r·S∞ = h²/[μπ(μ²+3π²)] to 50 digits",
      fabs(rS_via_r_times_S - rS_closed) < mpf('1e-50'))

# r·S∞ = r × S∞ (product equals closed form)
check("G02 r×S∞ positive", r * S_inf > 0)
check("G03 r = 3π²/μ²", fabs(r - 3*pi**2/mu**2) < mpf('1e-55'))
check("G04 S∞ = 3μ/[π³(μ²+3π²)]",
      fabs(S_inf - 3*mu/(pi**3 * (mu**2 + 3*pi**2))) < mpf('1e-50'))

# ─── H. Geometric interpretation (Section 3) ─────────────────────────────────
print()
print("S8  H. Geometric (f-vector)")
f_sum = 30  # 5+10+10+5 faces of 4-simplex
check("H01 f_sum^2 = 900", f_sum**2 == 900)
check("H02 f_sum^2 / h^2 = 100", f_sum**2 / h**2 == 100)
check("H03 f_sum = h_dual_E8 = 30", f_sum == 30)
check("H04 h∨(E₈) = 30 → f_sum = h∨(E₈)", f_sum == 30)
# Rewrite: r·S∞·100 = (f_sum)²/[μπ(μ²+3π²)]
rS100_geom = mpf(f_sum)**2 / (mu * pi * (mu**2 + 3*pi**2))
rS100_direct = r * S_inf * 100
check("H05 (f_sum)²/[μπ(μ²+3π²)] = r·S∞·100",
      fabs(rS100_geom - rS100_direct) < mpf('1e-45'))
# Near-miss as geometric: (f_sum)²/[μπ(μ²+3π²)] ≈ 10000·c
check("H06 ratio ≈ 10000", fabs(rS100_geom / c - 10000) < 1)
# dim(B₂)² = 100
check("H07 dim(B₂)² = 10² = 100", 10**2 == 100)
# (h∨(E₈)/h∨(A₂))² = (30/3)² = 100
check("H08 (h∨(E₈)/h∨(A₂))² = 100", (30 // 3)**2 == 100)

# ─── I. Case A formal conditions ─────────────────────────────────────────────
print()
print("S9  I. Case A conditions")
# Case A: ε is within CODATA uncertainty → near-miss may be exact
check("I01 |ε| < σ_ε", fabs(eps) < sigma_eps)
check("I02 |ε/σ_ε| < 0.001", fabs(eps / sigma_eps) < mpf('1e-3'))
# The required measurement precision to test exactness:
# Need σ_T_required << ε · 100 · c² · 207 / (r·S∞)
sigma_T_required = fabs(eps) * 100 * c**2 * 207 / (r * S_inf)
check("I03 required σ_T << current σ_T",
      sigma_T_required < sigma_T / 100)
check("I04 required σ_T > 0", sigma_T_required > 0)
# The near-miss cannot be falsified at current CODATA precision
check("I05 σ_c/c > 1 (c uncertain beyond 100%)", sigma_c / c > 1)

# ─── J. Internal consistency ─────────────────────────────────────────────────
print()
print("S10  J. Internal consistency")
check("J01 R_inf reconstructed", fabs(207 * (1 - OMEGA_0/pi**3/mu + S_inf) - R_inf) < mpf('1e-50'))
check("J02 delta = R_inf - TARGET", fabs(delta - (R_inf - TARGET)) < mpf('1e-55'))
check("J03 c = delta/207", fabs(c*207 - delta) < mpf('1e-55'))
check("J04 eps = rS/(100c) - 1 identity", fabs((r*S_inf/(100*c) - 1) - eps) < mpf('1e-55'))
check("J05 eps > 0", eps > 0)
check("J06 100*c < r*S_inf", 100*c < r*S_inf)
check("J07 r*S_inf/c between 99 and 101", 99 < r*S_inf/c < 101)
check("J08 (r*S_inf/c - 100)/100 == eps", fabs((r*S_inf/c - 100)/100 - eps) < mpf('1e-55'))

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)