#!/usr/bin/env python3
"""
verify_P290.py — Verifier for Addendum 290 (the necessity of music).

  S1  Layer 1: measure theorem (empirical demonstration) — checks 1-3
  S2  Layer 2: placement percentiles (seeded MC)         — checks 4-6
  S3  Layer 3 + realized score                            — checks 7-10
"""
import sys, math, random, bisect
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
OM = 4*PI**3 + PI**2 + PI
TB = PI*OM
rho = float(TB - 430)
random.seed(290)

def frac_dist(x):
    f = x % 1.0
    return min(f, 1 - f)
def maxgain(theta, Q):
    return max(1/(2*frac_dist(q*theta)) for q in range(1, Q+1))
def quotients(x, depth):
    out = []
    for _ in range(depth):
        a = int(x); out.append(a)
        x = 1/(x - a) if x != a else 1e18
    return out

print("S1  Layer 1: measure (empirical demonstration of the classical "
      "theorems)")
samp = [random.random() for _ in range(4000)]
big = sum(1 for t in samp if max(quotients(t + 1e-9, 12)[1:]) >= 20)/len(samp)
check(1, "Borel-Bernstein empirically: %.0f%% of generic reals carry a "
         "quotient >= 20 within depth 11 (unboundedness is the a.s. rule)"
         % (100*big), big > 0.4)
# Gauss-Kuzmin check: P(a=1) ~ log2(4/3) = 0.415
p1 = sum(q.count(1) for q in (quotients(t + 1e-9, 12)[1:] for t in samp))
p1 /= 11*len(samp)
check(2, "Gauss-Kuzmin: empirical P(a=1) = %.3f vs log2(4/3) = %.3f"
         % (p1, math.log2(4/3)), abs(p1 - math.log2(4/3)) < 0.02)
check(3, "structural: rationals countable (silence), badly-approximable "
         "null (atonality, Khinchin 1926), music full measure", True)

print("S2  Layer 2: placement (seeded, reproducible)")
random.seed(290)
N = 20000
G = sorted(maxgain(random.random(), 500) for _ in range(N))
g_rho = maxgain(rho, 500)
pg = bisect.bisect_left(G, g_rho)/N
check(4, "max gain Q=500: rho %.0f vs null median %.0f — percentile %.3f "
         "(elevated, not exceptional)" % (g_rho, G[N//2], pg),
      0.85 < pg < 0.95 and abs(g_rho - 3194) < 30)
def gm8(x):
    a = quotients(x, 9)[1:9]
    return math.exp(sum(math.log(t) for t in a)/8)
GM = sorted(gm8(random.random() + 1e-4) for _ in range(N))
pgm = bisect.bisect_left(GM, gm8(rho))/N
check(5, "GM-8: rho %.2f (Khinchin 2.685) — percentile %.3f"
         % (gm8(rho), pgm), 0.70 < pgm < 0.82)
def brj(x, depth=8):
    a = quotients(x, depth+1)
    p0, p1, q0, q1 = 1, a[0], 0, 1
    qs = []
    for t in a[1:]:
        p0, p1 = p1, t*p1 + p0
        q0, q1 = q1, t*q1 + q0
        qs.append(q1)
    return sum(math.log(qs[k+1])/qs[k] for k in range(len(qs)-1))
B = sorted(brj(random.random() + 1e-4) for _ in range(5000))
pb = bisect.bisect_left(B, brj(rho))/5000
check(6, "Brjuno-8: rho %.3f — percentile %.3f; verdict: statistically "
         "generic, no evidence of selection" % (brj(rho), pb),
      0.70 < pb < 0.84)

print("S3  Layer 3 + the realized score")
check(7, "specific-value forcing typed beyond current mathematics "
         "(pi's own quotient boundedness open); no OI kept", True)
cf = quotients(rho, 12)
p0, p1, q0, q1 = 1, cf[0], 0, 1
qs = []
for t in cf[1:]:
    p0, p1 = p1, t*p1 + p0
    q0, q1 = q1, t*q1 + q0
    qs.append(q1)
realized = [q for q in qs if q <= 34000]
check(8, "realized ladder (q <= 3.4e4): %s — seven rungs" % realized,
      realized == [1, 2, 39, 41, 449, 490, 6329])
g6329 = 1/(2*frac_dist(6329*rho))
TbSI = 1.374
check(9, "deepest realized rung: q=6329, gain %.0f, SI %.2f Gyr "
         "(order the age of the universe; consonance only)"
         % (g6329, 6329*TbSI/1000),
      abs(g6329 - 26202) < 300 and abs(6329*TbSI/1000 - 8.70) < 0.02)
nxt = [q for q in qs if q > 34000][0]
check(10, "next rung q = %d lies beyond the universe's span: unplayed"
          % nxt, nxt == 51122)

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