#!/usr/bin/env python3
"""
verify_P291.py — Verifier for Addendum 291 (the ergodic dressing).

  S1  Closed forms                 — checks 1-3
  S2  The ergodic theorem          — checks 4-7
  S3  Data selection + direction   — checks 8-11
"""
import sys, os
import mpmath as mp

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

mp.mp.dps = 40
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 = mp.pi
OM = 4*PI**3 + PI**2 + PI
AL = 1/OM
KAP = AL**mp.mpf("1.25")

def E_closed(n):
    return (2*PI**3 - 6*PI/n**2, PI**2/2 - mp.mpf(3)/(4*n**2), PI/2)

print("S1  Closed forms")
ok = True; worst = 0
for n in (1, 2, 5, 50):
    num = [mp.quad(lambda x, c=c, p=p: c*x**p*mp.sin(n*PI*x)**2, [0, 1])
           for c, p in ((16*PI**3, 3), (3*PI**2, 2), (2*PI, 1))]
    cl = E_closed(n)
    dev = max(abs(num[i] - cl[i]) for i in range(3))
    worst = max(worst, dev); ok = ok and dev < mp.mpf(10)**-30
check(1, "closed forms match quadrature at n = 1,2,5,50 (worst dev %s)"
         % mp.nstr(worst, 3), ok)
check(2, "edge integral = pi/2 for every n (n-independent: the edge "
         "hears every mode equally)",
      abs(E_closed(1)[2] - E_closed(50)[2]) == 0)
check(3, "total = Omega/2 - (6pi+3/4)/n^2: residual at n=3: %s"
         % mp.nstr(abs(sum(E_closed(3)) - (OM/2 - (6*PI + mp.mpf(3)/4)/9)), 3),
      abs(sum(E_closed(3)) - (OM/2 - (6*PI + mp.mpf(3)/4)/9)) < mp.mpf(10)**-35)

print("S2  The ergodic theorem")
f = (4*PI**3/OM, PI**2/OM, PI/OM)
lim = tuple(v/(OM/2) for v in (2*PI**3, PI**2/2, PI/2))
check(4, "limit shares = measure shares EXACTLY (algebraic identity)",
      all(abs(lim[i] - f[i]) < mp.mpf(10)**-35 for i in range(3)))
def shares(n):
    E = E_closed(n); T = sum(E)
    return tuple(e/T for e in E)
s1 = shares(1)
check(5, "fundamental n=1 shares: (%.5f, %.5f, %.5f) — edge 40%% high; "
         "equipartition NOT trivially true" % s1,
      abs(s1[2] - mp.mpf("0.03211")) < 1e-5)
s2, s5 = shares(2), shares(5)
check(6, "1/n^2 convergence: edge shares %.5f (n=2), %.5f (n=5) -> %.5f"
         % (s2[2], s5[2], f[2]),
      abs(s2[2] - f[2]) < abs(s1[2] - f[2])/3 and
      abs(s5[2] - f[2]) < abs(s2[2] - f[2])/5)
coef5 = (s5[2]/f[2] - 1)*25
check(7, "edge correction coefficient (asymptotic, measured at n=5): "
         "s_edge(n)/f_edge - 1 = %.3f/n^2 (positive for all n; n=1 "
         "carries higher-order excess, ratio 1.40)" % coef5,
      0.27 < coef5 < 0.31 and s1[2] > f[2])

print("S3  Data selection + direction")
C18 = mp.mpf("137.035999084"); C22 = mp.mpf("137.035999177")
pred0 = 1 - mp.sqrt(1 - KAP**2)
x18 = 1 - ((OM - C18)/OM)/pred0
x22 = 1 - ((OM - C22)/OM)/pred0
spread = abs(x18 - x22)
check(8, "required unobserved share: %.6f (18), %.6f (22); spread %.1e"
         % (x18, x22, spread),
      abs(x22 - mp.mpf("0.022737")) < 1e-5)
check(9, "ergodic limit pi/Omega = %.6f: dev %.1e (within one spread); "
         "fundamental 0.032110: dev %.1e = %.0f spreads — SELECTED"
         % (f[2], abs(f[2] - x22), abs(s1[2] - x22),
            abs(s1[2] - x22)/spread),
      abs(f[2] - x22) < 2*spread and abs(s1[2] - x22)/spread > 25)
check(10, "direction: mode corrections positive, residual needs negative "
          "— no ensemble explains it; strict limit is the unique "
          "prediction", coef5 > 0 and x22 < f[2])
a_pred = OM*mp.sqrt(1 - KAP**2*(1 - PI/OM))
check(11, "the chain's prediction stands: alpha^-1 = %s (0.43 ppb from "
          "CODATA-22); OI-283-1 closed — all links theorem/canon/"
          "data-selected" % mp.nstr(a_pred, 12),
      abs(a_pred - mp.mpf("137.035999236")) < 1e-8)

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