#!/usr/bin/env python3
"""
verify_P047.py -- Addendum 47: rho(x) as sector trace generator.

This verifier checks the negative spectral-density result and the positive
harmonic-diagonal generating-function theorem in
47_Addendum_SpectralDensityRoute.tex.

Result expected: clean. The Peirce spectrum, spectral trace 9*mu0, Gauss
quadrature nodes, harmonic diagonals D_k, rational generating function F(z),
and open OP-B/OP-D status all reproduce.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path

import numpy as np



PASS = FAIL = 0
_N = 0

def record(label, ok, computed="", claimed="", detail=""):
    """Modern-format check line; behavior-preserving port of verify_common."""
    global PASS, FAIL, _N
    _N += 1
    ok = bool(ok)
    desc = label
    if ok:
        PASS += 1
    else:
        FAIL += 1
        if "Expected" in detail:
            i = detail.find("Expected")
            desc = f"{label} -- {detail[i:]}"
            detail = detail[:i].rstrip().rstrip(";")
    print(f"  [{'PASS' if ok else 'FAIL'}] {_N:>2}. {desc}")
    if computed != "" or claimed != "":
        print(f"        computed: {computed}")
        print(f"        claimed : {claimed}")
    if detail:
        print(f"        {detail}")
    return ok

def check(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 record(label, ok, computed, claimed, err_detail + (f"; {detail}" if detail else ""))

print("P047 -- Spectral Density Route")

ROOT = Path(__file__).resolve().parents[1]
TEX = (ROOT / "47_Addendum_SpectralDensityRoute.tex").read_text()

PI = math.pi
ALPHA = [PI, PI**2, 4.0 * PI**3]
N = [1, 2, 3]
MU0 = sum(ALPHA)
MU1 = 16.0 * PI**3 / 5.0 + 3.0 * PI**2 / 4.0 + 2.0 * PI / 3.0
MU = MU1 / MU0


def mu(k: int) -> float:
    return 16.0 * PI**3 / (k + 4.0) + 3.0 * PI**2 / (k + 3.0) + 2.0 * PI / (k + 2.0)


def d_k(k: int) -> list[float]:
    return [2.0 / (k + 2.0), 3.0 / (k + 3.0), 4.0 / (k + 4.0)]


def trace_xsec_d(k: int) -> float:
    return sum(a * d for a, d in zip(ALPHA, d_k(k)))


def F(z: float) -> float:
    return 2.0 * PI / (z + 2.0) + 3.0 * PI**2 / (z + 3.0) + 16.0 * PI**3 / (z + 4.0)


peirce_values = [
    ALPHA[0],
    ALPHA[1],
    ALPHA[2],
    (ALPHA[0] + ALPHA[1]) / 2.0,
    (ALPHA[0] + ALPHA[2]) / 2.0,
    (ALPHA[1] + ALPHA[2]) / 2.0,
]
peirce_weights = [1, 1, 1, 8, 8, 8]
spectral_trace = sum(w * x for w, x in zip(peirce_weights, peirce_values))

# Three-node Gauss rule nodes: roots of the monic degree-3 polynomial
# orthogonal to 1, x, x^2 for the weight rho(x) on [0, 1].
gram = np.array(
    [
        [mu(2), mu(1), mu(0)],
        [mu(3), mu(2), mu(1)],
        [mu(4), mu(3), mu(2)],
    ],
    dtype=float,
)
rhs = -np.array([mu(3), mu(4), mu(5)], dtype=float)
a2, a1, a0 = np.linalg.solve(gram, rhs)
gauss_nodes = sorted(float(x.real) for x in np.roots([1.0, a2, a1, a0]))

check("mu0", MU0, 137.036304, rel=2e-9)
check("mu1", MU1, 108.716684, rel=3e-9)
check("MU", MU, 0.793342, rel=3e-7)

check("Peirce eigenvalue alpha1", peirce_values[0], 3.14159, rel=1e-6)
check("Peirce eigenvalue alpha2", peirce_values[1], 9.86960, rel=5e-7)
check("Peirce eigenvalue alpha3", peirce_values[2], 124.025, rel=1e-6)
check("Peirce eigenvalue alpha12", peirce_values[3], 6.5056, rel=5e-6)
check("Peirce eigenvalue alpha13", peirce_values[4], 63.583, rel=6e-6)
check("Peirce eigenvalue alpha23", peirce_values[5], 66.947, rel=6e-6)
check("Peirce multiplicity total", sum(peirce_weights), 27, rel=0)
record(
    "Peirce spectrum support is disjoint from [0,1]",
    min(peirce_values) > 1.0,
    computed=f"min spectrum = {min(peirce_values):.6f}",
    claimed="spec(L_Xsec) subset [pi, 4pi^3]",
)
check("spectral trace", spectral_trace, 9.0 * MU0, rel=1e-12)
check("spectral mean", spectral_trace / 27.0, MU0 / 3.0, rel=1e-12)
record(
    "spectral trace is not the Jordan trace",
    abs(spectral_trace - MU0) > 1000.0,
    computed=f"Tr_spec={spectral_trace:.6f}, Tr_JO={MU0:.6f}",
    claimed="Tr_spec differs from mu0 by factor 9",
)

expected_nodes = [0.33521, 0.68662, 0.93548]
for idx, (node, expected) in enumerate(zip(gauss_nodes, expected_nodes), start=1):
    check(f"three-node Gauss node {idx}", node, expected, rel=1e-5)

normalized_peirce = sorted([ALPHA[0] / MU0, ALPHA[1] / MU0, ALPHA[2] / MU0])
check("normalized Peirce alpha1/mu0", normalized_peirce[0], 0.0229, rel=2e-3)
check("normalized Peirce alpha2/mu0", normalized_peirce[1], 0.0720, rel=1e-3)
check("normalized Peirce alpha3/mu0", normalized_peirce[2], 0.9051, rel=7e-5)
record(
    "Gauss nodes do not match normalized Peirce eigenvalues",
    max(abs(a - b) for a, b in zip(gauss_nodes, normalized_peirce)) > 0.03,
    computed=f"Gauss={gauss_nodes}, Peirce/mu0={normalized_peirce}",
    claimed="no spectral identification by mu0 normalization",
)
record(
    "Gauss nodes do not match Peirce centers",
    max(abs(a - b) for a, b in zip(gauss_nodes, [2.0 / 3.0, 3.0 / 4.0, 4.0 / 5.0])) > 0.05,
    computed=f"Gauss={gauss_nodes}, centers={[2.0/3.0, 3.0/4.0, 4.0/5.0]}",
    claimed="nodes do not match Peirce center eigenvalues",
)

expected_moments = [137.036304, 108.716684, 90.175963, 77.062929, 67.289581]
for k, expected in enumerate(expected_moments):
    check(f"mu_{k} density moment", mu(k), expected, rel=1e-8)
    check(f"Tr(Xsec o D_{k})", trace_xsec_d(k), mu(k), rel=1e-12)
    check(f"F({k})", F(float(k)), mu(k), rel=1e-12)

record("D0 is identity", d_k(0) == [1.0, 1.0, 1.0], computed=d_k(0), claimed="[1,1,1]")
record(
    "D1 is Peirce center K",
    all(abs(a - b) < 1e-15 for a, b in zip(d_k(1), [2.0 / 3.0, 3.0 / 4.0, 4.0 / 5.0])),
    computed=d_k(1),
    claimed="[2/3, 3/4, 4/5]",
)
check("F pole residue at -2", 2.0 * PI, 2.0 * PI, rel=1e-12)
check("F pole residue at -3", 3.0 * PI**2, 3.0 * PI**2, rel=1e-12)
check("F pole residue at -4", 16.0 * PI**3, 16.0 * PI**3, rel=1e-12)
record("F(z) decays to zero", abs(F(1.0e6)) < 1.0e-3, computed=F(1.0e6), claimed="F(z)->0")

g0_g1_gap = -math.log(MU) / MU
record(
    "G0-G1 Stieltjes ratio formula is internally consistent",
    abs(g0_g1_gap - (-math.log(F(1.0) / F(0.0)) / (F(1.0) / F(0.0)))) < 1e-12,
    computed=g0_g1_gap,
    claimed="-ln(MU)/MU = -ln(F(1)/F(0))/(F(1)/F(0))",
)
record(
    "OP-B and OP-D remain open",
    "OP-B open" in TEX and "OP-D open" in TEX,
    computed="open status found",
    claimed="Conjectures B and D are not closed",
)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)
