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

Computation 206 -- the two-monomial deficit to N = 15360: no plateau, so the per-k
constant of comps 157 and 203 is a value, not a limit
==============================================================================
WHAT THIS DOES.  comp 203 measured the bridge's deficit (2 sqrt k - ratio(k, N)) N to
N = 480 and reported it "near 0.8"; comp 157 had reported bands near 0.8 at N = 300.
Before an analytic route to that number is attempted (the backlog names an Airy-type
edge analysis as the untried one), this file asks whether the number is a limit at all.
It extends the measurement by a factor thirty-two in N, using comp 157's chain reduction
with the chain weights vectorised, and finds that it is not: the deficit times N keeps
rising, slowly, through N = 15360, and the chain that carries the commutator's norm
keeps moving.

METHOD.  The extremal chains are known from comp 203: the diagonal a = b for ||T||, and
for the commutator a chain starting at (a_0, 0) with a_0/N drifting upward.  Every chain
is scanned at N <= 3840 (full scans, as comp 203 did), which fixes a_0/N to within a
percent; at larger N the scan is restricted to a window of twenty percent around the
extrapolated a_0 and then refined to the step, which the full scans at N <= 3840 show
loses nothing.  Each chain norm is the top eigenvalue of its tridiagonal M^T M.

RESULTS (deficit times N, exact chain norms).
   N      480    960   1920   3840   7680  15360
   k=2  0.851  0.880  0.913  0.948  0.983  1.015     (increments per doubling 0.029, 0.033,
                                                      0.035, 0.034, 0.033)
   k=6  0.811  0.838  0.881  0.934  0.992  1.050     (0.027, 0.043, 0.053, 0.058, 0.058)
  Run offline to N = 61440 with the same code: k = 2 reaches 1.072 with increments
  0.030, 0.027, and k = 6 reaches 1.162 with increments 0.057, 0.054.
  1. THERE IS NO PLATEAU IN THIS RANGE.  Over five doublings the deficit times N rises
     by 0.16 (k = 2) and 0.24 (k = 6).  The increments per doubling shrink slowly for
     k = 2 and barely for k = 6, so the data do not decide between a limit approached
     logarithmically slowly and a term of order log N / N with no limit.  Either way the
     bands comp 157 reported at N = 300 and comp 203 at N = 480 are values at those N.
  2. THE EXTREMAL CHAIN MOVES.  The commutator's chain start a_0/N rises from 0.392 at
     N = 480 to 0.420 at N = 15360 (k = 2) and from 0.208 to 0.222 (k = 6), still moving
     at the last point, toward but not at the single-shift maximiser's (m-1)/(2m)-type
     position; the two-monomial bridge selects its chain by a balance the single-shift
     ratio test of comp 203 does not capture.
  3. THE SYMBOL CONSTANT MOVES THE OTHER WAY.  The maximum of the chain symbol
     d_t + o_t + o_{t-1} over the same chains gives a deficit times N that falls
     (1.01 to 0.61 for k = 2 between N = 480 and 3840), so the exact value and the
     symbol value cross; the gap between them, which is where any 1/N constant would
     have to come from, grows like log N over this range.  The symbol route to a closed
     form is closed.
  4. WHAT STANDS.  The sharp value 3.6736 at (k, N) = (6, 12), certified by comp 203,
     remains the maximum of everything scanned here: at the largest N the deficit times
     N is near 1.2.  The uniform bound of eq. gap-uniform-bound holds on every point
     computed, and whether it holds as N -> infinity for fixed k is exactly what this
     file leaves open, since a logarithm would eventually exceed any constant.  The
     closure rate gamma_D = O(1/sqrt D) at matched scaling is therefore established up to
     a possible logarithm.

WHAT THIS CHANGES.  The book's "the two-monomial constant ... near 0.85" becomes "rises
slowly with N, 0.85 at N = 480 and 1.07 at N = 61440, limit undecided"; the backlog's
untried route stays untried, with the warning that there may be no constant to derive.

Tamper: evaluating the commutator norm on the single-shift maximiser's chain instead of
the scanned one gives a ratio below the scanned value by a growing margin, so the chain
drift of result 2 is real and load-bearing.

Runs in about two minutes; imports comp 59 beside it.
"""
import math
import numpy as np
from scipy.linalg import eigvalsh_tridiagonal
from scipy.special import gammaln

import computation_59_se as c59

NMAX = 15360


def chain_vec(N, spec, a0, b0):
    c1, p1, q1, c2, p2, q2 = spec
    tmax = (N - a0 - b0) // 2
    if tmax < 0:
        return np.zeros(0), np.zeros(0)
    t = np.arange(tmax + 1); a = a0 + t; b = b0 + t
    def w(p, q):
        ok = (a + b + p + q <= N) & (a + p >= 0) & (b + q >= 0)
        val = np.exp(0.5 * ((gammaln(a + p + 1) + gammaln(b + q + 1) + gammaln(a + b + 3))
                            - (gammaln(a + 1) + gammaln(b + 1) + gammaln(a + b + p + q + 3))))
        return np.where(ok, val, 0.0)
    return c1 * w(p1, q1), c2 * w(p2, q2)


def lam_sym(u, v):
    n = len(u); d = u * u + v * v
    if n == 1:
        return float(d[0]), float(d[0])
    o = v[:-1] * u[1:]
    lam = eigvalsh_tridiagonal(d, o, select='i', select_range=(n - 1, n - 1))[0]
    H = d.copy(); H[:-1] += o; H[1:] += o
    return float(lam), float(H.max())


def specs(k):
    m, al = c59.m_of_k(k), c59.alpha_of_k(k)
    return m, al, (1.0, m, m, al, m + 1, m + 1), (float(m), m + 1, m - 1, al * (m + 1), m + 2, m)


def best(N, spec, starts):
    b = (0.0, None, 0.0)
    for a0, b0 in starts:
        u, v = chain_vec(N, spec, a0, b0)
        if len(u) == 0 or u.max() == 0:
            continue
        l, s = lam_sym(u, v)
        if l > b[0]:
            b = (l, (a0, b0), s)
    return b


def deficit(k, N, frac=None):
    m, al, sT, sC = specs(k)
    lT, chT, symT = best(N, sT, [(0, 0), (1, 0), (2, 0), (0, 1)])
    if frac is None:
        starts = [(a, 0) for a in range(N + 1)] + [(0, b) for b in range(1, 4)]
    else:
        lo, hi = int(0.8 * frac * N), int(1.2 * frac * N)
        step = max(1, (hi - lo) // 400)
        starts = [(a, 0) for a in range(lo, hi, step)]
    lC, chC, symC = best(N, sC, starts)
    if frac is not None:
        a0 = chC[0]; step = max(1, (int(1.2 * frac * N) - int(0.8 * frac * N)) // 400)
        lC, chC, symC = best(N, sC, [(a, 0) for a in range(max(0, a0 - step), a0 + step + 1)])
    tgt = 2 * math.sqrt(k)
    r = math.sqrt(lC / lT); rs = math.sqrt(symC / symT)
    return (tgt - r) * N, (tgt - rs) * N, chC[0] / N


def main():
    print("=" * 100)
    print("Computation 206: the two-monomial deficit to N = 15360, no plateau")
    print("=" * 100)
    ok_all = True
    rows = {}
    for k in (2, 6):
        print(f"  k = {k}:")
        frac = None; prev = None; incs = []
        for N in (480, 960, 1920, 3840, 7680, 15360):
            if N > NMAX:
                break
            g, gs, fr = deficit(k, N, frac if N > 3840 else None)
            frac = fr
            inc = None if prev is None else g - prev
            if inc is not None:
                incs.append(inc)
            prev = g
            rows[(k, N)] = (g, gs, fr)
            print(f"      N={N:6d}: deficit x N = {g:+.4f}   symbol {gs:+.4f}   chain start/N = {fr:.4f}" + (f"   increment {inc:+.4f}" if inc is not None else ""))
        ok_all &= all(i > 0 for i in incs) and (rows[(k, 15360)][0] - rows[(k, 480)][0]) > 0.1
        print(f"      rises monotonically by {rows[(k, 15360)][0] - rows[(k, 480)][0]:.3f} over five doublings; no plateau")
    # the symbol constant moves the other way
    ok_sym = rows[(2, 3840)][1] < rows[(2, 480)][1] and rows[(2, 3840)][0] > rows[(2, 3840)][1]
    print(f"  the symbol value falls while the exact value rises; they cross for k = 2 between N = 480 and 3840: {ok_sym}")
    # the sharp value stands
    ok_sharp = max(v[0] for v in rows.values()) < 3.6736
    print(f"  every value here is below the certified sharp constant 3.6736 at (6, 12): {ok_sharp}")
    # tamper: the single-shift maximiser's chain is not the bridge's
    m, al, sT, sC = specs(2); N = 3840; s = N - 2 * m
    a_star = ((m + 1) * s - (m - 1)) // (2 * m) + 1; b_star = s - a_star
    lT, _, _ = best(N, sT, [(0, 0)]); lC_scan, chC, _ = best(N, sC, [(a, 0) for a in range(0, N + 1, 4)])
    lC_single, _, _ = best(N, sC, [(a_star - b_star, 0)])
    ok_t = lC_single < lC_scan
    print(f"  tamper, k = 2, N = {N}: commutator norm on the single-shift maximiser's chain {math.sqrt(lC_single):.6f} against the scanned chain {math.sqrt(lC_scan):.6f} (start {chC[0]} versus {a_star - b_star}): the drift is real: {ok_t}")
    print()
    ok = ok_all and ok_sym and ok_sharp and ok_t
    print("SUMMARY")
    print(f"  [{'PASS' if ok_all else 'FAIL'}] the deficit times N rises monotonically through N = 15360 for k = 2 and k = 6")
    print(f"  [{'PASS' if ok_sym else 'FAIL'}] the symbol value moves the other way; the symbol route to a closed form is closed")
    print(f"  [{'PASS' if ok_sharp else 'FAIL'}] the certified sharp value 3.6736 stands as the maximum")
    print(f"  [{'PASS' if ok_t else 'FAIL'}] the extremal chain's drift is real")
    print()
    print("  CONCLUSION: the per-k constant comps 157 and 203 reported is a value at N of a few hundred,")
    print("  not a limit; to N = 15360 the deficit times N rises slowly and the data do not decide between")
    print("  a slowly approached limit and a logarithm.  The sharp value and the uniform bound stand on")
    print("  every point computed; their status as N -> infinity for fixed k is what remains open.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()