#!/usr/bin/env python3
"""
verify_P002.py -- Paper 02: Geometric Spectral Theory.

This verifier checks the main numerical chain in
toe/02_Paper_GeometricSpectralTheory.tex: cubic-density moments, self-lensing
energy, dynamic beta comparison, QED beta ratio, the critical point, and the
mass-hierarchy conjecture.

Most basic moment and energy arithmetic is reproducible. The audit flags:
  * the QED beta ratio 70,205.483 is obtained only if the experimental alpha is
    used, while the displayed formula uses the geometric m0 value;
  * the comparator-density E_self values in the optimal-density table are not
    the values produced by the stated normalization;
  * the wave equation cannot have rho_cubic as a Neumann equilibrium as stated;
  * the critical-point proof has the wrong s -> 0 limit;
  * the zeta-prime numeric value is right, but the displayed log formula is not;
  * M_Pl / exp(C) with C ~= 1362 is not 2.7e16 GeV.
"""

import math
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))

PASS = FAIL = 0
_N = 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}")
    return ok


class Verifier:
    """Output adapter: identical check semantics, modern [PASS]/[FAIL] format."""

    def __init__(self, name):
        print(name)

    def check(self, label, computed, claimed, *, rel=1e-3, abs_tol=None, detail=""):
        if abs_tol is not None:
            ok = abs(computed - claimed) <= abs_tol
            err_detail = f"abs err={abs(computed - claimed):.6g}, tol={abs_tol:.6g}"
        else:
            if claimed == 0:
                ok = abs(computed) <= (rel or 1e-12)
                err_detail = f"abs value={abs(computed):.6g}, tol={rel:.6g}"
            else:
                err = (computed - claimed) / abs(claimed)
                ok = abs(err) <= (rel or 0)
                err_detail = f"rel err={100 * err:+.6g}%, tol={100 * (rel or 0):.6g}%"
        return self.record(label, ok, computed, claimed, detail, err_detail)

    def record(self, label, ok, computed="", claimed="", detail="", err_detail=""):
        global _N
        _N += 1
        desc = label + (f" -- {detail}" if detail else "")
        check(_N, desc, ok)
        if computed != "" or claimed != "":
            print(f"        computed: {computed}")
            print(f"        claimed : {claimed}")
        if err_detail:
            print(f"        {err_detail}")
        return ok

    def summary(self):
        print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
        return 1 if FAIL else 0


v = Verifier("P002 -- Geometric Spectral Theory")

PI = math.pi
GAMMA = 0.5772156649015329
PHI_GOLDEN = (1 + math.sqrt(5)) / 2

A3 = 16 * PI**3
A2 = 3 * PI**2
A1 = 2 * PI

M0 = 4 * PI**3 + PI**2 + PI
M1 = A3 / 5 + A2 / 4 + A1 / 3
ALPHA_GEOM = 1 / M0
ALPHA_EXP_INV = 137.035999084
ALPHA_EXP = 1 / ALPHA_EXP_INV


def rho_prime(x: float) -> float:
    return 48 * PI**3 * x**2 + 6 * PI**2 * x + 2 * PI


def rho_second(x: float) -> float:
    return 96 * PI**3 * x + 6 * PI**2


def phi_laplace(s: float) -> float:
    e = math.exp(-s)
    return (
        A3 * (6 - e * (6 + 6 * s + 3 * s**2 + s**3)) / s**4
        + A2 * (2 - e * (2 + 2 * s + s**2)) / s**3
        + A1 * (1 - e * (1 + s)) / s**2
    )


def solve_critical() -> float:
    target = M0**2 / M1
    lo, hi = 0.52, 0.53
    for _ in range(100):
        mid = (lo + hi) / 2
        if phi_laplace(mid) / mid > target:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


energy = 0.5 * (
    (48 * PI**3) ** 2 / 5
    + (6 * PI**2) ** 2 / 3
    + (2 * PI) ** 2
    + 2 * (48 * PI**3) * (6 * PI**2) / 4
    + 2 * (48 * PI**3) * (2 * PI) / 3
    + 2 * (6 * PI**2) * (2 * PI) / 2
)
e_self = energy / M0**2
delta_e = e_self - 4 * PI
beta_geom = M1 / M0
beta_dyn = 2 * PI / (3 * PHI_GOLDEN**2) - GAMMA / 100
delta_beta_over_beta = abs(beta_geom - beta_dyn) / beta_dyn
kappa_beta = 2 * delta_beta_over_beta
kappa_alpha = ALPHA_GEOM ** (5 / 4)

v.check("m0 = 4*pi^3 + pi^2 + pi", M0, 137.036303776, rel=1e-12)
v.check(
    "geometric alpha inverse vs CODATA relative percent",
    100 * (M0 - ALPHA_EXP_INV) / ALPHA_EXP_INV,
    0.000222,
    rel=2e-3,
)
v.check("m1 moment", M1, 108.716683780, rel=1e-11)
v.check("beta_geom = m1/m0", beta_geom, 0.793342207849, rel=2e-10)

v.check("Dirichlet energy E[rho]", energy, 247444.81, rel=5e-8)
v.check("self-lensing energy E_self", e_self, 13.17671, rel=5e-7)
v.check("floor 4*pi", 4 * PI, 12.56637, rel=5e-7)
v.check("gap Delta E", delta_e, 0.61034, rel=1e-4)
v.check("fractional gap Delta E / E_self", delta_e / e_self, 0.04633, rel=3e-4)

v.check("beta_dyn = 2*pi/(3 phi^2) - gamma/100", beta_dyn, 0.794215778, rel=3e-7)
v.check("static-dynamic relative difference", delta_beta_over_beta, 0.0011, rel=5e-4)
v.check("kappa from beta difference", kappa_beta, 0.0022, rel=5e-4)
v.check("kappa = alpha^(5/4) comparison ratio", kappa_beta / kappa_alpha, 1.031, rel=1e-3)

v.check(
    "Delta E / E_self agreement with 22*kappa is within 1.3 percent",
    100 * abs(delta_e / e_self - 22 * kappa_beta) / (22 * kappa_beta),
    1.3,
    rel=1e-2,
    detail="Expected fail: the stated quantities differ by about 4.27%, as the proof later acknowledges.",
)
v.check(
    "Delta E / E_self agreement with 22*kappa is within 5 percent",
    abs(delta_e / e_self - 22 * kappa_beta) / (22 * kappa_beta),
    0.0,
    abs_tol=0.05,
)

# The optimal-density table uses normalized densities with integral m0.
linear_e_self = 2.0
quadratic_e_self = 6.0
v.check("optimal-density table: uniform E_self", 0.0, 0.0, abs_tol=1e-12)
v.check(
    "optimal-density table: linear E_self",
    linear_e_self,
    0.171,
    rel=1e-3,
    detail="Expected fail: for rho=2*m0*x, E/m0^2 = 2, not 0.171.",
)
v.check(
    "optimal-density table: quadratic E_self",
    quadratic_e_self,
    0.383,
    rel=1e-3,
    detail="Expected fail: for rho=3*m0*x^2, E/m0^2 = 6, not 0.383.",
)
v.record(
    "optimal-density proof exhausts polynomial densities",
    False,
    computed="four monomial test cases",
    claimed="minimum among polynomial densities",
    detail="Expected fail: checking uniform/linear/quadratic/cubic monomials is not an optimization proof over polynomial coefficients.",
)

beta_qed_geom = 2 * ALPHA_GEOM**2 / (3 * PI)
ratio_geom_alpha = beta_geom / beta_qed_geom
beta_qed_exp = 2 * ALPHA_EXP**2 / (3 * PI)
ratio_exp_alpha = beta_geom / beta_qed_exp
log_planck_e = math.log(1.220890e19 / (5.10999e-4))
c_emp_exp = ratio_exp_alpha / log_planck_e
c_formula = 10 * M0**3 / (M0**2 + M1)

v.check("QED beta using displayed geometric m0 formula", beta_qed_geom, 1.130024e-5, rel=5e-7)
v.check(
    "beta_geom/beta_QED using displayed geometric m0 formula",
    ratio_geom_alpha,
    70205.483,
    rel=1e-6,
    detail="Expected fail: 70,205.483 uses experimental alpha, not beta_QED=2/(3*pi*m0^2).",
)
v.check("beta_geom/beta_QED using experimental alpha", ratio_exp_alpha, 70205.483, rel=1e-6)
v.check("log(M_Pl/m_e) from stated masses", log_planck_e, 51.527840, rel=1e-8)
v.check("empirical C from experimental-alpha ratio", c_emp_exp, 1362.477, rel=5e-7)
v.check(
    "exact C formula value",
    c_formula,
    1362.482800,
    abs_tol=0.001,
    detail="Expected fail: exact m0,m1 give about 1362.47528; the quoted value is a small arithmetic mismatch.",
)

s_star = solve_critical()
v.check("critical point s*", s_star, 0.525159133415, rel=1e-11)
v.check("critical point proximity to pi/6 percent", 100 * (s_star - PI / 6) / s_star, 0.297, rel=2e-3)
v.record(
    "critical-point proof limit as s -> 0+",
    False,
    computed="Phi(s)/s -> +infinity",
    claimed="Phi(s)/s -> Phi'(0) = -m1",
    detail="Expected fail: the root value is correct, but this proof step is not.",
)

rho_prime_0 = rho_prime(0.0)
rho_prime_1 = rho_prime(1.0)
v.check(
    "rho_cubic satisfies Neumann boundary at x=0",
    rho_prime_0,
    0.0,
    abs_tol=1e-12,
    detail="Expected fail: rho'(0)=2*pi.",
)
v.check(
    "rho_cubic satisfies Neumann boundary at x=1",
    rho_prime_1,
    0.0,
    abs_tol=1e-12,
    detail="Expected fail: rho'(1) is large and nonzero.",
)
v.check(
    "rho_cubic is static equilibrium of stated wave equation at x=1/2",
    rho_second(0.5),
    0.0,
    abs_tol=1e-12,
    detail="Expected fail: at rho=rho_cubic, the RHS leaves rho_cubic'' instead of zero.",
)
v.record(
    "linearization reduces to eta_tt = (1-kappa) eta_xx",
    False,
    computed="eta_xx - kappa*eta*rho_cubic'' plus a nonzero zeroth-order residual",
    claimed="(1-kappa) eta_xx",
    detail="Expected fail: differentiating -kappa(rho-rho_cubic)rho'' gives a multiplicative rho_cubic'' term, not -kappa eta''.",
)

rhs_log = (3 * PI / 20) * M1 / (1 - M1 * ALPHA_GEOM**2)
v.check("geometric-series mass-log RHS", rhs_log, 51.529851, rel=5e-9)
v.check("geometric-series relative error percent", 100 * abs(rhs_log - log_planck_e) / log_planck_e, 0.0039, rel=2e-3)
v.check("series convergence parameter m1*alpha^2", M1 * ALPHA_GEOM**2, 0.005789285, rel=2e-7)
v.check("leading mass-log term", (3 * PI / 20) * M1, 51.231530, rel=5e-8)
v.check("alpha^2 correction term", (3 * PI / 20) * M1**2 * ALPHA_GEOM**2, 0.2966, rel=1e-4)

zeta_prime_integral = A3 / 16 + A2 / 9 + A1 / 4
zeta_prime_printed_log_formula = A3 * math.log(1 / 4) + A2 * math.log(1 / 3) + A1 * math.log(1 / 2)
v.check("zeta'(0) from integral definition", zeta_prime_integral, 35.867, rel=2e-6)
v.check(
    "displayed zeta'(0) log formula",
    zeta_prime_printed_log_formula,
    35.867,
    rel=1e-3,
    detail="Expected fail: the printed sum A_k log(1/k) is negative and not the derivative of the integral zeta function.",
)
v.check("zeta'(0)/zeta(0)", zeta_prime_integral / M0, 0.262, rel=2e-3)
v.check("one-loop action -1/2 zeta'(0)", -0.5 * zeta_prime_integral, -17.93, rel=2e-4)

log10_e_geom = math.log10(1.220890e19) - 1362.477 / math.log(10)
v.check(
    "log10(M_Pl / exp(C)) with C=1362.477",
    log10_e_geom,
    math.log10(2.7e16),
    abs_tol=0.1,
    detail="Expected fail: exp(1362) is enormous; the scale is around 10^-573 GeV, not 2.7e16 GeV.",
)

sys.exit(v.summary())
