"""
verify_P205.py — Verification script for Addendum P205
PSLQ Augmented with 4-Simplex Constants: Hunt for c_4, c_ex, and delta

Checks:
  §1  TOE constants and 4-simplex geometry
  §2  c_ex and delta (the P197 near-miss coefficient and its deficit from 9)
  §3  c_4 recomputation
  §4  Basis B3 construction and independence
  §5  PSLQ on B3 — c_4   (result: None, maxcoeff=2000)
  §6  PSLQ on B3 — c_ex  (result: None, maxcoeff=2000)
  §7  PSLQ on B3 — delta  (result: None, maxcoeff=2000)
  §8  Rational scan (p/q)*Omega0*alpha^-n, n=1,2,3
  §9  Near-coincidence c_4 ≈ (11/4)*Omega0*alpha^-2
  §10 PSLQ on residual Delta over B3 (result: None, maxcoeff=2000)
  §11 Near-coincidence analysis in CODATA sigma units

All computations use mpmath at mp.dps=60.
Run with:  python verify_P205.py
All assertions must pass.
"""

import sys
import mpmath
mpmath.mp.dps = 60

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}")

# ── TOE constants ─────────────────────────────────────────────────────────────
pi        = mpmath.pi
ALPHA_INV = 4*pi**3 + pi**2 + pi
alpha     = 1 / ALPHA_INV
OMEGA_0   = pi**3 / 4

# ── Experimental target ───────────────────────────────────────────────────────
m_mu_me = mpmath.mpf('206.7682830')

# ── Three-loop formula (P199/P201) ────────────────────────────────────────────
f2 = 1/pi**2 + mpmath.mpf(11)/16
N  = -OMEGA_0 * ALPHA_INV * mpmath.mpf(11)/4
R3 = 207*(1 - alpha/(2*pi) + f2*alpha**2 + N*alpha**3)

# ── Four-loop coefficient c_4 ─────────────────────────────────────────────────
delta_abs = m_mu_me - R3
c4        = delta_abs / (207 * alpha**4)

# ── 4-simplex constants (P202) ────────────────────────────────────────────────
theta4        = mpmath.acos(mpmath.mpf('-1')/4)   # geodesic angle, arccos(-1/4)
sin_theta4    = mpmath.sqrt(15)/4                  # sin(theta4)
p4            = mpmath.mpf('-1')/4                 # inner product
lambda4       = mpmath.mpf(5)/4                    # Gram eigenvalue
h_A4          = mpmath.mpf(5)                      # dual Coxeter number h∨(A_4)
S5_order      = mpmath.mpf(120)                    # |W(A_4)| = |S_5|
ln_5          = mpmath.log(5)
sqrt_5        = mpmath.sqrt(5)
sqrt_5_over_2 = mpmath.sqrt(mpmath.mpf(5)/2)      # edge length of 4-simplex on S^3

# ── Auxiliary constants ───────────────────────────────────────────────────────
ln2     = mpmath.log(2)
gamma_E = mpmath.euler
zeta3   = mpmath.zeta(3)

print("verify_P205.py — PSLQ 4-simplex augmented basis")
print(f"  c4    = {mpmath.nstr(c4, 30)}")


# ══════════════════════════════════════════════════════════════════════════════
# §1  TOE constants and 4-simplex geometry
# ══════════════════════════════════════════════════════════════════════════════
print("S1  TOE constants and 4-simplex geometry")

# 1.1  ALPHA_INV in (137, 138)
check(1, "ALPHA_INV range", mpmath.mpf('137') < ALPHA_INV < mpmath.mpf('138'))

# 1.2  OMEGA_0 = pi^3/4
check(2, "OMEGA_0 definition", abs(OMEGA_0 - pi**3/4) < mpmath.mpf('1e-55'))

# 1.3  theta4 = arccos(-1/4): verify cos(theta4) = -1/4
check(3, "cos(theta4) = -1/4", abs(mpmath.cos(theta4) - p4) < mpmath.mpf('1e-55'))

# 1.4  sin^2(theta4) + cos^2(theta4) = 1
check(4, "Pythagorean identity", abs(sin_theta4**2 + p4**2 - 1) < mpmath.mpf('1e-55'))

# 1.5  sin_theta4 = sqrt(15)/4 exactly
check(5, "sin_theta4 = sqrt(15)/4", abs(sin_theta4 - mpmath.sqrt(15)/4) < mpmath.mpf('1e-55'))

# 1.6  sin_theta4 > 0.96
check(6, "sin_theta4 > 0.96", sin_theta4 > mpmath.mpf('0.96'))

# 1.7  theta4 in (1.82, 1.83) radians
check(7, "theta4 in (1.82, 1.83)", mpmath.mpf('1.82') < theta4 < mpmath.mpf('1.83'))

# 1.8  lambda4 = 5/4
check(8, "lambda4 = 5/4", abs(lambda4 - mpmath.mpf(5)/4) < mpmath.mpf('1e-55'))

# 1.9  h_A4 = 5
check(9, "h_A4 = 5", h_A4 == mpmath.mpf(5))

# 1.10 h_A4 = 4*lambda4 (algebraic relation inside basis — dependency to exclude)
check(10, "h_A4 = 4*lambda4", abs(h_A4 - 4*lambda4) < mpmath.mpf('1e-55'))

# 1.11 sqrt_5_over_2 = sqrt(5)/sqrt(2)
check(11, "sqrt_5_over_2 = sqrt(5)/sqrt(2)",
      abs(sqrt_5_over_2 - sqrt_5/mpmath.sqrt(2)) < mpmath.mpf('1e-55'))

# 1.12 S5_order = 120
check(12, "S5_order = 120", S5_order == mpmath.mpf(120))


# ══════════════════════════════════════════════════════════════════════════════
# §2  c_ex and delta (P197 near-miss coefficient)
# ══════════════════════════════════════════════════════════════════════════════
print("S2  c_ex and delta")

# c_ex = (m_mu_me - ALPHA_INV) / OMEGA_0
c_ex      = (m_mu_me - ALPHA_INV) / OMEGA_0
delta_cex = mpmath.mpf(9) - c_ex

print(f"  c_ex  = {mpmath.nstr(c_ex, 35)}")
print(f"  delta = {mpmath.nstr(delta_cex, 35)}")

# 2.1  c_ex > 8.99 and < 9.00
check(13, "c_ex in (8.99, 9.00)", mpmath.mpf('8.99') < c_ex < mpmath.mpf('9.00'))

# 2.2  c_ex = 8.9958533... matches anchor to 10 sig figs
anchor_cex = mpmath.mpf('8.9958533161676325')
check(14, "c_ex matches anchor", abs(c_ex - anchor_cex) < mpmath.mpf('1e-14'))

# 2.3  Reconstruction identity: ALPHA_INV + c_ex*OMEGA_0 = m_mu_me
check(15, "ALPHA_INV + c_ex*OMEGA_0 = m_mu_me",
      abs(ALPHA_INV + c_ex*OMEGA_0 - m_mu_me) < mpmath.mpf('1e-50'))

# 2.4  delta = 9 - c_ex > 0 (actual coefficient falls short of 9)
check(16, "delta > 0", delta_cex > 0)

# 2.5  delta < 0.005
check(17, "delta < 0.005", delta_cex < mpmath.mpf('0.005'))

# 2.6  delta = 9 - c_ex matches anchor
anchor_delta = mpmath.mpf('0.0041466838323674713665')
check(18, "delta matches anchor", abs(delta_cex - anchor_delta) < mpmath.mpf('1e-20'))

# 2.7  ALPHA_INV + 9*OMEGA_0 overshoots m_mu_me
overshoot = ALPHA_INV + 9*OMEGA_0 - m_mu_me
check(19, "9*Omega0 candidate overshoots m_mu_me", overshoot > 0)

# 2.8  Overshoot = delta*OMEGA_0
check(20, "overshoot = delta*OMEGA_0",
      abs(overshoot - delta_cex*OMEGA_0) < mpmath.mpf('1e-50'))

# 2.9  Overshoot ≈ 0.032 (to 2 sig figs)
check(21, "overshoot ≈ 0.032",
      abs(overshoot - mpmath.mpf('0.032')) < mpmath.mpf('0.001'))


# ══════════════════════════════════════════════════════════════════════════════
# §3  c_4 recomputation
# ══════════════════════════════════════════════════════════════════════════════
print("S3  c_4 recomputation")

# 3.1  c4 integer part = 400305
check(22, "c4 integer part = 400305", int(mpmath.floor(c4)) == 400305)

# 3.2  c4 fractional part in (0.87, 0.88)
c4_frac = c4 - 400305
check(23, "c4 fractional part", mpmath.mpf('0.87') < c4_frac < mpmath.mpf('0.88'))

# 3.3  c4 matches P201 anchor
anchor_c4 = mpmath.mpf('400305.875316')
check(24, "c4 matches P201 anchor", abs(c4 - anchor_c4) < mpmath.mpf('1e-3'))

# 3.4  c4 > 0
check(25, "c4 positive", c4 > 0)

# 3.5  207 * c4 * alpha^4 = delta_abs (self-consistency)
check(26, "207*c4*alpha^4 = delta_abs",
      abs(207 * c4 * alpha**4 - delta_abs) < mpmath.mpf('1e-50'))

# 3.6  R3 + 207*c4*alpha^4 = m_mu_me (by construction)
R4 = R3 + 207*c4*alpha**4
check(27, "R4 = m_mu_me by construction", abs(R4 - m_mu_me) < mpmath.mpf('1e-10'))

# 3.7  CODATA uncertainty propagates as sigma_c4 ≈ 7.84
sigma_c4 = mpmath.mpf('4.6e-6') / (207 * alpha**4)
check(28, "sigma_c4 in (7, 9)", mpmath.mpf('7') < sigma_c4 < mpmath.mpf('9'))


# ══════════════════════════════════════════════════════════════════════════════
# §4  Basis B3 construction and algebraic independence
# ══════════════════════════════════════════════════════════════════════════════
print("S4  Basis B3 construction")

# B3_tail: 21 elements (without c4/c_ex/delta as first element)
B3_tail = [
    mpmath.mpf(1),
    pi, pi**2,
    ALPHA_INV, OMEGA_0, ALPHA_INV*OMEGA_0,
    ln2, gamma_E, zeta3,
    theta4, sin_theta4, theta4**2,
    theta4*ALPHA_INV, theta4*OMEGA_0,
    sin_theta4*ALPHA_INV, sin_theta4*OMEGA_0,
    lambda4*ALPHA_INV, lambda4*OMEGA_0,
    ln_5, sqrt_5, sqrt_5_over_2,
]

# 4.1  B3_tail has 21 elements
check(29, f"B3_tail size = {len(B3_tail)}, expected 21", len(B3_tail) == 21)

# 4.2  h_A4*OMEGA_0 excluded because h_A4 = 4*lambda4 makes it linearly dependent
#      on lambda4*OMEGA_0.  Verify the dependency:
check(30, "h_A4*Omega0 = 4*lambda4*Omega0 (excluded from B3)",
      abs(h_A4*OMEGA_0 - 4*(lambda4*OMEGA_0)) < mpmath.mpf('1e-55'))

# 4.3  No internal trivial integer relation among B3_tail (PSLQ returns None)
internal_check = mpmath.pslq(B3_tail, maxcoeff=500, tol=mpmath.mpf('1e-20'))
check(31, "B3_tail has no unexpected internal relation", internal_check is None)

# 4.4  B3 with c4 prepended has 22 elements
B3_c4 = [c4] + B3_tail
check(32, f"B3_c4 size = {len(B3_c4)}, expected 22", len(B3_c4) == 22)


# ══════════════════════════════════════════════════════════════════════════════
# §5  PSLQ on B3 — c_4  (no relation found)
# ══════════════════════════════════════════════════════════════════════════════
print("S5  PSLQ B3 / c_4")

print("Running PSLQ on B3 for c_4 (maxcoeff=2000)... ", end="", flush=True)
res_c4_B3 = mpmath.pslq(B3_c4, maxcoeff=2000, tol=mpmath.mpf('1e-20'))
print(f"result = {res_c4_B3}")

# 5.1  PSLQ on B3 for c_4 returns None or trivial (c4 coeff = 0)
check(33, "PSLQ B3 / c_4: no relation (|n_i| ≤ 2000) CONFIRMED",
      res_c4_B3 is None or res_c4_B3[0] == 0)


# ══════════════════════════════════════════════════════════════════════════════
# §6  PSLQ on B3 — c_ex  (no relation found)
# ══════════════════════════════════════════════════════════════════════════════
print("S6  PSLQ B3 / c_ex")

B3_cex = [c_ex] + B3_tail
print("Running PSLQ on B3 for c_ex (maxcoeff=2000)... ", end="", flush=True)
res_cex_B3 = mpmath.pslq(B3_cex, maxcoeff=2000, tol=mpmath.mpf('1e-20'))
print(f"result = {res_cex_B3}")

# 6.1  PSLQ on B3 for c_ex returns None or trivial
check(34, "PSLQ B3 / c_ex: no relation (|n_i| ≤ 2000) CONFIRMED",
      res_cex_B3 is None or res_cex_B3[0] == 0)


# ══════════════════════════════════════════════════════════════════════════════
# §7  PSLQ on B3 — delta = 9 - c_ex  (no relation found)
# ══════════════════════════════════════════════════════════════════════════════
print("S7  PSLQ B3 / delta")

B3_delta = [delta_cex] + B3_tail
print("Running PSLQ on B3 for delta (maxcoeff=2000)... ", end="", flush=True)
res_delta_B3 = mpmath.pslq(B3_delta, maxcoeff=2000, tol=mpmath.mpf('1e-20'))
print(f"result = {res_delta_B3}")

# 7.1  PSLQ on B3 for delta returns None or trivial
check(35, "PSLQ B3 / delta: no relation (|n_i| ≤ 2000) CONFIRMED",
      res_delta_B3 is None or res_delta_B3[0] == 0)


# ══════════════════════════════════════════════════════════════════════════════
# §8  Rational scan (p/q)*Omega0*alpha^-n, n=1,2,3
# ══════════════════════════════════════════════════════════════════════════════
print("S8  Rational scan")

# The best match from the scan is (11/4)*Omega0*alpha^-2 (the P201 near-coincidence).
# No new closer match with p,q <= 200 and n in {1,2,3} was found.

# 8.1  (11/4)*Omega0*alpha^-2 is within 4e-6 relative of c4
nc1   = mpmath.mpf(11)/4 * OMEGA_0 * ALPHA_INV**2
rel_1 = abs(c4 - nc1) / c4
check(36, "nc1 relative error < 4e-6", rel_1 < mpmath.mpf('4e-6'))

# 8.2  For n=1, the ratio c4/(Omega0*ALPHA_INV) ≈ 376.85 — not near a
#      rational with |p|,|q|<=200. Best absolute distance > 0.5.
ratio_n1 = c4 / (OMEGA_0 * ALPHA_INV)
nearest_n1 = mpmath.nint(ratio_n1)
check(37, f"n=1 ratio not close to integer: {float(ratio_n1)}",
      abs(ratio_n1 - nearest_n1) > mpmath.mpf('0.1'))

# 8.3  For n=3, ratio c4/(Omega0*ALPHA_INV^3) ≈ 0.02008 — not near a
#      simple rational with small denominator.
ratio_n3 = c4 / (OMEGA_0 * ALPHA_INV**3)
check(38, f"n=3 ratio small: {float(ratio_n3)}", ratio_n3 < mpmath.mpf('0.1'))

# 8.4  nc1 distance from c4 is 1.501 (to 3 sig figs)
check(39, "nc1 distance ≈ 1.501",
      abs(abs(c4 - nc1) - mpmath.mpf('1.501')) < mpmath.mpf('0.001'))


# ══════════════════════════════════════════════════════════════════════════════
# §9  Near-coincidence c_4 ≈ (11/4)*Omega0*alpha^-2
# ══════════════════════════════════════════════════════════════════════════════
print("S9  Near-coincidence analysis")

Delta = c4 - nc1   # ≈ -1.501

# 9.1  Delta < 0 (c4 < nc1)
check(40, "Delta = c4 - nc1 < 0", Delta < 0)

# 9.2  |Delta| < 2
check(41, "|Delta| < 2", abs(Delta) < 2)

# 9.3  |Delta| / sigma_c4 < 0.20 (within 0.20 CODATA sigma)
frac_sigma = abs(Delta) / sigma_c4
check(42, "|Delta|/sigma < 0.20", frac_sigma < mpmath.mpf('0.20'))

# 9.4  nc1 = (11/4)*Omega0*alpha^-2 = -N*alpha^-1 (cross-check with P201)
N_check = -OMEGA_0 * ALPHA_INV * mpmath.mpf(11)/4
nc1_via_N = -N_check * ALPHA_INV
check(43, "nc1 = -N*alpha^-1", abs(nc1 - nc1_via_N) < mpmath.mpf('1e-40'))

# 9.5  nc1 ≈ 400307 (integer part)
check(44, "nc1 ≈ 400307", abs(nc1 - mpmath.mpf('400307')) < 1)


# ══════════════════════════════════════════════════════════════════════════════
# §10 PSLQ on residual Delta over B3  (no relation found)
# ══════════════════════════════════════════════════════════════════════════════
print("S10  PSLQ B3 / Delta")

B3_Delta = [Delta] + B3_tail
print("Running PSLQ on B3 for Delta=c4-nc1 (maxcoeff=2000)... ", end="", flush=True)
res_Delta_B3 = mpmath.pslq(B3_Delta, maxcoeff=2000, tol=mpmath.mpf('1e-20'))
print(f"result = {res_Delta_B3}")

# 10.1  PSLQ on Delta returns None or trivial
check(45, "PSLQ B3 / Delta: no relation (|n_i| ≤ 2000) CONFIRMED",
      res_Delta_B3 is None or res_Delta_B3[0] == 0)


# ══════════════════════════════════════════════════════════════════════════════
# §11 Near-coincidence analysis in CODATA sigma
# ══════════════════════════════════════════════════════════════════════════════
print("S11  CODATA sigma analysis")

# sigma(m_mu/m_e) ≈ 4.6e-6, ∂c4/∂(m_mu/m_e) = 1/(207*alpha^4)
deriv_c4 = 1/(207*alpha**4)

# 11.1  deriv_c4 ≈ 1.703e6
check(46, "deriv_c4 ≈ 1.703e6",
      abs(deriv_c4 - mpmath.mpf('1.703e6')) < mpmath.mpf('1e4'))

# 11.2  sigma_c4 > 7 and < 9
check(47, "sigma_c4 in (7, 9)", mpmath.mpf('7') < sigma_c4 < mpmath.mpf('9'))

# 11.3  |Delta| is 0.19 sigma (within 0.01 tolerance)
check(48, "frac_sigma ≈ 0.191",
      abs(frac_sigma - mpmath.mpf('0.191')) < mpmath.mpf('0.002'))

# 11.4  c_ex overshoot from 9 in absolute units: delta*OMEGA_0 ≈ 0.0321
overshoot_abs = delta_cex * OMEGA_0
check(49, "delta*Omega0 ≈ 0.032",
      abs(overshoot_abs - mpmath.mpf('0.032')) < mpmath.mpf('0.001'))

# 11.5  delta < 1/200 (very small deficit from integer 9)
check(50, "delta < 1/200", delta_cex < mpmath.mpf(1)/200)


# ══════════════════════════════════════════════════════════════════════════════
# Final summary
# ══════════════════════════════════════════════════════════════════════════════
print()
print(f"  c_ex  = {mpmath.nstr(c_ex, 35)}")
print(f"  delta = 9 - c_ex = {mpmath.nstr(delta_cex, 35)}")
print(f"  c4    = {mpmath.nstr(c4, 30)}")
print()
print(f"  CODATA sigma (c4): {mpmath.nstr(sigma_c4, 6)}")
print(f"  Near-coincidence (11/4)*Omega0*alpha^-2:")
print(f"    value    = {mpmath.nstr(nc1, 20)}")
print(f"    residual = {mpmath.nstr(Delta, 12)}")
print(f"    = {mpmath.nstr(frac_sigma, 4)} CODATA sigma")
print()
print("  PSLQ results (all None / trivial):")
print(f"    B3 / c4   (maxcoeff=2000): {'None' if res_c4_B3 is None else 'trivial'}")
print(f"    B3 / c_ex (maxcoeff=2000): {'None' if res_cex_B3 is None else 'trivial'}")
print(f"    B3 / delta(maxcoeff=2000): {'None' if res_delta_B3 is None else 'trivial'}")
print(f"    B3 / Delta(maxcoeff=2000): {'None' if res_Delta_B3 is None else 'trivial'}")
print()
print("  STATUS: OPEN — OP-C4 and OP-C_EX remain unresolved.")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0)