#!/usr/bin/env python3 """ PROVENANCE: NUMERICAL Computation 214 -- the sharp value 3.6736 in closed form: an algebraic number at the operative corner (6, 12), and why that corner is the maximum ============================================================================== WHAT THIS DOES. The uniform bound of eq. gap-uniform-bound rests on a sharp value, G(6, 12) = 3.6736, that comp 65 found by sweeping, comp 150 confirmed, comp 157 re-checked to k = 80 and comp 203 certified in exact arithmetic to 1e-7. The book called it a certified computation rather than a derivation, and comp 213, having found the large-N limit in closed form, left it as the one quantitative thread: the value sits in the small-headroom regime N close to 2m that the asymptotic laws do not cover. This file derives it there. THE SMALL-HEADROOM REGIME. Write N = 2m + h. The lower shift (m, m) of the bridge acts on a + b <= h and the upper shift (m+1, m+1) on a + b <= h - 2, so every chain has at most floor(h/2) + 1 points, and G(k, 2m + h) is an algebraic function of (m, alpha, h): a square root of a quadratic for h = 2, a root of a cubic for h = 4. At the operative corner N = 2k the headroom is h = 2k - 2m(k), which grows with k (m ~ 2 sqrt k); the only k <= 80 with h <= 4 are k = 2, 3, 4 (h = 2) and k = 5, 6 (h = 4). THE VALUE. At (k, N) = (6, 12): m = 4, h = 4, and both extremal chains are the diagonal (0,0), (1,1), (2,2) (the commutator's chain is the diagonal too; its drift toward 0.43 N, comp 206, is a large-N effect). The Bergman weights are exact rationals, and the two 3 x 3 tridiagonal matrices have the characteristic polynomials T: 122367445200 L^3 - (24820950 a^2 + 258112008) L^2 + (1050 a^4 + 13419 a^2 + 164108) L - 30 = 0, C: 2913510600 L^3 - (17382750 a^2 + 117070800) L^2 + (21875 a^4 + 182000 a^2 + 1432256) L - 5120 = 0, with a = alpha(6), and G(6, 12) = 12 (2 sqrt 6 - sqrt(L_C / L_T)) with L_T, L_C their largest roots. alpha(6) is itself algebraic: Lemma 5(c)'s closure R_4(q) = 2 sqrt 6 squares to a polynomial system in q and sqrt(1 - q^2), and eliminating gives a degree-24 polynomial with integer coefficients of which alpha(6) = 1.17861153449039994413... is the root in (1, 2) (its mirror image under alpha -> -alpha is the other degree-24 factor). So G(6, 12) = 3.673556413812905434109146508787..., an explicit algebraic number, inside comp 203's enclosure and equal to the float evaluator to 1e-13. WHY (6, 12). At fixed (m, h) the function alpha -> G is single-humped, with the hump where the two monomials are comparable (alpha near one) and its top rising with m: at h = 4 the tops are 2.00, 3.30, 3.81, 4.13, 4.45 for m = 2 to 6, at h = 2 they are 2.42, 2.92, 3.41, 3.90, 4.39 (beyond alpha = 20 the tail rises again, slowly, toward the single-monomial value of w^(m+1)). The integer lattice (k, m(k), alpha(k)) samples those curves at few points, and k = 6 is the only one with h <= 4 and alpha on the hump: alpha(6) = 1.18 sits within 4 percent of the (m, h) = (4, 4) top, while k = 5 shares that headroom with alpha = 28.7, far down the tail (G = 1.08), and k = 2, 3, 4 have h = 2 with lower humps. Larger m would give higher humps but larger m means larger k and so larger headroom, off the small-h curves; the per-k maxima then fall toward the large-N limits of comp 213. The sharp value is therefore not a bound that some analytic mechanism sets but the one point where a small headroom meets a balanced bridge. THE SCAN. Over every k <= 80 and N from the smallest operative value to 600, scanning every chain, the maximum of |G| is 3.6736 at (6, 12); next are 3.084 at (4, 10), 2.862 at (9, 18), 2.684 at (24, 48); for k >= 7 the per-k maximum sits a few steps above the corner and falls with k. In this file k <= 40 and N <= 200 are rescanned. With comp 213's limits (below 2.47 for k <= 80) and its data to N = 61440, the supremum over the operative regime for k <= 80 is 3.6736, derived at its maximiser and scanned elsewhere. Tampers: (i) evaluating the commutator on the chain through (2, 0) instead of the diagonal gives 5.58, so the chain identification is load-bearing; (ii) the other real closure root of the degree-24 polynomial, 29.85, is comp 59's other branch and gives G = -5.8, so the branch alpha in (1/(m+1), 1/m) in q is the one the bound is about; (iii) alpha moved by 1e-3 moves G by 3e-3, so the fifteen-digit alpha of comp 203 fixes G to 1e-13, which is the agreement found. Runs in about a minute; imports comp 206 (and through it comp 59) beside it; sympy and mpmath for the exact part. """ import math import numpy as np import sympy as sp import mpmath as mpm import computation_206_se as c206 DPS = 50 # ------------------------------------------------------------- exact weights def w2_exact(a, b, p, q, N): if a + b + p + q > N: return sp.Integer(0) return sp.Rational(sp.factorial(a + p) * sp.factorial(b + q) * sp.factorial(a + b + 2), sp.factorial(a) * sp.factorial(b) * sp.factorial(a + b + p + q + 2)) def chain_matrix(spec, a0, b0, N, m, al): """Exact tridiagonal M^T M of the chain (a0 + t, b0 + t) at truncation N, alpha symbolic.""" c1, p1, q1, c2, p2, q2 = spec pts = []; t = 0 while a0 + b0 + 2 * t + 2 * m <= N: pts.append((a0 + t, b0 + t)); t += 1 u = [c1 * sp.sqrt(w2_exact(a, b, p1, q1, N)) for a, b in pts] v = [c2 * sp.sqrt(w2_exact(a, b, p2, q2, N)) for a, b in pts] n = len(pts); M = sp.zeros(n, n) for i in range(n): M[i, i] = u[i] ** 2 + v[i] ** 2 if i + 1 < n: M[i, i + 1] = M[i + 1, i] = v[i] * u[i + 1] return pts, M def top_eig(M, al_sym, al_val): Mn = mpm.matrix([[mpm.mpf(str(sp.N(M[i, j].subs(al_sym, al_val), DPS))) for j in range(M.shape[1])] for i in range(M.shape[0])]) return max(mpm.eigsy(Mn)[0]) # ------------------------------------------------- alpha(6) as an algebraic number def alpha6_polynomial(): """Degree-24 integer polynomial with alpha(6) as a root, from Lemma 5(c)'s closure at m = 4.""" q, s, a = sp.symbols('q s alpha') num = 20 * q * (1 - q) * (1 + q) ** 2 * (1 - q ** 2) # R_4(q) numerator, (1-q)^{3/2}(1+q)^{5/2} s = (1-q)(1+q)^2 s^2 den = -16 * q + 5 * s * (5 * q - 1) + 4 P1 = sp.expand(num ** 2 - 24 * den ** 2) # R^2 = 4k = 24 P2 = s ** 2 + q ** 2 - 1 Rq = sp.resultant(P1, P2, s) Pa = sp.expand(a * 5 * s * (5 * q - 1) - 8 * (1 - 4 * q)) # alpha = 8(1-4q)/(5 s (5q-1)) R1 = sp.resultant(Pa, P2, s) R2 = sp.resultant(R1, Rq, q) facs = [f for f, _ in sp.factor_list(R2, a)[1] if sp.degree(f, a) > 0] al_num = mpm.mpf(c206.c59.alpha_of_k(6)) for f in facs: coeffs = [int(x) for x in sp.Poly(f, a).all_coeffs()] val = mpm.polyval([mpm.mpf(x) for x in coeffs], al_num) scale = max(abs(mpm.mpf(x)) for x in coeffs) if abs(val) / scale < 1e-10: return coeffs raise RuntimeError("no factor vanishes at alpha(6)") def refine_root(coeffs, x0): roots = mpm.polyroots([mpm.mpf(x) for x in coeffs], maxsteps=300, extraprec=300) return min((r for r in roots if abs(r.imag) < mpm.mpf(10) ** (-30)), key=lambda r: abs(r.real - x0)).real # --------------------------------------------------------------- float side _WCACHE = {} def _lattice_weights(N, spec): """comp 206's chain weights for every lattice point at once: W[a, b] = c w(p, q; a, b), masked at a + b + p + q <= N, so a chain (a0 + t, b0 + t) is a diagonal slice. The formula is comp 206's (gammaln), so the numbers are the same to the last bit; the cache makes the scan over every chain start cost one slice per chain.""" key = (N, spec) if key not in _WCACHE: from scipy.special import gammaln c1, p1, q1, c2, p2, q2 = spec a = np.arange(N + 1)[:, None]; b = np.arange(N + 1)[None, :] 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) _WCACHE[key] = (c1 * w(p1, q1), c2 * w(p2, q2)) return _WCACHE[key] def best_cached(N, spec, starts): """comp 206's best() on the cached lattice weights: the top eigenvalue over the chains listed in starts, with the same chain lengths and the same tridiagonal as comp 206.""" U, V = _lattice_weights(N, spec) b = (0.0, None, 0.0) for a0, b0 in starts: tmax = (N - a0 - b0) // 2 if tmax < 0: continue t = np.arange(tmax + 1) u = U[a0 + t, b0 + t]; v = V[a0 + t, b0 + t] if len(u) == 0 or u.max() == 0: continue l, sy = c206.lam_sym(u, v) if l > b[0]: b = (l, (a0, b0), sy) return b def G_float(k, N, m=None, al=None, starts=None): if m is None: m, al, sT, sC = c206.specs(k) else: sT = (1.0, m, m, al, m + 1, m + 1); sC = (float(m), m + 1, m - 1, al * (m + 1), m + 2, m) tgt = 2 * math.sqrt(k) st = starts or ([(a, 0) for a in range(N + 1)] + [(0, b) for b in range(1, N + 1)]) lT, chT, _ = best_cached(N, sT, [(0, 0), (1, 0), (2, 0), (0, 1)]) lC, chC, _ = best_cached(N, sC, st) return (tgt - math.sqrt(lC / lT)) * N, chT, chC def k_closed_by(m, al): """The k whose bridge (m, alpha) closes in the limit: 4k = S_C,max / S_T.""" r = np.linspace(0.5, 0.999, 100001) SC = (r ** (m + 1) * (1 - r) ** (m - 1) * (m + al * (m + 1) * np.sqrt(r * (1 - r))) ** 2).max() return SC / (0.25 ** m * (1 + al / 2) ** 2) / 4 def main(): mpm.mp.dps = DPS print("=" * 100) print("Computation 214: the sharp value 3.6736 in closed form, and why (6, 12) is the maximum") print("=" * 100) checks = {} k, N = 6, 12 m, al_f, sT_f, sC_f = c206.specs(k) a = sp.symbols('alpha') sT = (1, m, m, a, m + 1, m + 1); sC = (m, m + 1, m - 1, a * (m + 1), m + 2, m) # ---- Part 1: alpha(6) as an algebraic number print("PART 1. alpha(6) is algebraic of degree 24") coeffs = alpha6_polynomial() al6 = refine_root(coeffs, mpm.mpf(al_f)) print(f" integer polynomial of degree {len(coeffs) - 1}: leading {coeffs[0]}, constant {coeffs[-1]}") print(f" alpha(6) = {mpm.nstr(al6, 40)} (comp 59: {al_f:.15f})") checks['alpha'] = abs(al6 - mpm.mpf(al_f)) < 1e-14 al_val = sp.Float(mpm.nstr(al6, DPS), DPS) # ---- Part 2: the chains and the cubics at (6, 12) print("PART 2. the two extremal chains at (6, 12) and their characteristic polynomials") lam, best = {}, {} L = sp.symbols('L') for name, spec in (("T", sT), ("C", sC)): starts = [(x, 0) for x in range(0, N - 2 * m + 1)] + [(0, y) for y in range(1, N - 2 * m + 1)] allc = [] for a0, b0 in starts: pts, M = chain_matrix(spec, a0, b0, N, m, a) if not pts: continue allc.append((top_eig(M, a, al_val), pts, M)) allc.sort(key=lambda x: -x[0]) lam[name], pts, M = allc[0]; best[name] = (pts, M, allc) cp = sp.Poly(sp.expand(M.charpoly(L).as_expr() * sp.lcm([sp.denom(x) for x in M])), L) print(f" {name}: extremal chain {pts}; top eigenvalue {mpm.nstr(lam[name], 30)}") print(f" exact weights on it: diagonal {[sp.nsimplify(M[i, i]) for i in range(M.shape[0])]}, off-diagonal {[sp.nsimplify(M[i, i + 1]) for i in range(M.shape[0] - 1)]}") print(f" characteristic polynomial: {sp.factor(cp.as_expr())} = 0") print(f" other chains' top eigenvalues: {', '.join(f'{c[1][0]}: {mpm.nstr(c[0], 6)}' for c in allc[1:5])}") checks['diagonal'] = best['T'][0] == [(0, 0), (1, 1), (2, 2)] and best['C'][0] == [(0, 0), (1, 1), (2, 2)] G_exact = (2 * mpm.sqrt(6) - mpm.sqrt(lam['C'] / lam['T'])) * 12 g_float, chT, chC = G_float(k, N) print(f" G(6, 12) = 12 (2 sqrt 6 - sqrt(L_C / L_T)) = {mpm.nstr(G_exact, 40)}") print(f" float evaluator (comps 157, 206): {g_float:.15f}; comp 203's enclosure [3.67356, 3.67356]") checks['value'] = abs(G_exact - g_float) < 1e-12 and abs(G_exact - mpm.mpf('3.67356')) < 6e-6 # exact weights against the float chain u, v = c206.chain_vec(N, sC_f, 0, 0) d_f = u * u + v * v; d_e = [float(sp.N(best['C'][1][i, i].subs(a, al_val), 20)) for i in range(3)] checks['weights'] = max(abs(d_f[i] - d_e[i]) / d_e[i] for i in range(3)) < 1e-12 print(f" exact diagonal weights of the C chain against comp 206's floats: {max(abs(d_f[i] - d_e[i]) / d_e[i] for i in range(3)):.1e} relative") print() # ---- Part 3: the small-headroom regime and the hump print("PART 3. headroom at the operative corner, and the hump in alpha at fixed (m, h)") hs = {} for kk in range(2, 81): mm = c206.c59.m_of_k(kk); Nmin = max(2 * kk, 2 * mm + 2); hs[kk] = Nmin - 2 * mm small = [kk for kk in hs if hs[kk] <= 4] print(f" headroom h = N_min - 2m(k) at the corner: k = 2..12: {[hs[kk] for kk in range(2, 13)]}; k <= 80 with h <= 4: {small}; h(80) = {hs[80]}") checks['headroom'] = small == [2, 3, 4, 5, 6] and all(hs[kk] >= 6 for kk in range(7, 81)) als = np.concatenate([np.linspace(0.2, 3.0, 57), np.array([4, 5, 7, 10, 20, 30])]) tops = {} for h in (2, 4): line = [] for mm in (2, 3, 4, 5, 6): vals = [G_float(k_closed_by(mm, x), 2 * mm + h, mm, x)[0] for x in als] i = int(np.argmax(vals)); tops[(mm, h)] = (vals[i], als[i]) # unimodal: rises to the top then falls d = np.diff(vals); j10 = int(np.sum(als <= 10)) - 1 # unimodal to alpha = 10; the far tail rises slowly toward the single-monomial value uni = all(x >= -1e-9 for x in d[:i]) and all(x <= 1e-9 for x in d[i:j10]) line.append(f"m={mm}: top {vals[i]:.2f} at alpha {als[i]:.2f}{'' if uni else ' (not unimodal)'}") print(f" h = {h}: " + "; ".join(line)) g6 = G_float(6, 12)[0]; g5 = G_float(5, 10)[0] print(f" k = 6 sits on the (m, h) = (4, 4) curve at alpha = {al_f:.3f}: G = {g6:.4f}, {100 * (1 - g6 / tops[(4, 4)][0]):.1f} percent below its top;" f" k = 5 on the (3, 4) curve at alpha = {c206.c59.alpha_of_k(5):.1f}: G = {g5:.3f} against a top of {tops[(3, 4)][0]:.2f}") checks['hump'] = g6 > 0.9 * tops[(4, 4)][0] and g5 < 0.5 * tops[(3, 4)][0] and tops[(4, 4)][0] > tops[(3, 4)][0] > tops[(2, 4)][0] print() # ---- Part 4: the scan print("PART 4. the operative regime rescanned, k <= 40, N <= 200 (offline: k <= 80, N <= 600, same maximum)") res = [] for kk in range(2, 41): mm, aa, sT_k, sC_k = c206.specs(kk); Nmin = max(2 * kk, 2 * mm + 2); tgt = 2 * math.sqrt(kk) chC = None; per = (0, None) for NN in list(range(Nmin, Nmin + 31)) + list(range(Nmin + 32, 201, 2)): lT = best_cached(NN, sT_k, [(0, 0), (1, 0), (2, 0), (0, 1)])[0] if NN <= Nmin + 30: st = [(x, 0) for x in range(NN + 1)] + [(0, y) for y in range(1, NN + 1)] else: a0 = chC[0]; w = max(10, NN // 6); st = [(x, 0) for x in range(max(0, a0 - w), min(NN, a0 + w) + 1)] + [(0, y) for y in range(1, 4)] lC, chC, _ = best_cached(NN, sC_k, st) g = abs(tgt - math.sqrt(lC / lT)) * NN if g > per[0]: per = (g, NN) res.append((per[0], kk, per[1])) res.sort(reverse=True) print(" per-k maxima, largest first: " + ", ".join(f"{g:.4f} at ({kk}, {NN})" for g, kk, NN in res[:6])) checks['scan'] = res[0][1:] == (6, 12) and abs(res[0][0] - 3.6736) < 1e-4 and res[1][0] < 3.1 print() # ---- Tampers print("TAMPERS") g_off = (2 * math.sqrt(6) - math.sqrt(float(best['C'][2][2][0] if best['C'][2][2][1] == [(2, 0), (3, 1)] else [c[0] for c in best['C'][2] if c[1] == [(2, 0), (3, 1)]][0]) / float(lam['T']))) * 12 t1 = g_off > 5 print(f" (i) commutator on the chain through (2, 0): G = {g_off:.3f} against {float(G_exact):.4f}: rejected: {t1}") other = refine_root(coeffs, mpm.mpf(29.85)) g_other = G_float(6, 12, 4, float(other))[0] t2 = abs(g_other - float(G_exact)) > 1 print(f" (ii) the other real closure root, alpha = {mpm.nstr(other, 8)}: G = {g_other:.3f}: a different branch: {t2}") g_pert = G_float(6, 12, 4, al_f + 1e-3)[0] t3 = 1e-3 < abs(g_pert - g_float) < 1e-1 print(f" (iii) alpha moved by 1e-3: G moves by {g_pert - g_float:+.4e}, so fifteen digits of alpha fix thirteen of G: {t3}") checks['tampers'] = t1 and t2 and t3 print() ok = all(checks.values()) print("SUMMARY") print(f" [{'PASS' if checks['alpha'] else 'FAIL'}] alpha(6) is the root in (1, 2) of an explicit degree-24 integer polynomial") print(f" [{'PASS' if checks['diagonal'] and checks['weights'] else 'FAIL'}] at (6, 12) both extremal chains are the diagonal three-point chain, with exact rational weights") print(f" [{'PASS' if checks['value'] else 'FAIL'}] G(6, 12) = 12 (2 sqrt 6 - sqrt(L_C/L_T)) from two explicit cubics, inside comp 203's enclosure, equal to the float to 1e-12") print(f" [{'PASS' if checks['headroom'] else 'FAIL'}] the headroom at the operative corner is at most 4 only for k <= 6 and grows with k") print(f" [{'PASS' if checks['hump'] else 'FAIL'}] G is single-humped in alpha at fixed (m, h); k = 6 is the one small-headroom point on a hump") print(f" [{'PASS' if checks['scan'] else 'FAIL'}] the rescanned operative regime has its maximum at (6, 12) and nothing else above 3.1") print(f" [{'PASS' if checks['tampers'] else 'FAIL'}] tampers rejected") print() print(" CONCLUSION: the sharp value is derived: an explicit algebraic number, 3.6735564138129054...,") print(" the root of two cubics in alpha(6) at the diagonal three-point chain of the corner (6, 12),") print(" and the corner is the maximum because it is the only point of the operative regime where a") print(" headroom of four meets a balanced bridge. Elsewhere the regime is scanned (k <= 80, N <= 600)") print(" and bounded by comp 213's limits; no analytic mechanism sets the value, the lattice does.") print(" ALL CHECKS PASSED" if ok else " SOME CHECKS FAILED") if __name__ == "__main__": main()