"""
verify_P239.py  — Addendum 239: Gravity / Cosmology Map
Numerical verifier using mpmath (dps=60).

Checks:
  C01  Fine-structure constant α = 1/(4π³+π²+π)
  C02  Dark-energy corner residual Λ₀ = 1 − π²/32 vs Planck 2018
  C03  Cosmic breathing period T_breath = π·α⁻¹ ≈ 430.51 (within 2 of 432)
  C03b T_breath exact formula check |T − (4π⁴+π³+π²)| < 1e-20
  C04  Gravity coupling β = 3π/20 value
  C05a First moment μ₁ = 16π³/5 + 3π²/4 + 2π/3 ≈ 108.7167
  C05b μ₁ via direct integration agrees with closed form to 1e-20
  C06  P17 prediction log(M_Pl/m_e) ≈ 51.5299 vs observed 51.5271; |Δ| < 0.01
  C07  Observed log(M_Pl/m_e) from SI masses ≈ 51.527; verify in [51.5, 51.6]
  C08  γ = 3/4 exact
  C09  Λ₀ vs updated Planck 2020 Ω_Λ = 0.6889 ± 0.0056; within 1σ
  C10  ε = π/α⁻¹ (Hopf edge fraction) ≈ 0.0229
  C11  Moment hierarchy: μ₀ > μ₁ > μ₂ > μ₃ (monotone decay)
  C12  μ₁·α² ≈ 5.79e-3 (self-lensing parameter small)
  C13  Self-lensing denominator 1 − μ₁·α² ∈ (0.99, 1.00)
  C14  β·μ₁ ≈ 51.22 (leading-order gravity approx)
  C15  V_4(B⁴) = π²/2; corner residual formula cross-check
"""

from mpmath import mp, mpf, pi, log, quad, fabs, sqrt
mp.dps = 60

PASS = 0
FAIL = 0

_N = 0

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

print("=" * 60)
print("verify_P239.py  —  Addendum 239: Gravity/Cosmology Map")
print(f"mpmath dps = {mp.dps}")
print("=" * 60)

# ── Constants ──────────────────────────────────────────────────────────────

# Fine-structure constant
mu0 = 4*pi**3 + pi**2 + pi          # = α⁻¹
alpha = 1 / mu0

# First moment (closed form from P17)
mu1_formula = mpf('16')*pi**3/5 + mpf('3')*pi**2/4 + 2*pi/3

# Second moment
mu2_formula = mpf('16')*pi**3/6 + mpf('3')*pi**2/5 + 2*pi/4

# Third moment
mu3_formula = mpf('16')*pi**3/7 + mpf('3')*pi**2/6 + 2*pi/5

# Gravity coupling
beta = 3*pi/20

# Dark-energy corner residual
Lambda0 = 1 - pi**2/32

# Breathing period
T_breath = pi * mu0

# γ (boundary/bulk ratio)
gamma = mpf('3') / 4

# Hopf edge fraction
epsilon = pi / mu0

# ── C01: Fine-structure constant ───────────────────────────────────────────
print("\n── C01  Fine-structure constant α ──────────────────────────")
alpha_codata = mpf('1') / mpf('137.035999084')
check("C01  |α_TOE − α_CODATA| < 1e-4",
      fabs(alpha - alpha_codata) < mpf('1e-4'),
      f"|Δα| = {fabs(alpha - alpha_codata)}")

# ── C02: Λ₀ vs Planck 2018 ────────────────────────────────────────────────
print("\n── C02  Λ₀ = 1 − π²/32 vs Planck 2018 ─────────────────────")
Omega_L_Planck2018 = mpf('0.6847')
sigma_2018 = mpf('0.0073')
check("C02  |Λ₀ − Ω_Λ(2018)| < 1σ",
      fabs(Lambda0 - Omega_L_Planck2018) < sigma_2018,
      f"|Δ| = {float(fabs(Lambda0 - Omega_L_Planck2018)):.6f}, σ = {float(sigma_2018)}")

# ── C03: Breathing period vs 432 ──────────────────────────────────────────
print("\n── C03  Cosmic breathing period T_breath ────────────────────")
check("C03  |T_breath − 432| < 2  (known 0.35% gap, within tolerance)",
      fabs(T_breath - 432) < 2,
      f"T_breath = {float(T_breath):.4f}")

check("C03b T_breath = 4π⁴+π³+π² (exact formula self-consistency)",
      fabs(T_breath - (4*pi**4 + pi**3 + pi**2)) < mpf('1e-40'),
      f"residual = {fabs(T_breath - (4*pi**4 + pi**3 + pi**2))}")

# ── C04: Gravity coupling β ───────────────────────────────────────────────
print("\n── C04  Gravity coupling β = 3π/20 ─────────────────────────")
beta_ref = mpf('0.471239')
check("C04  |β − 0.471239| < 1e-4",
      fabs(beta - beta_ref) < mpf('1e-4'),
      f"β = {float(beta):.6f}")

# ── C05a: First moment μ₁ (closed form) ──────────────────────────────────
print("\n── C05a μ₁ closed form ─────────────────────────────────────")
check("C05a |μ₁ − 108.717| < 0.01",
      fabs(mu1_formula - mpf('108.717')) < mpf('0.01'),
      f"μ₁ = {float(mu1_formula):.6f}")

# ── C05b: μ₁ via numerical integration ──────────────────────────────────
print("\n── C05b μ₁ via direct integration ─────────────────────────")
def rho(x):
    return 16*pi**3*x**3 + 3*pi**2*x**2 + 2*pi*x

def x_rho(x):
    return x * rho(x)

mu1_integral = quad(x_rho, [0, 1])
check("C05b |μ₁_integral − μ₁_formula| < 1e-20",
      fabs(mu1_integral - mu1_formula) < mpf('1e-20'),
      f"residual = {fabs(mu1_integral - mu1_formula)}")

# ── C06: P17 gravity prediction vs observed ──────────────────────────────
print("\n── C06  P17 gravity formula prediction ─────────────────────")
# Prediction
lhs_predicted = beta * mu1_formula / (1 - mu1_formula * alpha**2)

# Observed: ln(M_Pl / m_e) using P17 values (GeV)
# M_Pl = 1.22e19 GeV, m_e = 5.11e-4 GeV
M_Pl_GeV = mpf('1.22e19')
m_e_GeV  = mpf('5.11e-4')
lhs_observed_P17 = log(M_Pl_GeV / m_e_GeV)

check("C06  |predicted − observed(P17)| < 0.01",
      fabs(lhs_predicted - lhs_observed_P17) < mpf('0.01'),
      f"predicted={float(lhs_predicted):.4f}, observed={float(lhs_observed_P17):.4f}")

# ── C07: Observed ln(M_Pl/m_e) from SI masses ────────────────────────────
print("\n── C07  Observed ln(M_Pl/m_e) from SI masses ───────────────")
# M_Pl (reduced Planck mass in kg) = 2.176434e-8 kg
# m_e = 9.1093837015e-31 kg
M_Pl_SI = mpf('2.176434e-8')
m_e_SI  = mpf('9.1093837015e-31')
lhs_SI  = log(M_Pl_SI / m_e_SI)

check("C07  ln(M_Pl/m_e)_SI ∈ [51.50, 51.55]",
      mpf('51.50') < lhs_SI < mpf('51.55'),
      f"ln(M_Pl/m_e)_SI = {float(lhs_SI):.4f}")

check("C07b |predicted − ln(M_Pl/m_e)_SI| < 0.01",
      fabs(lhs_predicted - lhs_SI) < mpf('0.01'),
      f"|Δ| = {float(fabs(lhs_predicted - lhs_SI)):.5f}")

# ── C08: γ = 3/4 exact ────────────────────────────────────────────────────
print("\n── C08  γ = 3/4 (boundary/bulk ratio) ──────────────────────")
check("C08  γ == 3/4 exactly",
      gamma == mpf('3') / mpf('4'))

# ── C09: Λ₀ vs Planck 2020 ────────────────────────────────────────────────
print("\n── C09  Λ₀ vs Planck 2020 Ω_Λ = 0.6889 ± 0.0056 ──────────")
Omega_L_2020 = mpf('0.6889')
sigma_2020   = mpf('0.0056')
discrepancy  = fabs(Lambda0 - Omega_L_2020)
n_sigma      = discrepancy / sigma_2020
check("C09  |Λ₀ − Ω_Λ(2020)| < 1σ",
      n_sigma < 1,
      f"discrepancy = {float(discrepancy):.4f}, n_sigma = {float(n_sigma):.2f}")

# ── C10: Edge fraction ε ──────────────────────────────────────────────────
print("\n── C10  Hopf edge fraction ε = π/α⁻¹ ───────────────────────")
check("C10  |ε − 0.0229| < 1e-4",
      fabs(epsilon - mpf('0.0229')) < mpf('1e-4'),
      f"ε = {float(epsilon):.5f}")

# ── C11: Moment hierarchy monotone ────────────────────────────────────────
print("\n── C11  Moment hierarchy μ₀ > μ₁ > μ₂ > μ₃ ────────────────")
check("C11  μ₀ > μ₁",  mu0      > mu1_formula)
check("C11  μ₁ > μ₂",  mu1_formula > mu2_formula)
check("C11  μ₂ > μ₃",  mu2_formula > mu3_formula)

# ── C12: Self-lensing parameter ───────────────────────────────────────────
print("\n── C12  Self-lensing parameter μ₁α² ────────────────────────")
self_lens = mu1_formula * alpha**2
check("C12  μ₁·α² < 0.01  (small perturbation)",
      self_lens < mpf('0.01'),
      f"μ₁α² = {float(self_lens):.6f}")
check("C12b |μ₁·α² − 5.79e-3| < 1e-4",
      fabs(self_lens - mpf('5.79e-3')) < mpf('1e-4'),
      f"μ₁α² = {float(self_lens):.5e}")

# ── C13: Self-lensing denominator ─────────────────────────────────────────
print("\n── C13  Self-lensing denominator 1 − μ₁α² ─────────────────")
denom = 1 - mu1_formula * alpha**2
check("C13  1 − μ₁α² ∈ (0.99, 1.00)",
      mpf('0.99') < denom < mpf('1.00'),
      f"denominator = {float(denom):.6f}")

# ── C14: Leading-order gravity approximation ──────────────────────────────
print("\n── C14  Leading-order gravity approx β·μ₁ ───────────────────")
leading = beta * mu1_formula
check("C14  |β·μ₁ − 51.22| < 0.1",
      fabs(leading - mpf('51.22')) < mpf('0.1'),
      f"β·μ₁ = {float(leading):.4f}")

# ── C15: Volume of B⁴ and corner residual cross-check ────────────────────
print("\n── C15  V_4(B⁴) = π²/2 and corner residual ─────────────────")
V4_ball = pi**2 / 2
V4_cube = mpf('16')   # = 2^4
fraction_covered = V4_ball / V4_cube    # = π²/32
Lambda0_check    = 1 - fraction_covered
check("C15  V₄(B⁴) = π²/2 ≈ 4.935",
      fabs(V4_ball - mpf('4.935')) < mpf('0.001'),
      f"V₄ = {float(V4_ball):.5f}")
check("C15b |Λ₀_formula − Λ₀_check| < 1e-40",
      fabs(Lambda0 - Lambda0_check) < mpf('1e-40'))

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