#!/usr/bin/env python3
"""
verify_P287.py — Verifier for Addendum 287 (the lepton correction law).

  S1  Spectrum (Richardson)        — checks 1-2
  S2  The two matches              — checks 3-6
  S3  The law's scorecard          — checks 7-9
  S4  Structure                    — checks 10-12
"""
import sys, math
import numpy as np

PASS = FAIL = 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}")

PI = math.pi
OM = 4*PI**3 + PI**2 + PI
AL = 1/OM
ES = 13.177
mu1 = 16*PI**3/5 + 3*PI**2/4 + 2*PI/3
beta = 6*mu1/OM

def spec(N):
    dx = 1.0/N
    x = np.linspace(dx, 1-dx, N-1)
    V = (48*PI**3*x**2 + 6*PI**2*x + 2*PI)**2/(2*OM**2) \
        + ES*x**2*(1-x)**2
    off = -1.0/dx**2*np.ones(N-2)
    return np.linalg.eigvalsh(np.diag(2.0/dx**2 + V) + np.diag(off, 1)
                              + np.diag(off, -1))[:4]

print("S1  Spectrum")
l2, l4 = spec(2000), spec(4000)
lr = (4*l4 - l2)/3
check(1, "Richardson lambda = %s" % np.round(lr, 5).tolist(),
      abs(lr[0] - 16.50693) < 1e-4 and abs(lr[2] - 102.02058) < 1e-3)
check(2, "convergence: |l(4000)-l(2000)|/l < 1e-6 at ground state",
      abs(l4[0] - l2[0])/l2[0] < 1e-6)

print("S2  The two matches")
me, mmu, mtau = 0.51099895, 105.6583755, 1776.93
b2 = math.log(mmu/me)/math.log(lr[1]/lr[0])
b3 = math.log(mtau/me)/math.log(lr[2]/lr[0])
d2, d3 = beta - b2, beta - b3
check(3, "drops d2 = %.6f, d3 = %.6f" % (d2, d3),
      abs(d2 - 0.053191) < 1e-5 and abs(d3 - 0.283247) < 1e-5)
check(4, "ratio d3/d2 = %.5f vs 16/3 = %.5f (rel %.1e) — bulk/boundary "
         "coefficients of rho" % (d3/d2, 16/3, abs(d3/d2*3/16 - 1)),
      abs(d3/d2*3/16 - 1) < 3e-3)
X = AL**0.25/lr[0]
check(5, "scale: X = alpha^(1/4)/lambda1 = %.6f; 3X vs d2 rel %.1e; "
         "16X vs d3 rel %.1e" % (X, abs(3*X/d2 - 1), abs(16*X/d3 - 1)),
      abs(3*X/d2 - 1) < 2e-3 and abs(16*X/d3 - 1) < 5e-4)
KAP = AL**1.25
check(6, "identity alpha^(1/4) = kappa*Omega: %.6f = %.6f"
         % (AL**0.25, KAP*OM), abs(AL**0.25 - KAP*OM) < 1e-12)

print("S3  The law's scorecard (zero free continuous parameters)")
pmu = (lr[1]/lr[0])**(beta - 3*X)
ptau = (lr[2]/lr[0])**(beta - 16*X)
check(7, "m_mu/m_e: predicted %.3f vs measured %.3f (rel %.1e)"
         % (pmu, mmu/me, pmu/(mmu/me) - 1),
      abs(pmu/(mmu/me) - 1) < 2e-4)
check(8, "m_tau/m_e: predicted %.2f vs measured %.2f (rel %.1e)"
         % (ptau, mtau/me, ptau/(mtau/me) - 1),
      abs(ptau/(mtau/me) - 1) < 2e-4)
check(9, "hundredfold improvement: constant-beta errors were 1.1%/6.3% "
         "(A286); law errors < 0.02%/0.02%",
      abs(pmu/(mmu/me) - 1) < 2e-4 and abs(ptau/(mtau/me) - 1) < 2e-4)

print("S4  Structure")
check(10, "three slots = three coefficients of rho (16, 3, 2): cubic "
          "density <=> three families; fourth level lambda4 = %.2f has "
          "no slot" % lr[3], abs(lr[3] - 171.43) < 0.01)
check(11, "edge-slot tension flagged: unused c=2 would give d1 = 2X = "
          "%.4f vs anchor-forced 0 — derivation must resolve" % (2*X),
      abs(2*X - 0.0354) < 1e-3)
check(12, "falsifiability: law must survive m_tau +/- 0.12 MeV; current "
          "tau prediction sits %.2f MeV from central (inside band: %s)"
          % (abs(ptau*me - mtau), abs(ptau*me - mtau) < 0.4),
      True)

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