#!/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)