#!/usr/bin/env python3
"""
PROVENANCE: NUMERICAL

Computation 213 -- the two-monomial bridge constant is a limit after all: a closed form,
approached from below as N^(-1/3) through the diagonal chain's Airy edge
==============================================================================
WHAT THIS DOES.  comp 206 followed the bridge's deficit (2 sqrt k - ratio(k, N)) N to
N = 15360 and found no plateau: the value rose slowly, 0.85 to 1.02 at k = 2, and the
data could not decide between a limit approached slowly and a logarithm.  The backlog
named an Airy-type edge analysis as the one untried route and warned there might be no
constant to derive.  This file finds the constant, and finds why comp 206 could not see
it.  The deficit times N is the sum of three terms, each with its own asymptotic law,
and comp 206 was watching the difference of two of them decay.

THE DECOMPOSITION.  ratio(k, N) = sqrt(lambda_C / lambda_T) with lambda the top eigenvalue
of the extremal chain's tridiagonal M^T M (comp 157).  Each chain has a symbol
S_t = d_t + o_t + o_(t-1) (the Rayleigh quotient of a constant vector), whose maximum
over the lattice is the operator's would-be norm if the eigenvector could sit still.
Writing S_C and S_T for the two symbol maxima and taking logarithms,
   deficit x N  =  W(N) + Z_C(N) - Z_T(N)  +  O(1/N),
   W   = -2 sqrt k . N . log( sqrt(S_C / S_T) / (2 sqrt k) )      the weights' 1/N term
   Z_C =    sqrt k . N . log( S_C / lambda_C )                     the commutator chain's zero point
   Z_T =    sqrt k . N . log( S_T / lambda_T )                     the diagonal chain's zero point
Measured directly, the identity holds to 3e-4 at every (k, N) computed.

THE THREE LAWS.  The Bergman weight w(p, q)^2 = (a+p)!(b+q)!(a+b+2)! / (a! b! (a+b+p+q+2)!)
along a chain of degree s = a + b and slope r = a/s is r^p (1-r)^q (1 + e(p, q, r)/s + ...)
with e = p(p+1)/(2r) + q(q+1)/(2(1-r)) - (p+q)(p+q+5)/2, negative, so every weight
increases toward the truncation and the symbol's maximum sits on the top degree s = N.
  (W)  The weights' term converges to a constant, in closed form from the e's: with
       c = (u e_u + v e_v)/(u + v) for each operator's two shift weights, plus a drift
       term for the commutator, gamma_C = -(u v / S_C) (2 r* - 1) d/dr log(u/v), which is
       there because o_t = v_t u_(t+1) couples neighbouring lattice points whose weights
       differ at order 1/N (the diagonal chain has none, u/v being constant along it):
       W_inf = -sqrt k (c_C + gamma_C - c_T).  At k = 2: c_T = -6.8079, c_C = -6.8993,
       gamma_C = -0.1000, W_inf = 0.2707; the measured W is 0.2708 from N = 7680 on.
  (Z_C)  The commutator's symbol has an interior maximum in r, at r* = 0.7143 for k = 2
       (the root of the symbol's r-derivative), so the extremal chain is the one through
       the corner (r*, N), which starts at a_0 = (2 r* - 1) N: 0.4286 N at k = 2, the
       position comp 206 saw the chain drifting toward.  Along it the potential is
       quadratic with curvature S'' (2r* - 1)^2 / N^2 and the hopping is o = u v, so the
       top eigenvector is a harmonic ground state of width ~ sqrt N, and its zero point
       is Z_C,inf = sqrt k (2 r* - 1) sqrt(o S'' / 2) / S_C: 0.9556 at k = 2, 1.1784 at
       k = 6.  Checked without the wall by extending the chain past N (the truncation moved
       to 2N, the symbol's maximum then interior): the free zero point converges to the
       harmonic constant with an O(1/N) residual, (Z_free - Z_inf) N staying between 19
       and 21 at k = 2 and between 74 and 78 at k = 6 across three doublings.
       With the wall, the eigenvector keeps a distance ~ sqrt(N log N) from it, trading a
       linear loss of symbol against the exponential cost of the truncation, and the excess
       E(N) = Z_C(N) - Z_C,inf decays like N^(-1/2) sqrt(log N): its ratio per doubling from
       N = 7680 is 0.734, 0.732 at k = 2 (0.731 at the offline N = 61440) against the law's
       0.736, 0.734 (0.732), and 0.743, 0.741 at k = 6.
  (Z_T)  The diagonal chain a = b has r = 1/2 fixed, so its potential rises monotonically
       to the truncation with slope 2 S_T |c_T| / N^2 per step, and its top eigenvector is
       an Airy function pressed against the hard edge with a Dirichlet condition: the zero
       point is Z_T = A(k) N^(-1/3) with A = sqrt k . a_1 . (o_T slope^2)^(1/3) / S_T and
       a_1 = 2.33811 the first zero of Ai.  A(2) = 10.97; the measured Z_T N^(1/3) is
       11.26, 11.30, 11.29, 11.26, 11.22, 11.18, 11.15 for N = 480 to 30720, 1.6 percent
       above A and falling toward it.  This is the Airy edge the backlog named, and it
       sits on the bridge's own norm, not on the commutator's.

RESULT.  deficit x N = D_inf(k) - A(k) N^(-1/3) + E(N) + ..., with
   D_inf(k) = W_inf(k) + Z_C,inf(k)  =  1.2263 (k = 2),  1.4738 (k = 4),  1.5235 (k = 6),
all in closed form given r*.  The approach is from below at rate N^(-1/3): at N = 61440,
k = 2, the measured 1.072 is still 0.15 short of its limit, and the per-doubling increments
comp 206 reported (0.030 at N = 15360, k = 2) are the Airy term's decay minus the wall
excess's decay, 0.0928 - 0.0640 = 0.029.  There was never a logarithm.
   The closed form over k: every k from 2 to 1000 closes (ratio_inf = 2 sqrt k to 1e-10
with comp 59's alpha(k)); D_inf(k) is 2.4648 at most for k <= 80 (at k = 72) and below
the certified 3.6736 for every k <= 482; its envelope grows like k^0.22 (the harmonic
term is ~ sqrt(k/m)), crossing 3.6736 first at k = 483, where the asymptotic regime needs
N far beyond (A/D_inf)^3 ~ 1e8.  So the uniform bound of eq. gap-uniform-bound holds as
N -> infinity for every fixed k up to 482, with the limit well below the finite-N sharp
value, and the closure rate gamma_D = O(1/sqrt D) of Lemma 7 stands on that range without
the logarithm caveat; it is not uniform in k without bound, the limit constant growing
slowly with k, so beyond k of a few hundred, at substrate sizes no computation reaches,
the constant in Lemma 7 grows at most like k_D^(1/4) = D^(1/8).

Tampers: (i) the weights' constant without the drift term misses the measured W by 0.14
(k = 2) and 0.06 (k = 6) where the full form matches to 0.001; (ii) the Neumann Airy zero
(1.0188) in place of the Dirichlet one predicts a Z_T coefficient 2.3 times too small;
(iii) a harmonic well with the wall at its centre (zero point 3/2 h-bar omega) is rejected
by the free-chain check by a factor three.

Run offline at k = 12 and k = 20 to N = 61440 the identity, the W law (0.3609 against 0.3607;
0.3671 against 0.3666) and the Airy law (58.0 against 57.5; 90.2 against 89.5) hold the same way.

Runs in about half a minute; imports comp 206 (and through it comp 59) beside it.
"""
import math
import numpy as np
from scipy.linalg import eigvalsh_tridiagonal

import computation_206_se as c206

NS = (480, 960, 1920, 3840, 7680, 15360, 30720)
AIRY1 = 2.33810741045977          # first zero of Ai, Dirichlet edge
AIRY1_NEUMANN = 1.01879297164747  # first zero of Ai', the tamper


# ---------------------------------------------------------------- the chains
def chain(N, spec, a0, b0, mask_N=None):
    """Top eigenvalue and symbol of the chain (a0 + t, b0 + t) truncated at mask_N (default N)."""
    u, v = c206.chain_vec(N if mask_N is None else mask_N, spec, a0, b0)
    n = len(u)
    if n == 0 or u.max() == 0:
        return 0.0, np.zeros(0)
    d = u * u + v * v
    if n == 1:
        return float(d[0]), d
    o = v[:-1] * u[1:]
    lam = eigvalsh_tridiagonal(d, o, select='i', select_range=(n - 1, n - 1))[0]
    sym = d.copy(); sym[:-1] += o; sym[1:] += o
    return float(lam), sym


# ---------------------------------------------------------- the closed forms
def uv(spec, r):
    c1, p1, q1, c2, p2, q2 = spec
    return c1 * r ** (p1 / 2) * (1 - r) ** (q1 / 2), c2 * r ** (p2 / 2) * (1 - r) ** (q2 / 2)


def symbol(spec, r):
    u, v = uv(spec, r)
    return (u + v) ** 2


def e_coef(p, q, r):
    """First 1/s coefficient of the weight: w^2 = r^p (1-r)^q (1 + e/s + ...)."""
    return p * (p + 1) / (2 * r) + q * (q + 1) / (2 * (1 - r)) - (p + q) * (p + q + 5) / 2


def r_star(spec):
    """Maximiser of the commutator symbol on (1/2, 1): grid, then parabolic refinement."""
    r = np.linspace(0.5 + 1e-6, 1 - 1e-6, 20001)
    f = symbol(spec, r)
    i = int(np.argmax(f)); x = r[i]; h = r[1] - r[0]
    for _ in range(60):
        f0, fp, fm = symbol(spec, x), symbol(spec, x + h), symbol(spec, x - h)
        den = fp - 2 * f0 + fm
        if den >= 0:
            break
        x = x - h * (fp - fm) / (2 * den)
        h /= 2
        if h < 1e-13:
            break
    return float(x)


def constants(k, drift=True, airy_zero=AIRY1):
    m, al, sT, sC = c206.specs(k)
    rs = r_star(sC)
    SC, ST = symbol(sC, rs), symbol(sT, 0.5)
    c1, p1, q1, c2, p2, q2 = sC; u, v = uv(sC, rs)
    cC = (u * e_coef(p1, q1, rs) + v * e_coef(p2, q2, rs)) / (u + v)
    dlog_uv = (p1 - p2) / (2 * rs) - (q1 - q2) / (2 * (1 - rs))      # d/dr log(u/v)
    gC = -(u * v / SC) * (2 * rs - 1) * dlog_uv if drift else 0.0
    c3, p3, q3, c4, p4, q4 = sT; uT, vT = uv(sT, 0.5)
    cT = (uT * e_coef(p3, q3, 0.5) + vT * e_coef(p4, q4, 0.5)) / (uT + vT)
    W = -math.sqrt(k) * (cC + gC - cT)
    h = 1e-4; Spp = (symbol(sC, rs + h) - 2 * SC + symbol(sC, rs - h)) / h ** 2
    ZC = math.sqrt(k) * (2 * rs - 1) * math.sqrt(u * v * (-Spp) / 2) / SC
    slope = 2 * ST * abs(cT)
    A = math.sqrt(k) * airy_zero * (uT * vT * slope ** 2) ** (1 / 3) / ST
    return dict(m=m, al=al, rs=rs, SC=SC, ST=ST, cT=cT, cC=cC, gC=gC, W=W, ZC=ZC, A=A,
                D=W + ZC, ratio=math.sqrt(SC / ST))


# ---------------------------------------------------------- the measurements
def decompose(k, N, a0_prev):
    """D, W, Z_C, Z_T at (k, N); the commutator chain scanned fully at N = 480, in a window after."""
    m, al, sT, sC = c206.specs(k); tgt = 2 * math.sqrt(k)
    lamT, symT = 0.0, None
    for st in ((0, 0), (1, 0), (2, 0), (0, 1)):
        l, s = chain(N, sT, *st)
        if l > lamT:
            lamT, symT = l, s
    ST = float(symT.max())
    if a0_prev is None:
        starts = range(0, N)
    else:
        w = max(8, N // 64); starts = range(max(0, 2 * a0_prev - w), 2 * a0_prev + w + 1)
    lamC, a0, SC = 0.0, None, 0.0
    for a in starts:
        l, s = chain(N, sC, a, 0)
        if l > lamC:
            lamC, a0 = l, a
        if len(s) and s.max() > SC:
            SC = float(s.max())
    ratio = math.sqrt(lamC / lamT); ratio_sym = math.sqrt(SC / ST)
    return dict(N=N, a0=a0, D=(tgt - ratio) * N,
                W=-tgt * N * math.log(ratio_sym / tgt),
                ZC=math.sqrt(k) * N * (-math.log(lamC / SC)),
                ZT=math.sqrt(k) * N * (-math.log(lamT / ST)))


def free_zero_point(k, N, rs):
    """The commutator chain through the corner with the truncation moved to 2N: no wall."""
    m, al, sT, sC = c206.specs(k)
    a0 = int(round((2 * rs - 1) * N))
    lam, sym = chain(N, sC, a0, 0, mask_N=2 * N)
    return math.sqrt(k) * N * (-math.log(lam / sym.max()))


def wall_law_ratio(N):
    """Ratio E(2N)/E(N) for E ~ N^(-1/2) sqrt(log N)."""
    return math.sqrt(0.5 * math.log(2 * N) / math.log(N))


def main():
    print("=" * 100)
    print("Computation 213: the two-monomial bridge constant, in closed form, approached as N^(-1/3)")
    print("=" * 100)
    checks = {}
    rows = {}; cons = {}
    # ---- Part 1: the decomposition, measured
    print("PART 1.  deficit x N = W + Z_C - Z_T, measured on the extremal chains")
    for k in (2, 6):
        cons[k] = constants(k); a0 = None
        print(f"  k = {k} (m = {cons[k]['m']}, alpha = {cons[k]['al']:.4f}, r* = {cons[k]['rs']:.4f}):")
        print(f"      N      a0/N     D        W       Z_C      Z_T    W+Z_C-Z_T   Z_T N^(1/3)   E = Z_C - Z_C,inf")
        for N in NS:
            r = decompose(k, N, a0); a0 = r['a0']; rows[(k, N)] = r
            r['E'] = r['ZC'] - cons[k]['ZC']
            print(f"  {N:7d}   {a0 / N:.4f}  {r['D']:.4f}  {r['W']:.4f}  {r['ZC']:.4f}  {r['ZT']:.4f}   {r['W'] + r['ZC'] - r['ZT']:.4f}     {r['ZT'] * N ** (1 / 3):.4f}        {r['E']:+.4f}")
    checks['identity'] = max(abs(r['D'] - (r['W'] + r['ZC'] - r['ZT'])) for r in rows.values()) < 2e-3
    print(f"  the identity holds to {max(abs(r['D'] - (r['W'] + r['ZC'] - r['ZT'])) for r in rows.values()):.1e} at every (k, N): {checks['identity']}")
    print()
    # ---- Part 2: the three laws against the closed forms
    print("PART 2.  the three terms against their closed forms")
    okW = okA = okF = okE = okI = True
    for k in (2, 6):
        c = cons[k]
        print(f"  k = {k}: c_T = {c['cT']:.4f}, c_C(r*) = {c['cC']:.4f}, gamma_C = {c['gC']:.4f}  ->  W_inf = {c['W']:.4f};"
              f"  harmonic Z_C,inf = {c['ZC']:.4f};  Airy A = {c['A']:.3f};  D_inf = W_inf + Z_C,inf = {c['D']:.4f}")
        dW = max(abs(rows[(k, N)]['W'] - c['W']) for N in NS if N >= 7680)
        okW &= dW < 2e-3
        print(f"      (W)   measured W at N >= 7680 within {dW:.4f} of W_inf")
        q = [rows[(k, N)]['ZT'] * N ** (1 / 3) / c['A'] for N in NS]
        okA &= 1.0 < q[-1] < 1.03 and q[-1] < q[-2] < q[-3]
        print(f"      (Z_T) Z_T N^(1/3) / A = {', '.join(f'{x:.4f}' for x in q)}: above one and falling toward it")
        fr = [(N, free_zero_point(k, N, c['rs'])) for N in (1920, 3840, 7680, 15360)]
        resid = [(z - c['ZC']) * N for N, z in fr]
        okF &= abs(fr[-1][1] - c['ZC']) < 0.01 and max(abs(x) for x in resid) / min(abs(x) for x in resid) < 1.3
        print(f"      (Z_C) free chain (no wall): Z_free = {', '.join(f'{z:.4f}' for _, z in fr)} at N = 1920 to 15360;"
              f" (Z_free - Z_C,inf) N = {', '.join(f'{x:.1f}' for x in resid)}, an O(1/N) residual")
        E = [rows[(k, N)]['E'] for N in NS]
        rat = [E[i + 1] / E[i] for i in range(len(NS) - 1)]
        law = [wall_law_ratio(N) for N in NS[:-1]]
        okE &= all(abs(rat[i] - law[i]) < 0.02 for i in range(len(rat) - 2, len(rat)))
        print(f"      (wall) E ratio per doubling from N = 7680: {', '.join(f'{x:.3f}' for x in rat[-2:])} against N^(-1/2) sqrt(log N)'s {', '.join(f'{x:.3f}' for x in law[-2:])}")
        for i in (len(NS) - 3, len(NS) - 2):
            N = NS[i]; r0, r1 = rows[(k, N)], rows[(k, NS[i + 1])]
            pred = r0['E'] * (wall_law_ratio(N) - 1) + r0['ZT'] * (1 - 2 ** (-1 / 3))
            meas = r1['D'] - r0['D']
            okI &= abs(pred - meas) < 0.015
            print(f"      increment N = {N} -> {2 * N}: measured {meas:+.4f}, the two laws give {pred:+.4f} (Airy decay {r0['ZT'] * (1 - 2 ** (-1 / 3)):+.4f}, wall decay {r0['E'] * (wall_law_ratio(N) - 1):+.4f})")
        below = all(rows[(k, N)]['D'] < c['D'] for N in NS)
        checks[f'below_{k}'] = below
        print(f"      every measured deficit x N lies below D_inf = {c['D']:.4f}: {below};  at N = 30720 the shortfall is {c['D'] - rows[(k, 30720)]['D']:.4f}")
    checks['W'] = okW; checks['airy'] = okA; checks['free'] = okF; checks['wall'] = okE; checks['increments'] = okI
    print()
    # ---- Part 3: the limit over k
    print("PART 3.  the limit constant D_inf(k) = W_inf + Z_C,inf over k")
    ks = list(range(2, 1001)); D = []; closed = True
    for k in ks:
        c = constants(k); D.append(c['D'])
        closed &= abs(c['ratio'] - 2 * math.sqrt(k)) < 1e-8
        if k in (2, 3, 4, 6, 8, 10, 12, 16, 20, 40, 80):
            print(f"  k = {k:3d}: m = {c['m']:2d}  alpha = {c['al']:9.4f}  r* = {c['rs']:.4f}  W_inf = {c['W']:.4f}  Z_C,inf = {c['ZC']:.4f}  D_inf = {c['D']:.4f}  A = {c['A']:7.2f}")
    D = np.array(D); ks = np.array(ks); env = np.maximum.accumulate(D)
    first = int(ks[np.argmax(D > 3.6736)]) if (D > 3.6736).any() else None
    sel = ks >= 20; expo = np.polyfit(np.log(ks[sel]), np.log(env[sel]), 1)[0]
    print(f"  every k from 2 to 1000 closes, ratio_inf = 2 sqrt k to 1e-8: {closed}")
    print(f"  max D_inf over k <= 80: {D[ks <= 80].max():.4f} at k = {int(ks[ks <= 80][D[ks <= 80].argmax()])}; median {np.median(D[ks <= 80]):.3f}")
    print(f"  first k with D_inf above the certified 3.6736: {first}; envelope ~ k^{expo:.2f} for k >= 20;"
          f" envelope at k = 160, 320, 640, 1000: {', '.join(f'{env[ks == x][0]:.3f}' for x in (160, 320, 640, 1000))}")
    checks['closed'] = closed
    checks['bound80'] = D[ks <= 80].max() < 3.6736
    checks['first'] = first is not None and 400 < first < 600 and bool((D[ks < first] < 3.6736).all())
    checks['envelope'] = 0.15 < expo < 0.30
    print()
    # ---- Tampers
    print("TAMPERS")
    cW = [constants(k, drift=False)['W'] for k in (2, 6)]
    t1 = all(abs(cW[i] - rows[(k, 30720)]['W']) > 0.05 for i, k in enumerate((2, 6)))
    print(f"  (i)   W_inf without the drift term: {cW[0]:.4f}, {cW[1]:.4f} against the measured {rows[(2, 30720)]['W']:.4f}, {rows[(6, 30720)]['W']:.4f}: rejected: {t1}")
    AN = constants(2, airy_zero=AIRY1_NEUMANN)['A']
    t2 = rows[(2, 30720)]['ZT'] * 30720 ** (1 / 3) / AN > 2
    print(f"  (ii)  Neumann Airy zero: A = {AN:.3f} against the measured Z_T N^(1/3) = {rows[(2, 30720)]['ZT'] * 30720 ** (1 / 3):.3f}: rejected: {t2}")
    zf = free_zero_point(2, 15360, cons[2]['rs'])
    t3 = abs(zf / (3 * cons[2]['ZC']) - 1) > 0.5
    print(f"  (iii) half-well zero point 3 Z_C,inf = {3 * cons[2]['ZC']:.4f} against the free chain's {zf:.4f}: rejected: {t3}")
    checks['tampers'] = t1 and t2 and t3
    print()
    ok = all(checks.values())
    print("SUMMARY")
    print(f"  [{'PASS' if checks['identity'] else 'FAIL'}] deficit x N = W + Z_C - Z_T at every (k, N) computed")
    print(f"  [{'PASS' if checks['W'] else 'FAIL'}] the weights' term converges to its closed form W_inf, the drift term included")
    print(f"  [{'PASS' if checks['airy'] else 'FAIL'}] the diagonal chain's zero point is the Airy edge, A(k) N^(-1/3)")
    print(f"  [{'PASS' if checks['free'] else 'FAIL'}] the commutator chain's zero point is the harmonic constant once the wall is removed")
    print(f"  [{'PASS' if checks['wall'] and checks['increments'] else 'FAIL'}] the wall excess decays like N^(-1/2) sqrt(log N); comp 206's increments are the difference of the two decays")
    print(f"  [{'PASS' if checks['below_2'] and checks['below_6'] else 'FAIL'}] every measured value lies below its limit: the approach is from below")
    print(f"  [{'PASS' if checks['closed'] and checks['bound80'] and checks['first'] and checks['envelope'] else 'FAIL'}] D_inf(k) in closed form for k <= 1000: below 2.47 for k <= 80, below 3.6736 for k <= 482, envelope ~ k^0.22")
    print(f"  [{'PASS' if checks['tampers'] else 'FAIL'}] tampers rejected")
    print()
    print("  CONCLUSION: the bridge's 1/N term is a constant after all, D_inf(k) = W_inf + Z_C,inf in closed")
    print("  form (1.2263 at k = 2, 1.5235 at k = 6), approached from below as A(k) N^(-1/3) through the")
    print("  Airy edge of the diagonal chain, which is why no plateau was visible at N = 15360.  The")
    print("  uniform bound holds as N -> infinity for every fixed k up to 482, so the closure rate of Lemma 7")
    print("  stands without a logarithm; the limit constant grows slowly with k, so the bound is not uniform")
    print("  in k without end, at substrate sizes no computation reaches.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()