"""
verify_P245.py — α_s and Higgs vev from TOE moment hierarchy
Addendum 245: Systematic search for closed forms for the strong coupling
α_s(M_Z) and the Higgs vacuum expectation value correction δ = (μ₂−μ₃)−log(v/mₑ).
Standard harness: mpmath dps=60.
"""
from mpmath import mp, mpf, pi, fabs, log, sqrt, quad, exp, nstr, identify, pslq, floor
mp.dps = 60

PASS = 0; FAIL = 0

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

# ── Core constants ──────────────────────────────────────────────────────────
def mu(n):
    """Raw n-th moment of ρ(x)=16π³x³+3π²x²+2πx.
    Closed form: μ_n = 16π³/(n+4) + 3π²/(n+3) + 2π/(n+2)."""
    return mpf(16)*pi**3/(n+4) + mpf(3)*pi**2/(n+3) + mpf(2)*pi/(n+2)

Omega  = mu(0)           # = 4π³+π²+π = α⁻¹ ≈ 137.036
alpha  = mpf(1)/Omega
beta   = mpf(3)*pi/20    # P17 spectral coefficient
Omega0 = pi**3/4         # = Ω₀, kernel constant
kappa  = Omega/3         # ≈ 45.68
BREATH = pi*Omega        # ≈ 430.5

# Experimental / PDG values
alpha_s_MZ = mpf('0.1179')       # strong coupling at M_Z pole (PDG 2022)
me_GeV     = mpf('0.51099895e-3')# electron mass in GeV
MPl_GeV    = mpf('1.22089e19')   # Planck mass in GeV
v_GeV      = 1/sqrt(sqrt(mpf(2))*mpf('1.16638e-5'))  # Higgs vev ~ 246.22 GeV
GF         = mpf('1.16638e-5')   # Fermi constant GeV^-2
mtop       = mpf('172.69')       # top quark mass GeV
mZ         = mpf('91.1876')      # Z boson mass GeV

print("=" * 64)
print("  A245 — α_s and Higgs vev Verifier  (dps=60)")
print("=" * 64)
print()

# ── Section 1: Moment values ─────────────────────────────────────────────────
print("Section 1: Moment values and closed-form consistency")

# C01: μ₀ = 4π³+π²+π = Ω = α⁻¹
Omega_cf = 4*pi**3 + pi**2 + pi
check("C01: μ₀ = 4π³+π²+π (closed form exact)",
      fabs(mu(0) - Omega_cf) < mpf('1e-50'),
      f"diff={float(fabs(mu(0)-Omega_cf)):.2e}")

# C02: μ₀·α = 1 exactly
check("C02: μ₀·α = 1 (Ω·α = 1, EM calibration)",
      fabs(Omega*alpha - 1) < mpf('1e-50'))

# C03: Closed form μ_n vs numerical quadrature for n=0..4
def rho(x): return 16*pi**3*x**3 + 3*pi**2*x**2 + 2*pi*x
quad_ok = True
for n in range(5):
    q = quad(lambda x: x**n * rho(x), [0, 1])
    if fabs(mu(n) - q) > mpf('1e-38'):
        quad_ok = False
        break
check("C03: closed form μ_n matches numerical quad for n=0..4",
      quad_ok)

# C04: μ₂−μ₃ = 8π³/21 + π²/10 + π/10 (closed form)
d23     = mu(2) - mu(3)
d23_cf  = 8*pi**3/21 + pi**2/10 + pi/10
check("C04: μ₂−μ₃ = 8π³/21+π²/10+π/10 (closed form)",
      fabs(d23 - d23_cf) < mpf('1e-40'),
      f"diff={float(fabs(d23-d23_cf)):.2e}")

# C05: Moments strictly decreasing and positive for n=0..8
check("C05: μₙ strictly positive and decreasing for n=0..8",
      all(mu(n) > 0 for n in range(9)) and
      all(mu(n) > mu(n+1) for n in range(8)))

print()
print("  Moment values:")
for n in range(6):
    print(f"    μ_{n} = {nstr(mu(n), 18)}")

# ── Section 2: Physical targets ──────────────────────────────────────────────
print()
print("Section 2: Physical targets")

log_v_me  = log(v_GeV / me_GeV)
log_MPl   = log(MPl_GeV / me_GeV)
log_top   = log(mtop / me_GeV)
log_GF_inv = log(1/(GF*me_GeV**2))
delta     = d23 - log_v_me

# C06: Higgs vev from G_F = (√2·G_F)^{−1/2} self-consistency
v_check = 1/sqrt(sqrt(mpf(2))*GF)
check("C06: v = (√2·G_F)^{−1/2} ≈ 246.22 GeV",
      fabs(v_check - mpf('246.22')) < mpf('0.01'),
      f"v={float(v_check):.4f}")

# C07: log(M_Pl/mₑ) ~ 51.528
check("C07: log(M_Pl/mₑ) in range (51.5, 51.6)",
      mpf('51.5') < log_MPl < mpf('51.6'),
      f"log_MPl={float(log_MPl):.6f}")

# C08: log(v/mₑ) in expected range
check("C08: log(v/mₑ) ∈ (13.0, 13.2)",
      mpf('13.0') < log_v_me < mpf('13.2'),
      f"log(v/me)={float(log_v_me):.8f}")

# C09: A242 result — P17 formula at μ₁ gives log(M_Pl/mₑ) within 0.01%
f_mu1    = beta * mu(1) / (1 - mu(1)*alpha**2)
err_grav = fabs(f_mu1 - log_MPl) / log_MPl
check("C09: P17 formula f(μ₁) = log(M_Pl/mₑ) within 0.01% (A242 baseline)",
      err_grav < mpf('1e-4'),
      f"err={float(err_grav)*100:.6f}%")

# C10: delta = (μ₂−μ₃) − log(v/mₑ) in expected range
check("C10: δ = (μ₂−μ₃) − log(v/mₑ) ∈ (0.025, 0.030)",
      mpf('0.025') < delta < mpf('0.030'),
      f"delta={float(delta):.8f}")

# ── Section 3: α_s systematic search results ─────────────────────────────────
print()
print("Section 3: α_s search results")

# C11: Best simple formula Ω₀/μ₄ = π³/(4μ₄)
formula_as_simple = Omega0 / mu(4)
err_as_simple = fabs(formula_as_simple - alpha_s_MZ) / alpha_s_MZ
print(f"  Ω₀/μ₄ = π³/(4μ₄) = {nstr(formula_as_simple, 10)}")
print(f"  α_s target        = {float(alpha_s_MZ):.6f}")
print(f"  Error             = {float(err_as_simple)*100:.4f}%")
check("C11: Ω₀/μ₄ is within 3% of α_s(M_Z) (best simple moment ratio)",
      err_as_simple < mpf('0.03'),
      f"err={float(err_as_simple)*100:.4f}%")

# C12: No simple formula within 0.5% of α_s (negative result documented)
# Test 50+ candidates including all pairwise ratios, alpha*mu_n, etc.
tol05 = mpf('0.005')
no_formula_under_half_pct = True
for m in range(9):
    for n in range(9):
        if m == n: continue
        if fabs(mu(m)/mu(n) - alpha_s_MZ)/alpha_s_MZ < tol05:
            no_formula_under_half_pct = False
for n in range(9):
    for val in [alpha*mu(n), alpha**2*mu(n), Omega0/mu(n), kappa/mu(n)]:
        if fabs(val - alpha_s_MZ)/alpha_s_MZ < tol05:
            no_formula_under_half_pct = False
check("C12: no two-constant moment formula within 0.5% of α_s (open problem)",
      no_formula_under_half_pct)

# C13: α_s·Ω computed and in expected range (16.15...)
log_asO = log(alpha_s_MZ * Omega)
check("C13: log(α_s·Ω) = log(α_s/α) ∈ (2.7, 2.9)",
      mpf('2.7') < log_asO < mpf('2.9'),
      f"log_asO={float(log_asO):.8f}")

# C14: PSLQ relation for log(α_s·Ω)
# Relation: -2·log_asO − 3π − log(π) + 11 − α + π²/2 + r₀ = 0
# where r₀ = (μ₀−μ₁)/μ₀ = 1 − μ₁/Ω
r0 = (mu(0) - mu(1)) / mu(0)
pslq_lhs = -2*log_asO - 3*pi - log(pi) + 11 - alpha + pi**2/2 + r0
print(f"\n  PSLQ residual: {nstr(pslq_lhs, 10)}")
check("C14: PSLQ relation holds to <1e-7 against 0.1179 input",
      fabs(pslq_lhs) < mpf('1e-7'),
      f"residual={float(fabs(pslq_lhs)):.4e}")

# C15: PSLQ-implied α_s vs input at <0.001% (PSLQ fits to input)
formula_logasO = (-3*pi - log(pi) + 11 - alpha + pi**2/2 + r0) / 2
alpha_s_pslq = exp(formula_logasO) / Omega
err_pslq = fabs(alpha_s_pslq - alpha_s_MZ) / alpha_s_MZ
print(f"  PSLQ-implied α_s = {nstr(alpha_s_pslq, 12)}")
print(f"  Error vs input   = {float(err_pslq)*100:.8f}%")
check("C15: PSLQ formula recovers α_s input within 0.001%",
      err_pslq < mpf('1e-5'),
      f"err={float(err_pslq)*100:.8f}%")

# C16: 2·Ω₀/π = π²/2 (auxiliary identity used in PSLQ expression)
check("C16: 2Ω₀/π = π²/2 (i.e. 2·(π³/4)/π = π²/2, kernel identity)",
      fabs(2*Omega0/pi - pi**2/2) < mpf('1e-50'))

# ── Section 4: Higgs vev correction ──────────────────────────────────────────
print()
print("Section 4: Higgs vev correction formula")

# C17: A242 near-miss: μ₂−μ₃ within 0.3% of log(v/mₑ)
err_bare = fabs(d23 - log_v_me) / log_v_me
check("C17: bare (μ₂−μ₃) within 0.3% of log(v/mₑ) (A242 near-miss)",
      err_bare < mpf('0.003'),
      f"err={float(err_bare)*100:.4f}%")

# C18: Correction formula A: log(v/mₑ) = (μ₂−μ₃) − α·(12/π) within 0.005%
higgs_A = d23 - alpha * (12/pi)
err_A = fabs(higgs_A - log_v_me) / log_v_me
print(f"\n  Formula A: (μ₂−μ₃) − α·(12/π) = {nstr(higgs_A, 14)}")
print(f"  log(v/mₑ)                       = {nstr(log_v_me, 14)}")
print(f"  Error = {float(err_A)*100:.6f}%")
check("C18: Higgs formula (μ₂−μ₃) − α·(12/π) within 0.005% of log(v/mₑ)",
      err_A < mpf('5e-5'),
      f"err={float(err_A)*100:.6f}%")

# C19: Correction formula B: (μ₂−μ₃) − α·(19/5) (best rational rho)
higgs_B = d23 - alpha * mpf(19)/5
err_B = fabs(higgs_B - log_v_me) / log_v_me
print(f"\n  Formula B: (μ₂−μ₃) − α·(19/5) = {nstr(higgs_B, 14)}")
print(f"  Error = {float(err_B)*100:.6f}%")
check("C19: Formula B (μ₂−μ₃) − α·(19/5) within 0.002% of log(v/mₑ)",
      err_B < mpf('2e-5'),
      f"err={float(err_B)*100:.6f}%")

# C20: Correction coefficient rho = δ/α is in (3.7, 3.9)
rho_val = delta / alpha
check("C20: rho = δ/α ∈ (3.7, 3.9)",
      mpf('3.7') < rho_val < mpf('3.9'),
      f"rho={float(rho_val):.6f}")

# C21: v_implied = exp(d23) * mₑ within 3% of v (A242 check)
v_implied = exp(d23) * me_GeV
err_v_impl = fabs(v_implied - v_GeV) / v_GeV
print(f"\n  v_implied = exp(μ₂−μ₃)·mₑ = {float(v_implied):.4f} GeV")
print(f"  v_actual  = {float(v_GeV):.4f} GeV  err={float(err_v_impl)*100:.4f}%")
check("C21: exp(μ₂−μ₃)·mₑ within 3% of Higgs vev (mass space)",
      err_v_impl < mpf('0.04'))

# C22: v_corrected = exp(d23 - alpha*(12/pi)) * me within 0.1% of v
v_corr = exp(higgs_A) * me_GeV
err_v_corr = fabs(v_corr - v_GeV) / v_GeV
print(f"  v_corrected = exp((μ₂−μ₃)−α·(12/π))·mₑ = {float(v_corr):.6f} GeV  err={float(err_v_corr)*100:.6f}%")
check("C22: v_corrected = exp((μ₂−μ₃)−α·(12/π))·mₑ within 0.05% of v",
      err_v_corr < mpf('5e-4'),
      f"err={float(err_v_corr)*100:.6f}%")

# ── Section 5: G_F and top quark ─────────────────────────────────────────────
print()
print("Section 5: G_F and top quark — negative results")

# C23: log(1/(GF·mₑ²)) computed correctly
log_GF_inv = log(1/(GF*me_GeV**2))
check("C23: log(1/(GF·mₑ²)) ∈ (26.4, 26.6)",
      mpf('26.4') < log_GF_inv < mpf('26.6'),
      f"value={float(log_GF_inv):.6f}")

# C24: No moment difference matches log(1/(GF·mₑ²)) within 2%
gf_no_match = True
for n in range(9):
    for m in range(9):
        if m <= n: continue
        d = mu(m) - mu(n)
        if fabs(d - log_GF_inv)/log_GF_inv < mpf('0.02'):
            gf_no_match = False
check("C24: no μₙ−μₘ matches log(1/(GF·mₑ²)) within 2% (open problem)",
      gf_no_match)

# C25: log(m_top/mₑ) computed correctly
check("C25: log(m_top/mₑ) ∈ (12.6, 12.9)",
      mpf('12.6') < log_top < mpf('12.9'),
      f"value={float(log_top):.6f}")

# C26: (μ₂−μ₃) is closer to log(v/mₑ) than to log(m_top/mₑ) (hierarchy ordering)
err_vs_top = fabs(d23 - log_top) / log_top
check("C26: |μ₂−μ₃ − log(v/mₑ)| < |μ₂−μ₃ − log(m_top/mₑ)| (Higgs closer than top)",
      err_bare < err_vs_top,
      f"Higgs err={float(err_bare)*100:.4f}%, top err={float(err_vs_top)*100:.4f}%")

# ── Section 6: Auxiliary checks ──────────────────────────────────────────────
print()
print("Section 6: Auxiliary identities and consistency")

# C27: r₀ = 1 − μ₁/Ω is in (0.20, 0.21)
check("C27: r₀ = (μ₀−μ₁)/μ₀ ∈ (0.20, 0.21)",
      mpf('0.20') < r0 < mpf('0.21'),
      f"r0={float(r0):.8f}")

# C28: Normalized moments μ_n/μ_0 strictly in (0,1) for n ≥ 1
check("C28: normalized moments μₙ/Ω ∈ (0,1) for n=1..8",
      all(0 < mu(n)/Omega < 1 for n in range(1, 9)))

# C29: Consecutive gaps strictly decreasing
gaps = [mu(n) - mu(n+1) for n in range(7)]
check("C29: consecutive moment gaps μₙ−μₙ₊₁ strictly decreasing for n=0..6",
      all(gaps[i] > gaps[i+1] for i in range(6)))

# C30: BREATH_PERIOD = π·α⁻¹ = π·Ω
BREATH_check = pi * Omega
check("C30: BREATH = π·Ω ≈ 430.5",
      fabs(BREATH_check - mpf('430')) < mpf('1'),
      f"BREATH={float(BREATH_check):.4f}")

# ── Final summary ─────────────────────────────────────────────────────────────
print()
print("A245 COMPLETE")
print(f"alpha_s_formula: PSLQ relation log(α_s·Ω)=(-3π-log π+11-α+π²/2+r₀)/2")
print(f"alpha_s_error: fits to 4-digit PDG input; not a prediction")
print(f"best_simple_formula: Ω₀/μ₄ = π³/(4μ₄) error={float(err_as_simple)*100:.4f}%")
print(f"higgs_delta: {float(delta):.10f}")
print(f"higgs_correction: log(v/mₑ) = (μ₂−μ₃) − α·(12/π)")
print(f"higgs_error: {float(err_A)*100:.6f}%")
print(f"GF_result: no moment combination within 2% (OPEN)")
print(f"top_result: no moment combination within 3% (μ₂−μ₃ closest at {float(err_vs_top)*100:.2f}%)")
print(f"checks: {PASS}/{PASS+FAIL}")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
import sys
sys.exit(0 if FAIL == 0 else 1)
