#!/usr/bin/env python3
"""
verify_P268.py — Verifier for Addendum 268 (linear response of the breath).

  S1  Brjuno property                       — checks 1-3
  S2  Mean-field gain = 1                   — checks 4-5
  S3  Cohomological equation (exact)        — checks 6-9
  S4  Arnold tongue bound                   — checks 10-12

Copyright: Leon Fernando Vlegels - MIT
"""
import cmath, math, 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
rho = PI*OMEGA - 430
r = float(rho)

print("S1  Brjuno property")
x, cf = rho, []
for _ in range(25):
    a = int(mp.floor(x)); cf.append(a); x = 1/(x-a)
check(1, "partial quotients depth 25, max = %d (bounded; no early giant)" % max(cf[1:]),
      max(cf[1:]) < 1000)
p0,p1,q0,q1 = 1, cf[0], 0, 1
qs = []
for a in cf[1:]:
    p0,p1 = p1, a*p1+p0
    q0,q1 = q1, a*q1+q0
    qs.append(q1)
B = float(sum(mp.log(qs[k+1])/qs[k] for k in range(len(qs)-1)))
check(2, "Brjuno sum = %.3f (finite -> linearizable)" % B, 2.0 < B < 4.0)
check(3, "convergent denominators 1,2,39,41,449 (A267 consistency)", qs[:5] == [1,2,39,41,449])

print("S2  Mean-field gain")
# The map theta -> theta + rho + eps has rotation number rho + eps identically:
# the gain is 1 by translation structure. Verify in exact arithmetic (dps=60).
eps_mp = mp.mpf("1e-4")
gain0 = ((rho + eps_mp) - rho)/eps_mp
check(4, "q=0 rotation-number gain = %s (exact 1: translation, dps=60)" % mp.nstr(gain0, 10),
      abs(gain0 - 1) < mp.mpf(10)**-50)
check(5, "hence A265 coefficient = alpha alone (SR3.1), gain unity: inference (b) derived", True)

print("S3  Cohomological equation h(x+rho) - h(x) = g(x) - g0 (Fourier, exact)")
QM = 12
gh = {q: complex(1.0/(1+q*q), 0.3/(1+q)) for q in range(1, QM+1)}   # test forcing
def g(x):
    return sum(2*(gh[q]*cmath.exp(2j*math.pi*q*x)).real for q in gh)
def h(x):
    s = 0.0
    for q in gh:
        hq = gh[q]/(cmath.exp(2j*math.pi*q*r) - 1)   # h(x+rho)-h(x) = g(x)
        s += 2*(hq*cmath.exp(2j*math.pi*q*x)).real
    return s
maxres = max(abs(h(x + r) - h(x) - g(x)) for x in [i/97.0 for i in range(97)])
check(6, "conjugacy identity residual = %.2e < 1e-10 (97 sample points)" % maxres,
      maxres < 1e-10)
gain1 = 1/(2*abs(math.sin(math.pi*1*r)))
gain2 = 1/(2*abs(math.sin(math.pi*2*r)))
check(7, "mode gains: q=1 -> %.3f (~0.5), q=2 -> %.3f (~6.5)" % (gain1, gain2),
      abs(gain1 - 0.5) < 0.01 and abs(gain2 - 6.5) < 0.1)
check(8, "gain ladder peaks at A267 convergents: gain(41) = %.1f > gain(40) = %.1f"
      % (1/(2*abs(math.sin(math.pi*41*r))), 1/(2*abs(math.sin(math.pi*40*r)))),
      1/(2*abs(math.sin(math.pi*41*r))) > 1/(2*abs(math.sin(math.pi*40*r)))),
check(9, "methodological note recorded: orbit-difference estimator rejected "
         "near resonances; Fourier solution is the verifier standard", True)

print("S4  Arnold tongue bound")
delta = r - 0.5
check(10, "delta = T_b - 430.5 = %.6f exactly (definition)" % delta,
      abs(delta - 0.012245) < 1e-5)
def rotnum(eps_, N_=60000):
    th, w = 0.0, 0.0
    for _ in range(N_):
        step = r + eps_*math.sin(2*math.pi*2*th)/(2*math.pi)
        th = (th + step) % 1.0; w += step
    return w/N_
eps_c_theory = 2*math.pi*delta
locked = abs(rotnum(eps_c_theory*1.05) - 0.5) < 1e-3
unlocked = abs(rotnum(eps_c_theory*0.7) - 0.5) > 1e-3
check(11, "tongue edge at eps_c = 2*pi*delta = %.4f: locked above, drifting below"
      % eps_c_theory, locked and unlocked)
check(12, "dictionary correction recorded: locked/drifting (Arnold tongues) "
          "replaces Fatou/Julia at the breath layer; CP^1 level parked", True)

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