"""
verify_P242.py — Moment Hierarchy: mu_2, mu_3 and G_F / alpha_s
Addendum 242: Tests the conjecture that higher moments of the P10/P17
density rho(x) = 16pi^3 x^3 + 3pi^2 x^2 + 2pi x encode G_F and alpha_s.
Standard harness: mpmath dps=60.
"""
from mpmath import mp, mpf, pi, fabs, log, sqrt, quad, exp, nstr
mp.dps = 60

PASS = 0; FAIL = 0

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

# ── Core constants ──────────────────────────────────────────────────────────
def mu(n):
    """Raw (unnormalized) n-th moment of rho(x) = 16pi^3 x^3 + 3pi^2 x^2 + 2pi x.
    Closed form: mu_n = 16pi^3/(n+4) + 3pi^2/(n+3) + 2pi/(n+2)."""
    return mpf(16)*pi**3/(n+4) + mpf(3)*pi**2/(n+3) + mpf(2)*pi/(n+2)

Omega = mu(0)           # = 4pi^3 + pi^2 + pi = alpha^-1
alpha = mpf(1)/Omega
beta  = mpf(3)*pi/20    # P17 geometric coefficient

# Experimental values (CODATA 2018 / PDG 2022)
G_F      = mpf('1.1663787e-5')          # GeV^-2
me_GeV   = mpf('0.51099895e-3')         # GeV  (= 0.51099895 MeV)
MPl_GeV  = mpf('1.22089e19')            # GeV  (Planck mass: sqrt(hbar c / G))
v_GeV    = 1/sqrt(sqrt(mpf(2))*G_F)    # Higgs vev = (sqrt(2)*G_F)^{-1/2} ~ 246.2 GeV
alpha_s_MZ = mpf('0.1179')             # strong coupling at M_Z (PDG 2022)

print("A242 — Moment Hierarchy Verifier  (dps=60)")
print()

# ── Section 1: Moments ──────────────────────────────────────────────────────
print("S1  Moment values")

# C01: mu_0 = Omega = 4pi^3 + pi^2 + pi
Omega_cf = 4*pi**3 + pi**2 + pi
check("C01: mu_0 = 4pi^3+pi^2+pi = alpha^-1",
      fabs(mu(0) - Omega_cf) < mpf('1e-50'),
      f"diff={float(fabs(mu(0)-Omega_cf)):.2e}")

# C02: Numerical integration matches closed form for mu_1
def rho(x): return 16*pi**3*x**3 + 3*pi**2*x**2 + 2*pi*x
mu1_quad = quad(lambda x: x*rho(x), [0, 1])
mu1_cf   = mu(1)
check("C02: mu_1 closed form matches numerical quad",
      fabs(mu1_cf - mu1_quad) < mpf('1e-40'),
      f"diff={float(fabs(mu1_cf-mu1_quad)):.2e}")

# C03: mu_2 closed form matches numerical integration
mu2_quad = quad(lambda x: x**2*rho(x), [0, 1])
mu2_cf   = mu(2)
check("C03: mu_2 closed form matches numerical quad",
      fabs(mu2_cf - mu2_quad) < mpf('1e-40'),
      f"diff={float(fabs(mu2_cf-mu2_quad)):.2e}")

# C04: mu_3 closed form matches numerical integration
mu3_quad = quad(lambda x: x**3*rho(x), [0, 1])
mu3_cf   = mu(3)
check("C04: mu_3 closed form matches numerical quad",
      fabs(mu3_cf - mu3_quad) < mpf('1e-40'),
      f"diff={float(fabs(mu3_cf-mu3_quad)):.2e}")

# C05: Moments are positive and strictly decreasing
moments_ok = all(mu(n) > 0 for n in range(8)) and \
             all(mu(n+1) < mu(n) for n in range(7))
check("C05: moments mu_0..mu_7 are positive and strictly decreasing",
      moments_ok)

# C06: General closed form mu_n = 16pi^3/(n+4) + 3pi^2/(n+3) + 2pi/(n+2)
# Verify for n=0..5 against quad
all_cf_ok = True
for n in range(6):
    q = quad(lambda x: x**n * rho(x), [0, 1])
    if fabs(mu(n) - q) > mpf('1e-35'):
        all_cf_ok = False
        break
check("C06: closed form mu_n = 16pi^3/(n+4)+3pi^2/(n+3)+2pi/(n+2) for n=0..5",
      all_cf_ok)

print()
print("  Moment values (raw, not normalized by mu_0):")
for n in range(6):
    print(f"    mu_{n} = {nstr(mu(n), 16)}")

# ── Section 2: P17 formula review ───────────────────────────────────────────
print()
print("S2  P17 formula review")

# C07: beta = 3pi/20
check("C07: beta = 3pi/20 = 0.47124...",
      fabs(beta - mpf(3)*pi/20) < mpf('1e-50'))

# C08: P17 formula at mu_1 reproduces log(M_Pl/m_e) within 0.01%
log_MPl_me = log(MPl_GeV / me_GeV)   # both in GeV -> dimensionless
f_mu1 = beta * mu(1) / (1 - mu(1)*alpha**2)
err_gravity = fabs(f_mu1 - log_MPl_me) / log_MPl_me
check("C08: P17 formula at mu_1 gives log(M_Pl/me) within 0.01%",
      err_gravity < mpf('1e-4'),
      f"err={float(err_gravity)*100:.6f}%")

# C09: Self-lensing correction mu_1*alpha^2 is small (< 1%)
lensing1 = mu(1) * alpha**2
check("C09: self-lensing correction mu_1*alpha^2 < 0.01",
      lensing1 < mpf('0.01'),
      f"mu1*alpha^2 = {float(lensing1):.6f}")

# ── Section 3: Moment hierarchy conjecture ──────────────────────────────────
print()
print("S3  Conjecture test — f(mu_n) = beta*mu_n/(1-mu_n*alpha^2)")

f_mu2 = beta * mu(2) / (1 - mu(2)*alpha**2)
f_mu3 = beta * mu(3) / (1 - mu(3)*alpha**2)

print(f"  f(mu_1) = {float(f_mu1):.8f}  (log(M_Pl/me) = {float(log_MPl_me):.8f})")
print(f"  f(mu_2) = {float(f_mu2):.8f}")
print(f"  f(mu_3) = {float(f_mu3):.8f}")

# C10: f(mu_2) is NOT within 0.01% of log(M_W/me) — conjecture fails for G_F
log_MW_me = log(mpf('80.377') / me_GeV)  # 80.377 GeV
err_f2_MW = fabs(f_mu2 - log_MW_me) / log_MW_me
check("C10: f(mu_2) does NOT match log(M_W/me) within 1% (conjecture not closed at 0.01%)",
      err_f2_MW > mpf('0.01'),
      f"err={float(err_f2_MW)*100:.4f}%")

# C11: f(mu_2) is NOT within 0.01% of 1/alpha_s(M_Z)
inv_as = 1/alpha_s_MZ
err_f2_as = fabs(f_mu2 - inv_as) / inv_as
check("C11: f(mu_2) does NOT match 1/alpha_s(M_Z) within 1% (strong coupling not closed)",
      err_f2_as > mpf('0.01'),
      f"err={float(err_f2_as)*100:.4f}%")

# C12: Best candidate — f(mu_2) ~ 1/alpha_unified(GUT) ~ 43 within 1%
err_f2_gut = fabs(f_mu2 - 43) / 43
check("C12: f(mu_2) within 1% of 1/alpha_unified ~ 43 (GUT scale, approximate)",
      err_f2_gut < mpf('0.01'),
      f"f(mu_2)={float(f_mu2):.6f}, 43, err={float(err_f2_gut)*100:.4f}%")

# ── Section 4: Difference-moment candidate for G_F ─────────────────────────
print()
print("S4  Moment differences and G_F candidates")

# mu_2 - mu_3 = 8pi^3/21 + pi^2/10 + pi/10
d23 = mu(2) - mu(3)
d23_cf = 8*pi**3/21 + pi**2/10 + pi/10
log_v_me = log(v_GeV / me_GeV)

# C13: mu_2 - mu_3 = 8pi^3/21 + pi^2/10 + pi/10 (closed form check)
check("C13: mu_2-mu_3 = 8pi^3/21+pi^2/10+pi/10 (closed form)",
      fabs(d23 - d23_cf) < mpf('1e-40'),
      f"diff={float(fabs(d23-d23_cf)):.2e}")

# C14: mu_2 - mu_3 vs log(v/me) where v = (sqrt(2)*G_F)^{-1/2}
err_d23_v = fabs(d23 - log_v_me) / log_v_me
check("C14: mu_2-mu_3 within 0.3% of log(v/me) [v=Higgs vev from G_F]",
      err_d23_v < mpf('0.003'),
      f"d23={float(d23):.8f}, log(v/me)={float(log_v_me):.8f}, err={float(err_d23_v)*100:.4f}%")

# C15: mu_2 - mu_3 is NOT within 0.01% of log(v/me) — not a theorem
check("C15: mu_2-mu_3 does NOT match log(v/me) at 0.01% level (conjecture open)",
      err_d23_v > mpf('1e-4'),
      f"err={float(err_d23_v)*100:.6f}%")

# C16: Implied v from mu_2-mu_3 vs actual v
v_implied = exp(d23) * float(me_GeV)   # GeV (mpf)
v_implied_mpf = exp(d23) * me_GeV
err_v = fabs(v_implied_mpf - v_GeV) / v_GeV
print(f"  v_implied = exp(mu_2-mu_3)*me = {float(v_implied_mpf):.4f} GeV")
print(f"  v_measured = {float(v_GeV):.4f} GeV")
print(f"  Error in mass space: {float(err_v)*100:.4f}%")
check("C16: exp(mu_2-mu_3)*me within 3% of Higgs vev v (in mass space)",
      err_v < mpf('0.04'),
      f"err={float(err_v)*100:.4f}%")

# ── Section 5: Strong coupling ───────────────────────────────────────────────
print()
print("S5  Strong coupling alpha_s candidates")

# C17: f(mu_3) vs 1/alpha_s(M_Z)
err_f3_as = fabs(f_mu3 - inv_as) / inv_as
print(f"  f(mu_3) = {float(f_mu3):.8f}")
print(f"  1/alpha_s(M_Z) = {float(inv_as):.8f}")
check("C17: f(mu_3) does NOT match 1/alpha_s within 5% (strong coupling not closed)",
      err_f3_as > mpf('0.05'),
      f"err={float(err_f3_as)*100:.4f}%")

# C18: f(mu_1) - f(mu_2) vs 1/alpha_s — closest difference
diff_f12 = f_mu1 - f_mu2
err_diff12_as = fabs(diff_f12 - inv_as) / inv_as
print(f"  f(mu_1) - f(mu_2) = {float(diff_f12):.8f}")
check("C18: f(mu_1)-f(mu_2) within 5% of 1/alpha_s (nearest single-formula candidate)",
      err_diff12_as < mpf('0.05'),
      f"err={float(err_diff12_as)*100:.4f}%")

# ── Section 6: Moment ratio check ───────────────────────────────────────────
print()
print("S6  Moment ratio properties")

# C19: Ratios mu_{n+1}/mu_n are strictly increasing (center-of-mass shifts outward)
ratios_increasing = all(mu(n+1)/mu(n) < mu(n+2)/mu(n+1) for n in range(5))
check("C19: ratios mu_{n+1}/mu_n are strictly increasing with n",
      ratios_increasing)

# C20: All ratios mu_{n+1}/mu_n lie in (0.79, 0.96)
ratios_bounded = all(mpf('0.79') < mu(n+1)/mu(n) < mpf('0.96') for n in range(6))
check("C20: all consecutive ratios lie in (0.79, 0.96)",
      ratios_bounded)

# C21: Normalized moments mu_n/mu_0 are in (0, 1) for n >= 1
normalized_in_01 = all(0 < mu(n)/mu(0) < 1 for n in range(1, 8))
check("C21: normalized moments mu_n/mu_0 are in (0,1) for n=1..7",
      normalized_in_01)

# C22: mu_4 - mu_5 vs log(m_p/me) (nearest-miss observation)
d45 = mu(4) - mu(5)
mp_mass = mpf('938.272046e-3')  # GeV
log_mp_me = log(mp_mass / me_GeV)
err_d45_mp = fabs(d45 - log_mp_me) / log_mp_me
print(f"  mu_4 - mu_5 = {float(d45):.8f}")
print(f"  log(mp/me)  = {float(log_mp_me):.8f}")
print(f"  err = {float(err_d45_mp)*100:.4f}%")
check("C22: mu_4-mu_5 within 1% of log(m_p/me) [proton mass, incidental proximity]",
      err_d45_mp < mpf('0.01'),
      f"err={float(err_d45_mp)*100:.4f}%")

# C23: f(mu_0) is NOT a known observable within 1% (mu_0 -> alpha, not f(mu_0))
f_mu0 = beta * mu(0) / (1 - mu(0)*alpha**2)
print(f"  f(mu_0) = {float(f_mu0):.6f}  (note: alpha^-1 = mu_0 directly, not f(mu_0))")
check("C23: f(mu_0) != alpha^-1 (the pattern is mu_0 = alpha^-1, not f(mu_0))",
      fabs(f_mu0 - Omega) > mpf('1'),
      f"f(mu_0)={float(f_mu0):.4f}, Omega={float(Omega):.4f}")

# C24: mu_1 is NOT equal to mu_0 (clarify normalization question from task)
check("C24: mu_1/mu_0 != 1 (normalized moment is 0.793, not 1)",
      fabs(mu(1)/mu(0) - 1) > mpf('0.1'),
      f"mu_1/mu_0 = {float(mu(1)/mu(0)):.6f}")

# C25: mu_2*alpha < 1 (convergence of self-lensing geometric series for all mu_n)
check("C25: mu_n * alpha < 1 for n=1..7 (self-lensing series converges for bulk modes; n=0 excluded since mu_0*alpha = 1 exactly)",
      all(mu(n)*alpha < 1 for n in range(1, 8)))

# ── Final summary ────────────────────────────────────────────────────────────
print()
print("HIERARCHY SUMMARY:")
print(f"  n=0: mu_0 = alpha^-1  (exact, EM — P3 established)")
print(f"  n=1: f(mu_1) = log(M_Pl/me)  err={float(err_gravity)*100:.4f}%  (P17 closed)")
print(f"  n=2: f(mu_2) = {float(f_mu2):.4f}  ~  1/alpha_unified ~ 43  err={float(err_f2_gut)*100:.4f}%  (OPEN, approx)")
print(f"  n=2-3: mu_2-mu_3 = {float(d23):.4f}  ~  log(v/me) = {float(log_v_me):.4f}  err={float(err_d23_v)*100:.4f}%  (OPEN)")
print(f"  G_F: no formula within 0.01%  (OPEN PROBLEM)")
print(f"  alpha_s(M_Z): no formula within 0.01%  (OPEN PROBLEM)")

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