#!/usr/bin/env python3
"""
verify_P232.py — Numerical companion to Addendum 232 (A232).

Topic: SR1 Closure via the Screening Axiom.

Verifies:
  C01. TOE constants: alpha_inv, kappa_P187, kappa_SA.
  C02. Yukawa ODE residual at kappa_P187 < 1e-50.
  C03. Yukawa ODE residual at kappa_SA   < 1e-50.
  C04. beta=2 stretched-exp ODE residual at kappa_P187 > 1e-3
       (evaluated at r ~ kappa_P187 where profiles differ).
  C05. Linear profile ODE residual at kappa_P187 > 1e-3.
  C06. Quadratic profile ODE residual at kappa_P187 > 1e-3.
  C07. Square-root profile ODE residual at kappa_P187 > 1e-3.
  C08. Sinusoidal profile ODE residual at kappa_P187 > 1e-3.
  C09. kappa_SA = alpha_inv/3 in (45, 46).
  C10. V_inf = E_self/mu_0^2.
  C11. Unit normalisation V_inf_SA * I = 1.0 to 1e-10.
  C12. Eigenvalue shift at kappa_SA < 1e-4 (well below S^3 gap of 5).
  C13. Differential shift Yukawa vs beta=2 at kappa_SA < 1e-5.
  C14. Fractional mass ratio change < 1e-4.
  C15. Yukawa eigenfunction: f'' = kappa^{-2} f, f(0)=1, f decays.
  C16. kappa_SA / kappa_P187 > 1e4 (distinct scales).
  C17. P187 bound V_inf*sqrt(kappa_P187) in (1e-5, 1e-4);
       delta_E1(kappa_SA) < P187 bound.

mp.dps = 60. All checks emit PASS/FAIL. Exits SystemExit(1) if FAIL > 0.

Copyright: Leon Fernando Vlegels. License: MIT. 2026-05-23.
"""

from mpmath import (mp, mpf, pi, exp, sqrt, fabs, quad, sin, nstr)
import sys

mp.dps = 60

SEP  = "=" * 72
SEP2 = "-" * 72

PASS_COUNT = 0
FAIL_COUNT = 0


def check(label, condition, detail=""):
    global PASS_COUNT, FAIL_COUNT
    n = PASS_COUNT + FAIL_COUNT + 1
    if condition:
        PASS_COUNT += 1
        print(f"  [PASS] {n:>2}. {label}")
    else:
        FAIL_COUNT += 1
        msg = f"  [FAIL] {n:>2}. {label}"
        if detail:
            msg += f"  [{detail}]"
        print(msg)


# ─────────────────────────────────────────────────────────────────────────────
print(SEP)
print("verify_P232 — SR1 Closure via the Screening Axiom (A232)")
print(SEP)

# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C01]  TOE CONSTANTS")
print(SEP2)

ALPHA_INV = 4*pi**3 + pi**2 + pi     # mu_0 = alpha^{-1}  (P03)
ALPHA     = mpf(1) / ALPHA_INV
E_SELF    = mpf("13.177")            # P03 Theorem 1
V_INF     = E_SELF / ALPHA_INV**2   # V_infinity = E_self / mu_0^2

# P187 screening length (saturation scale of stretched-exponential family)
KAPPA_P187 = ALPHA ** mpf("1.25")   # alpha^{5/4}

# SA screening length: kappa_SA = 1/(E_1*alpha) = alpha_inv/3
E1       = mpf(3)                   # S^3 eigenvalue at ell=1
KAPPA_SA = ALPHA_INV / E1           # = 1/(3*alpha)

MU1      = 16*pi**3/5 + 3*pi**2/4 + 2*pi/3   # int_0^1 r*rho(r) dr
MU_RATIO = MU1 / ALPHA_INV                    # ~ 0.7933

print(f"  alpha_inv  = mu_0 = {nstr(ALPHA_INV, 12)}  (4pi^3+pi^2+pi ~ 137.036)")
print(f"  alpha             = {nstr(ALPHA, 12)}")
print(f"  kappa_P187        = {nstr(KAPPA_P187, 8)}  (alpha^{{5/4}})")
print(f"  kappa_SA          = {nstr(KAPPA_SA, 8)}  (alpha_inv/3 ~ 45.7)")
print(f"  E_self            = {E_SELF}")
print(f"  V_inf             = {nstr(V_INF, 8)}")
print(f"  mu1/mu0           = {nstr(MU_RATIO, 8)}  (expected ~0.7933)")

check("C01a alpha_inv ~ 137.036",
      fabs(ALPHA_INV - mpf("137.036")) < mpf("0.001"))
check("C01b kappa_P187 in (2e-3, 3e-3)",
      mpf("2e-3") < KAPPA_P187 < mpf("3e-3"))
check("C01c kappa_SA in (45, 46)",
      mpf("45") < KAPPA_SA < mpf("46"))
check("C01d V_inf > 0", V_INF > 0)
check("C01e mu1/mu0 ~ 0.7933",
      fabs(MU_RATIO - mpf("0.7933")) < mpf("0.001"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C02-C03]  YUKAWA ODE RESIDUAL")
print("  V'' = kappa^{-2}(V - V_inf) must be satisfied exactly by Yukawa.")
print(SEP2)

def rho_toe(r):
    """P03 TOE density: 16*pi^3*r^3 + 3*pi^2*r^2 + 2*pi*r."""
    return 16*pi**3*r**3 + 3*pi**2*r**2 + 2*pi*r

def V_yukawa(r, kappa, Vinf=V_INF):
    return Vinf * (1 - exp(-r / kappa))

def V_yukawa_pp(r, kappa, Vinf=V_INF):
    """Analytic d²V/dr² for Yukawa."""
    return -Vinf / kappa**2 * exp(-r / kappa)

# Coarse grid: points spread across [0,1]
r_coarse = [mpf("0.01"), mpf("0.1"), mpf("0.3"), mpf("0.5"), mpf("0.7"), mpf("0.99")]

for lab, kappa in [("kappa_P187", KAPPA_P187), ("kappa_SA", KAPPA_SA)]:
    max_res = max(
        fabs(V_yukawa_pp(r, kappa) - (V_yukawa(r, kappa) - V_INF) / kappa**2)
        for r in r_coarse
    )
    print(f"  Yukawa ODE residual at {lab}: {nstr(max_res, 6)}")

check("C02 Yukawa ODE residual < 1e-50 at kappa_P187",
      max(fabs(V_yukawa_pp(r, KAPPA_P187) - (V_yukawa(r, KAPPA_P187) - V_INF) / KAPPA_P187**2)
          for r in r_coarse) < mpf("1e-50"))

check("C03 Yukawa ODE residual < 1e-50 at kappa_SA",
      max(fabs(V_yukawa_pp(r, KAPPA_SA) - (V_yukawa(r, KAPPA_SA) - V_INF) / KAPPA_SA**2)
          for r in r_coarse) < mpf("1e-50"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C04-C08]  ALTERNATIVE PROFILE ODE RESIDUALS")
print("  Evaluated with V_inf = V_INF and r-points spanning kappa_P187.")
print("  All alternatives must violate ODE by > 1e-3.")
print(SEP2)

kp = KAPPA_P187

# Points covering the saturation scale kappa_P187 and beyond
r_fine = [kp/2, kp, 2*kp, 5*kp, mpf("0.1"), mpf("0.5"), mpf("0.99")]

def ode_residual_vinf(V_func, V_pp_func, kappa, pts, Vinf=V_INF):
    """Max |V'' - kappa^{-2}(V - V_inf)| using the physical V_inf."""
    return max(fabs(V_pp_func(r) - (V_func(r) - Vinf) / kappa**2) for r in pts)

# --- C04: beta=2 stretched exponential ---
def V_b2(r):    return V_INF * (1 - exp(-(r/kp)**2))
def V_b2_pp(r):
    u = r / kp
    return V_INF * exp(-u**2) * (4*u**2 - 2) / kp**2

res_b2 = ode_residual_vinf(V_b2, V_b2_pp, kp, r_fine)
print(f"  beta=2 stretched-exp residual:      {nstr(res_b2, 6)}")
check("C04 beta=2 residual > 1e-3 (near kappa_P187)", res_b2 > mpf("1e-3"),
      f"got {nstr(res_b2, 4)}")

# --- C05: Linear V(r) = V_inf * r ---
def V_lin(r):    return V_INF * r
def V_lin_pp(r): return mpf(0)

res_lin = ode_residual_vinf(V_lin, V_lin_pp, kp, r_fine)
print(f"  Linear (V_inf*r) residual:          {nstr(res_lin, 6)}")
check("C05 linear residual > 1e-3", res_lin > mpf("1e-3"),
      f"got {nstr(res_lin, 4)}")

# --- C06: Quadratic V(r) = V_inf * r^2 ---
def V_qua(r):    return V_INF * r**2
def V_qua_pp(r): return 2 * V_INF

res_qua = ode_residual_vinf(V_qua, V_qua_pp, kp, r_fine)
print(f"  Quadratic (V_inf*r^2) residual:     {nstr(res_qua, 6)}")
check("C06 quadratic residual > 1e-3", res_qua > mpf("1e-3"),
      f"got {nstr(res_qua, 4)}")

# --- C07: Square-root V(r) = V_inf * sqrt(r), avoid r=0 ---
r_fine_pos = [kp, 2*kp, 5*kp, mpf("0.1"), mpf("0.5"), mpf("0.99")]
def V_sqr(r):    return V_INF * sqrt(r)
def V_sqr_pp(r): return V_INF * mpf("-0.25") * r**mpf("-1.5")

res_sqr = ode_residual_vinf(V_sqr, V_sqr_pp, kp, r_fine_pos)
print(f"  Square-root (V_inf*sqrt(r)) resid:  {nstr(res_sqr, 6)}")
check("C07 square-root residual > 1e-3", res_sqr > mpf("1e-3"),
      f"got {nstr(res_sqr, 4)}")

# --- C08: Sinusoidal V(r) = V_inf * sin(pi*r/2) ---
def V_sn(r):    return V_INF * sin(pi*r/2)
def V_sn_pp(r): return -V_INF * (pi/2)**2 * sin(pi*r/2)

res_sn = ode_residual_vinf(V_sn, V_sn_pp, kp, r_fine)
print(f"  Sinusoidal (V_inf*sin(pi*r/2)):     {nstr(res_sn, 6)}")
check("C08 sinusoidal residual > 1e-3", res_sn > mpf("1e-3"),
      f"got {nstr(res_sn, 4)}")


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C09]  kappa_SA = alpha_inv / 3")
print(SEP2)

print(f"  kappa_SA = {nstr(KAPPA_SA, 10)}  (expected ~45.679)")
check("C09 kappa_SA in (45.0, 46.0)",
      mpf("45") < KAPPA_SA < mpf("46"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C10]  V_inf = E_self / mu_0^2  (P187 value)")
print(SEP2)

V_inf_check = E_SELF / ALPHA_INV**2
print(f"  V_inf      = {nstr(V_inf_check, 10)}")
print(f"  V_inf/pi^2 = {nstr(V_inf_check / pi**2, 6)}")
check("C10 V_inf = E_self/mu_0^2",
      fabs(V_inf_check - V_INF) < mpf("1e-50"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C11]  UNIT NORMALISATION: V_inf_SA from integral = 1")
print("  V_inf_SA = 1 / integral_0^1 (1 - exp(-r/kappa_SA)) rho(r) dr")
print(SEP2)

I_norm = quad(lambda r: (1 - exp(-r / KAPPA_SA)) * rho_toe(r), [0, 1])
V_INF_SA = mpf(1) / I_norm
product  = V_INF_SA * I_norm

print(f"  Normalisation integral I = {nstr(I_norm, 10)}")
print(f"  V_inf_SA = 1/I           = {nstr(V_INF_SA, 10)}")
print(f"  Verification: V_inf_SA*I = {nstr(product, 12)}  (must be 1)")
check("C11 V_inf_SA * I = 1 to 1e-10",
      fabs(product - 1) < mpf("1e-10"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C12]  EIGENVALUE SHIFT at kappa_SA")
print("  delta_E_n = integral_0^1 V_self(r)*2*sin^2(n*pi*r) dr")
print(SEP2)

delta_E1 = quad(
    lambda r: V_yukawa(r, KAPPA_SA) * 2 * sin(pi * r)**2, [0, 1])
delta_E2 = quad(
    lambda r: V_yukawa(r, KAPPA_SA) * 2 * sin(2 * pi * r)**2, [0, 1])
analytic_approx = V_INF / (2 * KAPPA_SA)

print(f"  delta_E1 (n=1) = {nstr(delta_E1, 8)}")
print(f"  delta_E2 (n=2) = {nstr(delta_E2, 8)}")
print(f"  Approx V_inf/(2*kappa_SA) = {nstr(analytic_approx, 8)}")
print(f"  S^3 eigenvalue gap (E2-E1=5)")
print(f"  delta_E1 / gap = {nstr(delta_E1 / 5, 6)}")

check("C12a delta_E1 < 1e-4  (well below S^3 gap 5)",
      delta_E1 < mpf("1e-4"))
check("C12b delta_E1 within factor 2 of analytic approx",
      fabs(delta_E1 - analytic_approx) < analytic_approx)


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C13]  DIFFERENTIAL SHIFT: Yukawa vs beta=2 at kappa_SA")
print(SEP2)

def V_b2_kSA(r):
    return V_INF * (1 - exp(-(r / KAPPA_SA)**2))

dE1_b2 = quad(
    lambda r: V_b2_kSA(r) * 2 * sin(pi * r)**2, [0, 1])
dE2_b2 = quad(
    lambda r: V_b2_kSA(r) * 2 * sin(2 * pi * r)**2, [0, 1])

diff_E1  = fabs(delta_E1 - dE1_b2)
diff_gap = fabs((delta_E2 - delta_E1) - (dE2_b2 - dE1_b2))

print(f"  |delta_E1(Yukawa) - delta_E1(beta=2)| = {nstr(diff_E1, 6)}")
print(f"  |gap difference (E2-E1)|               = {nstr(diff_gap, 6)}")
check("C13 differential shift at kappa_SA < 1e-5",
      diff_E1 < mpf("1e-5"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C14]  FRACTIONAL MASS RATIO CHANGE at kappa_SA")
print(SEP2)

delta_gap_y  = delta_E2 - delta_E1
delta_gap_b2 = dE2_b2 - dE1_b2
frac_change  = MU_RATIO * fabs(delta_gap_y - delta_gap_b2)

print(f"  |gap(Yukawa) - gap(beta=2)| = {nstr(fabs(delta_gap_y - delta_gap_b2), 6)}")
print(f"  Frac. mass ratio change     = mu1/mu0 * |gap diff| = {nstr(frac_change, 6)}")
check("C14 fractional mass ratio change < 1e-4",
      frac_change < mpf("1e-4"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C15]  ODE UNIQUENESS: f''=kappa^{-2}*f, f(0)=1, f->0 forces exp(-r/kappa)")
print(SEP2)

kp2 = KAPPA_P187

def f_yuk(r):    return exp(-r / kp2)
def f_yuk_pp(r): return exp(-r / kp2) / kp2**2

max_res_f = max(fabs(f_yuk_pp(r) - f_yuk(r) / kp2**2) for r in r_coarse)
print(f"  max|f'' - kappa^{{-2}}*f| = {nstr(max_res_f, 6)}")
print(f"  f(0) = {nstr(f_yuk(mpf(0)), 4)}  (should be 1)")
print(f"  f(10/kappa) = {nstr(f_yuk(10*kp2), 6)}  (decays; should be < 1e-4)")

check("C15a f'' = kappa^{-2}*f to 1e-50", max_res_f < mpf("1e-50"))
check("C15b f(0) = 1", fabs(f_yuk(mpf(0)) - 1) < mpf("1e-50"))
check("C15c f decays: f(10/kappa) < 1e-4", f_yuk(10*kp2) < mpf("1e-4"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C16]  TWO DISTINCT PHYSICAL SCALES")
print(SEP2)

ratio_scales = KAPPA_SA / KAPPA_P187
print(f"  kappa_SA / kappa_P187 = {nstr(ratio_scales, 8)}")
check("C16 kappa_SA / kappa_P187 > 1e4",
      ratio_scales > mpf("1e4"))


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP2}")
print("[C17]  P187 BOUND: V_inf * sqrt(kappa_P187) ~ 3.2e-5")
print(SEP2)

p187_bound = V_INF * sqrt(KAPPA_P187)
print(f"  V_inf * sqrt(kappa_P187) = {nstr(p187_bound, 6)}  (expected ~3.2e-5)")
check("C17a P187 bound in (1e-5, 1e-4)",
      mpf("1e-5") < p187_bound < mpf("1e-4"))

print(f"  delta_E1(kappa_SA) = {nstr(delta_E1, 6)}")
print(f"  P187 bound         = {nstr(p187_bound, 6)}")
check("C17b delta_E1(kappa_SA) < P187 bound",
      delta_E1 < p187_bound)


# ─────────────────────────────────────────────────────────────────────────────
print(f"\n{SEP}")
print("SUMMARY")
print(SEP)
print(f"""
  A232 — SR1 Closure via the Screening Axiom

  C01: TOE constants consistent with P03.
  C02: Yukawa ODE residual < 1e-50 at kappa_P187.
  C03: Yukawa ODE residual < 1e-50 at kappa_SA.
  C04: beta=2 ODE residual > 1e-3 (at r ~ kappa_P187).
  C05: Linear profile ODE residual > 1e-3.
  C06: Quadratic profile ODE residual > 1e-3.
  C07: Square-root profile ODE residual > 1e-3.
  C08: Sinusoidal profile ODE residual > 1e-3.
  C09: kappa_SA = alpha_inv/3 in (45, 46).
  C10: V_inf = E_self/mu_0^2.
  C11: Unit normalisation: V_inf_SA * I = 1 to 1e-10.
  C12: Eigenvalue shift at kappa_SA < 1e-4.
  C13: Differential shift Yukawa vs beta=2 at kappa_SA < 1e-5.
  C14: Fractional mass ratio change < 1e-4.
  C15: Yukawa eigenfunction ODE uniqueness (f''=kappa^{{-2}}f, BCs).
  C16: kappa_SA / kappa_P187 > 1e4 (distinct scales).
  C17: P187 bound V_inf*sqrt(kappa_P187) ~ 3.2e-5; delta_E1 < bound.

  SR1 status: CLOSED (modulo kappa derivation from P09/P18).
  SA narrows R4-admissible family: inf-dim -> Yukawa (1-param residual).
  kappa_SA = alpha_inv/3 ~ 45.7  from S^3 mass gap at ell=1.

  Key values:
    alpha_inv    = {nstr(ALPHA_INV, 8)}
    kappa_P187   = {nstr(KAPPA_P187, 6)}  (alpha^{{5/4}})
    kappa_SA     = {nstr(KAPPA_SA, 8)}  (alpha_inv/3)
    V_inf        = {nstr(V_INF, 6)}
    V_inf_SA     = {nstr(V_INF_SA, 6)}  (unit-normalised at kappa_SA)
    delta_E1     = {nstr(delta_E1, 6)}  (eigenvalue shift at kappa_SA)
    P187 bound   = {nstr(p187_bound, 6)}
""")

print(f"\n{'='*60}\nRESULT: {PASS_COUNT} PASS / {FAIL_COUNT} FAIL")

if FAIL_COUNT > 0:
    raise SystemExit(1)

sys.exit(0)