#!/usr/bin/env python3
"""
verify_P289.py — Verifier for Addendum 289 (the score of the breath).

  S1  Instruments: discrepancy     — checks 1-3
  S2  The ladder                   — checks 4-8
  S3  The chord                    — checks 9-11
  S4  Thesis bookkeeping           — check  12
"""
import sys, math
import mpmath as mp

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
TB = PI*OM
rho = float(TB - 430)
phi = (math.sqrt(5) - 1)/2

def star_discrepancy(theta, N):
    pts = sorted((n*theta) % 1.0 for n in range(1, N+1))
    D = 0.0
    for i, x in enumerate(pts):
        D = max(D, abs((i+1)/N - x), abs(x - i/N))
    return D

def frac_dist(x):
    f = x % 1.0
    return min(f, 1 - f)

print("S1  Instruments: discrepancy (N = 1e4)")
N = 10000
Dr = star_discrepancy(rho, N); Dp = star_discrepancy(phi, N)
check(1, "D*(rho) = %.5f" % Dr, abs(Dr - 0.00075) < 5e-5)
check(2, "D*(phi) = %.5f (optimal-class)" % Dp, abs(Dp - 0.00026) < 5e-5)
check(3, "rho is ~3x worse-spread than phi (ratio %.2f): poor hash, "
         "for the tonal reason" % (Dr/Dp), 2.5 < Dr/Dp < 3.5)

print("S2  The ladder")
best = []; cur = 1.0
for q in range(1, 500):
    d = frac_dist(q*rho)
    if d < cur:
        best.append((q, d)); cur = d
qs = [q for q, d in best]
check(4, "record-gain rungs below 500: %s" % qs,
      qs == [1, 2, 39, 41, 449, 490])
gains = {q: 1/(2*d) for q, d in best}
check(5, "gains: q=2: %.1f, q=41: %.0f, q=490: %.0f — jumps, not creep"
         % (gains[2], gains[41], gains[490]),
      abs(gains[2] - 20.4) < 0.1 and abs(gains[41] - 243) < 2 and
      abs(gains[490] - 3194) < 30)
# 490 as sixth convergent
x = float(TB - 430); cf = []
for _ in range(7):
    a = int(math.floor(x)); cf.append(a); x = 1/(x - a)
p0, p1, q0, q1 = 1, cf[0], 0, 1
qden = []
for a in cf[1:]:
    p0, p1 = p1, a*p1 + p0
    q0, q1 = q1, a*q1 + q0
    qden.append(q1)
check(6, "q = 490 = 449 + 41 is the sixth convergent denominator: %s"
         % qden[:6], qden[:6] == [1, 2, 39, 41, 449, 490])
gphi = max(1/(2*frac_dist(q*phi)) for q in range(1, 500))
T = gains[490]/gphi
check(7, "tonality index T(rho,500) = %.1f (phi max gain %.0f)"
         % (T, gphi), 7.0 < T < 8.5 and abs(gphi - 421) < 5)
cents = lambda d: 1200*math.log2(1 + d)
check(8, "cents: q=2: %.2f (near-quarter-tone), q=41: %.2f "
         "(schisma-scale), q=490: %.2f; G1 = %.2f (432-comma)"
         % (cents(frac_dist(2*rho)), cents(frac_dist(41*rho)),
            cents(frac_dist(490*rho)), cents(float(432/TB - 1))),
      abs(cents(frac_dist(2*rho)) - 41.89) < 0.05 and
      abs(cents(frac_dist(41*rho)) - 3.55) < 0.02 and
      abs(cents(float(432/TB - 1)) - 5.97) < 0.02)

print("S3  The chord")
check(9, "(2,3,16): 3:2 = %.1f cents (perfect fifth), 16:3 = 2 octaves "
         "+ %.1f (perfect fourth), 16:2 = 3 octaves exactly"
         % (1200*math.log2(1.5), 1200*math.log2(16/3) - 2400),
      abs(1200*math.log2(1.5) - 702.0) < 0.1 and
      abs(1200*math.log2(16/3) - 2400 - 498.0) < 0.1 and
      16/2 == 8 and math.log2(8) == 3)
check(10, "A287 cross-reference: lepton corrections draw the chord tones "
          "(muon 3, tau 16, electron anchored at edge 2)", True)
check(11, "downstream echoes on file: 432 = 2^4 3^3; j(i) = 1728 = 4*432 "
          "(A266 flag, unpromoted)", 2**4*3**3 == 432 and 1728 == 4*432)

print("S4  Thesis bookkeeping")
TbSI = 1.374  # Myr, conditional (A265)
check(12, "SI rungs (conditional): q=2: %.1f Myr, q=41: %.0f Myr, "
          "q=490: %.0f Myr — stated without external identification; "
          "OI-289-1 filed (derive the tonality)"
          % (2*TbSI, 41*TbSI, 490*TbSI),
      abs(2*TbSI - 2.7) < 0.1 and abs(41*TbSI - 56) < 1 and
      abs(490*TbSI - 673) < 2)

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