#!/usr/bin/env python3
"""
verify_P267.py — Verifier for Addendum 267 (rotation number of the breath).

  S1  Rotation number + irrationality          — checks 1-3
  S2  Derivation of 432 (3-smooth scan)        — checks 4-6
  S3  Convergent ladder (corrected recurrence) — checks 7-11
  S4  Near-tritone + SI                         — checks 12-14

Copyright: Leon Fernando Vlegels - MIT
"""
import sys
import mpmath as mp

mp.mp.dps = 60
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
TB = PI*OMEGA
rho = TB - 430

print("S1  Rotation number")
check(1, "rho = pi*Omega - 430 = %s (dps=60)" % mp.nstr(rho, 10),
      abs(rho - mp.mpf("0.512245217")) < 1e-8)
check(2, "0 < rho < 1 and irrational (transcendence of pi; structural)", 0 < rho < 1)
check(3, "Necessity of Detuning: no q with q*rho integral (spot q<=10^4)",
      all(abs(q*rho - mp.nint(q*rho)) > 1e-6 for q in range(1, 10001)))

print("S2  Derivation of 432")
smooth = sorted(2**a*3**b for a in range(12) for b in range(8) if 200 < 2**a*3**b < 800)
dists = {s: abs(s - float(TB)) for s in smooth}
nearest = min(dists, key=dists.get)
check(4, "nearest 3-smooth to T_b is %d (neighbours 384, 486)" % nearest, nearest == 432)
check(5, "margin: d(432)=%.3f vs next-best d(%d)=%.1f — order of magnitude"
      % (dists[432], 486, dists[486]), dists[486]/dists[432] > 10)
check(6, "G1 = 432/T_b - 1 (A266 identity, re-confirmed)",
      abs((432/TB - 1) - mp.mpf("0.0034557781")) < 1e-9)

print("S3  Convergent ladder (standard recurrence p_k = a_k p_{k-1} + p_{k-2})")
x, cf = rho, []
for _ in range(9):
    a = int(mp.floor(x)); cf.append(a); x = 1/(x - a)
p0, p1, q0, q1 = 1, cf[0], 0, 1     # p_{-1}=1, p_0=a_0; q_{-1}=0, q_0=1
conv = []
for a in cf[1:]:
    p0, p1 = p1, a*p1 + p0
    q0, q1 = q1, a*q1 + q0
    conv.append((p1, q1))
check(7, "CF(rho) = %s" % cf, cf[:5] == [0, 1, 1, 19, 1])
qs = [q for _, q in conv]
check(8, "convergent denominators begin 1, 2, 39, 41, 449: %s" % qs[:5],
      qs[:5] == [1, 2, 39, 41, 449])
commas = [float(abs(q*rho - mp.nint(q*rho))) for q in qs[:5]]
check(9, "comma(q=2) = %.4f" % commas[1], abs(commas[1] - 0.0245) < 0.001)
check(10, "comma(q=41) = %.4f" % commas[3], abs(commas[3] - 0.0020) < 0.0005)
check(11, "commas strictly decrease along convergents (q=1->2->39->41->449)",
      all(commas[i+1] < commas[i] for i in range(4)))

print("S4  Near-tritone and SI")
check(12, "rho - 1/2 = %.4f: near-antipodal breath (tritone structure)" % float(rho - 0.5),
      0 < float(rho - 0.5) < 0.02)
TB_MYR = 1.373
check(13, "2-breath near-cycle = %.1f Myr (conditional, A265)" % (2*TB_MYR),
      2.6 < 2*TB_MYR < 2.9)
check(14, "41-breath near-closure = %.0f Myr (conditional)" % (41*TB_MYR),
      53 < 41*TB_MYR < 60)

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