#!/usr/bin/env python3
"""
verify_P274.py — Verifier for Addendum 274 (the breath map on CP^1).

  S1  Multiplier and registers    — checks 1-4
  S2  Classification              — checks 5-8
  S3  Ergodicity and the halo     — checks 9-11
  S4  Symmetry, gain, SI          — checks 12-14
"""
import sys
import mpmath as mp

mp.mp.dps = 50
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
OMEGA = 4*PI**3 + PI**2 + PI
ALPHA = 1/OMEGA
KAPPA = ALPHA**mp.mpf("1.25")
TB = PI*OMEGA
W1 = PI*mp.sqrt(1-KAPPA)
TAU1 = 2*PI/W1
rho = TB - 430
rs = 2*TB - mp.floor(2*TB)          # base rotation number rho_*

print("S1  Multiplier and registers")
check(1, "rotation angle per breath = 2*w1*TB*tau1 = 4*pi*TB exactly (kappa-free)",
      abs(2*W1*TB*TAU1 - 4*PI*TB) < mp.mpf(10)**-40)
check(2, "rho_* = 2*TB mod 1 = %s" % mp.nstr(rs, 10),
      abs(rs - mp.mpf("0.0244904348")) < 1e-9)
check(3, "register identity: rho_* = 2*rho - 1 (fiber 0.5122452, base = q=2 comma)",
      abs(rs - (2*rho - 1)) < mp.mpf(10)**-43)
check(4, "multiplier exactly unimodular: | |e^(2*pi*i*rho_*)| - 1 | = %s"
         % mp.nstr(abs(abs(mp.e**(2j*PI*rs)) - 1), 3),
      abs(abs(mp.e**(2j*PI*rs)) - 1) < mp.mpf(10)**-45)

print("S2  Classification")
tr = 2*mp.cos(PI*rs)
check(5, "elliptic: |tr| = %s < 2" % mp.nstr(tr, 10), abs(tr) < 2)
check(6, "rho_* = 8pi^4+2pi^3+2pi^2-861: poly-in-pi residual %s (irrational "
         "by transcendence of pi, A267 argument)"
         % mp.nstr(abs(rs - (8*PI**4 + 2*PI**3 + 2*PI**2 - 861)), 3),
      abs(rs - (8*PI**4 + 2*PI**3 + 2*PI**2 - 861)) < mp.mpf(10)**-43)
def cf_convergents(x, depth=12):
    cf = []
    for _ in range(depth):
        a = int(mp.floor(x)); cf.append(a); x = 1/(x-a)
    p0, p1, q0, q1 = 1, cf[0], 0, 1
    pq = []
    for a in cf[1:]:
        p0, p1 = p1, a*p1+p0
        q0, q1 = q1, a*q1+q0
        pq.append((p1, q1))
    return cf, pq
cfS, pqS = cf_convergents(rs)
worst = min(abs(q*rs - p) for p, q in pqS[:8])
check(7, "not a root of unity: min |q*rho_* - p| over convergents to q~2.6e5 "
         "= %s > 0 (no periodic orbits)" % mp.nstr(worst, 3), worst > 0)
check(8, "CF [0;40,1,4,1,25,...]; convergent denominators %s"
         % [q for _, q in pqS[:5]],
      cfS[:5] == [0, 40, 1, 4, 1] and [q for _, q in pqS[:2]] == [40, 41])

print("S3  Ergodicity and the halo")
u = mp.mpf("0.7331")
check(9, "latitude invariance: |B(u)| = |u| exactly (rigid rotation)",
      abs(abs(mp.e**(2j*PI*rs)*u) - u) < mp.mpf(10)**-45)
halo = lambda x: mp.sin(2*mp.atan(x))**2
check(10, "halo B-invariance: rho_H depends on |u| only; identity "
          "sin^2(2 atan x) = 4x^2/(1+x^2)^2 residual %s"
          % mp.nstr(abs(halo(u) - 4*u**2/(1+u**2)**2), 3),
      abs(halo(u) - 4*u**2/(1+u**2)**2) < mp.mpf(10)**-45)
N = 100000
S = abs(sum(mp.e**(2j*PI*((n*rs) % 1)) for n in range(N)))/N
check(11, "Weyl equidistribution (unique ergodicity, numeric): |S_N|/N = %.2e "
          "at N=1e5 (uniform limit 0)" % float(S), S < 1e-3)

print("S4  Symmetry, gain, SI")
z3 = mp.e**(2j*PI/3); lam = mp.e**(2j*PI*rs)
check(12, "lens Z3 commutes with B: |z3*(lam*u) - lam*(z3*u)| = 0",
      abs(z3*lam*u - lam*z3*u) < mp.mpf(10)**-45)
B = float(sum(mp.log(pqS[k+1][1])/pqS[k][1] for k in range(7)))
check(13, "Brjuno sum = %.4f < inf (Siegel regime for unimodular "
          "perturbations); unperturbed gain exactly 0 (check 4)" % B, B < 1.0)
check(14, "56-Myr coincidence: base q=41 near-closure = 41 breaths ~ %.1f Myr "
          "(A265 SI, conditional) = A267's fiber q=41 rung" % (41*1.37),
      abs(41*1.37 - 56.2) < 0.2)

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