#!/usr/bin/env python3
"""
verify_P286.py — Verifier for Addendum 286 (the tau residual).

  S1  Exclusion                    — checks 1-4
  S2  One-parameter families       — checks 5-6
  S3  Offset solve; pinned targets — checks 7-9
"""
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
ES = 13.177
mu1 = 16*PI**3/5 + 3*PI**2/4 + 2*PI/3
beta = 6*mu1/OM

print("S1  Exclusion")
N = 2000; 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)
l = np.linalg.eigvalsh(np.diag(2.0/dx**2 + V) + np.diag(off, 1)
                       + np.diag(off, -1))[:3]
check(1, "settled spectrum (N=2000): %s" % np.round(l, 4).tolist(),
      abs(l[0] - 16.5069) < 0.002 and abs(l[2] - 102.0204) < 0.01)
me, mmu, mtau = 0.51099895, 105.6583755, 1776.93
Lmu, Ltau = math.log(mmu/me), math.log(mtau/me)
b2 = Lmu/math.log(l[1]/l[0]); b3 = Ltau/math.log(l[2]/l[0])
check(2, "required exponents: beta_2 = %.4f, beta_3 = %.4f vs derived "
         "%.4f" % (b2, b3, beta),
      abs(b2 - 4.7069) < 0.001 and abs(b3 - 4.4768) < 0.001)
check(3, "constant beta excluded: drops %.2f%% and %.2f%% (monotone, "
         "level-increasing)" % (100*(1-b2/beta), 100*(1-b3/beta)),
      (1-b2/beta) > 0.01 and (1-b3/beta) > 0.055)
R = (beta - b3)/(beta - b2)
check(4, "drop ratio (beta-beta_3)/(beta-beta_2) = %.3f — the number any "
         "correction must reproduce" % R, 5.2 < R < 5.5)

print("S2  One-parameter families (predicted drop ratios vs %.2f)" % R)
fams = {
    "level-linear (n-1)": 2.0,
    "level-quadratic (n-1)^2": 4.0,
    "log-eigenvalue": math.log(l[2]/l[0])/math.log(l[1]/l[0]),
    "sqrt-eigenvalue": (math.sqrt(l[2])-math.sqrt(l[0]))
                       / (math.sqrt(l[1])-math.sqrt(l[0])),
    "eigenvalue-linear": (l[2]-l[0])/(l[1]-l[0]),
    "index-exponential e^(n-2)": math.e,
}
worst = min(abs(v - R) for v in fams.values())
for k, v in fams.items():
    print(f"        {k:<28s} {v:.3f}")
check(5, "no family reaches %.2f (nearest, level-quadratic, misses by "
         "%.2f = %.0f%% — after fitting the muon, its tau correction is "
         "off by a third); dressing (1+mu1 a^2)^n two orders too small"
         % (R, worst, 100*worst/R), worst > 1.0)
check(6, "verdict: tau residual requires new structure (derived beta(n) "
         "or V_obs correction); OI-286-1 filed", True)

print("S3  Offset solve; pinned targets")
f = lambda L0: math.log((l[1]-L0)/(l[0]-L0)) \
               / math.log((l[2]-L0)/(l[0]-L0)) - Lmu/Ltau
a, b = -50.0, l[0] - 0.05
fa = f(a)
for _ in range(300):
    m = (a + b)/2
    if f(m)*fa > 0:
        a = m; fa = f(a)
    else:
        b = m
L0 = (a + b)/2
bp = Lmu/math.log((l[1]-L0)/(l[0]-L0))
check(7, "offset family fits exactly (2 params, 2 targets) at lambda_0 = "
         "%.4f, beta' = %.4f — pinned pair for future proposals"
         % (L0, bp), abs(L0 - 5.203) < 0.01 and abs(bp - 3.797) < 0.005)
check(8, "neither pinned value is a corpus constant at bar: "
         "6(mu1/mu0)^2 = %.4f (off %.2f%%), lambda_1/pi = %.4f (off %.2f%%) "
         "— flags only" % (6*(mu1/OM)**2, 100*abs(6*(mu1/OM)**2/bp - 1),
                           l[0]/PI, 100*abs((l[0]/PI)/L0 - 1)),
      abs(6*(mu1/OM)**2/bp - 1) > 0.004 and abs((l[0]/PI)/L0 - 1) > 0.008)
check(9, "honest scorecard recorded: electron anchor, muon 1.1%, tau 6.0% "
         "in exponent; failure mode characterized", True)

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