"""bandheight.py -- does the Gaussian limit of a class's least height predict the last failing depth J_hi(r) of each rank's band, and its growth? QUESTION. J_hi = 30, 33, 58, 60, 88, 106, 122, 152, 168, 195 at ranks 5..14 (gap.py, bandgrowth.py) is the last depth where some vector off the champion's line beats the champion's height H_ch (J_ch = J_hi + 1 at every rank read). bandscale.py's limit of the dual bound gives J_c; does the same limit of the class's own least height give J_hi? THE LIMIT (by hand, from bandscale.py's). Write the cofactor of class s as (1 - x)^s c', c' = sum_k e_k (1 - x)^k, e_0 = c'(1) = 1. The coefficients of (1 + x)^J (1 - x)^n tend to 2^J sigma^(-1-n) He_n phi / sqrt(2 pi), sigma = sqrt(J)/2, phi = e^(-y^2/2), so b_k = e_k sigma^(-k) and the class's real least height over H_ch tends to sigma^d mu_s / M_(r-1), d = r - 1 - s, mu_s = min over real b of max |(He_s + sum_(m<=d) b_m He_(s+m)) phi|, M_n = max |He_n phi|. Class s beats the champion below J_hi*_s = 4 (M_(r-1) / mu_s)^(2/d); J_hi* is the largest over s <= r - 2. Parity (proved): a term He_(s+m) with m odd has the parity opposite to He_s, and max(|E + O|, |E - O|) = |E| + |O| >= |E| at y and -y, so mu_s is a problem in the even m only; at d = 1, mu_s = M_s. Holder: kappa*_s <= mu_s, so J_hi* <= J* is forced and reads nothing. The lattice: the integer least height is at least the real one, so the exact real relaxation's last crossing J_hi^R >= J_hi is forced. THE CONSTANT (by hand). At d = 2, in the leading form of bandscale.py's THE GROWTH, f ~ A_s delta y sin(w y - s pi/2) e^(-y^2/4), both minimisers alike, the oscillation reaching its envelope: mu ~ A_s delta sqrt(2) e^(-1/2), kappa ~ A_s delta pi sqrt(2 pi)/16, so mu/kappa -> 1.7428 and J_hi* ~ 0.5738 J* ~ 4.663 r^(3/2). At d = 4 the envelope He_2(y) e^(-y^2/4) gives a constant near 3.73, so d = 2 binds asymptotically. TRANSPLANT: the oscillation reaching its envelope needs sqrt(s) large; at s ~ 5 it is not. THE ARMS. (G) mu_s by an exchange LP on a fine grid, the even m only, every s <= r - 2, ranks 3..40: J_hi*_s, J_hi*, the binding order. (R) the exact real relaxation at finite J: per class an LP over e_1..e_d on the J + r coefficients of the columns (1 + x)^J (1 - x)^(s+k), exact in integers then scaled; J_hi^R the first depth, scanning down from J_c - 1, where some class beats H_ch, its beat rechecked exactly in rationals at the LP's e. Ranks 5..14. (A) mu_s at d = 1..4 to rank 801 from He_n/sqrt(n!), factorials in logs, by golden-section search on the even parameters. PREDICTIONS (frozen before the engine). P1 J_hi*/J_hi in [0.9, 1.2] at every rank 5..14 (TRANSPLANT from J*/J_c = 1.019..1.059; the lattice pushes J_hi below the real crossing). P2 J_hi^R/J_hi in [1.0, 1.15] and J_hi*/J_hi^R in [0.9, 1.1], ranks 5..14. P3 the binding order of J_hi* is d = 2 (s = r - 3) at ranks 5..14. P4 mu/kappa at d = 2 within 3% of 1.7428 at r = 800; J_hi*(800)/800^1.5 within 5% of 4.663. KILLS. K1 J_hi*/J_hi outside [0.8, 1.25] at some rank 5..14. K2 J_hi*(800)/800^1.5 outside [4.2, 5.1], or log2(J_hi*(800)/J_hi*(400)) outside [1.45, 1.55]. P1..P4 read, not gated. CONTROLS, read first. C1 the LP over every m (odd included) at d <= 4 equals the even-only LP to relative 1e-6, ranks 3..40 (at d = 1 that is mu_s = M_s). C2 at rank 5, J = 30, 36, ..., 60, the (R) LP's least height at d = 1, 2 against a golden-section search (nested at d = 2: the parity is the limit's, not exact at finite J, so both parameters are searched), relative 1e-6. C3 (R) over its limit at J = 16 r^2, s = r - 3, ranks 5..14, within 2% of 1. C4 (A) against (G) at r = 5..40, d = 1..4, relative 1e-4. C5 J_hi^R >= J_hi at every rank. C6 mu_s >= kappa*_s at every s read. INSTRUMENT, before the counted run. A first run was stopped past twice its estimate in (G), and the grid cut from 120001 to 30001 points (the cost is the LP, not the grid). A second read C1 4.81e-06 and C6 False: kappa by the normal equations loses the projection at large d, now by least squares; HiGHS's absolute tolerances on rows of size ~1e-2, now each column scaled to max 1 (HiGHS here returned no solution at tolerances 1e-9 and 1e-10). Its J_hi*, (R) and (A) rows equal the counted run's; (G) took 227.7 s of the counted run. FINDINGS (copied from the printed output). (1) Controls: C1 worst 2.56e-07; C2 1.11e-15; C3 1.0009..1.0035 at J = 16 r^2, ranks 5..14; C4 7.64e-07; C5 J_hi^R >= J_hi at every rank, every beat rechecked exactly in rationals; C6 holds. (2) K1 misses, P1 holds: J_hi* = 35.2, 36.8, 63.0, 69.5, 93.5, 112.4, 131.3, 161.6, 176.1, 201.2 at ranks 5..14 against J_hi = 30, 33, 58, 60, 88, 106, 122, 152, 168, 195: J_hi*/J_hi = 1.172, 1.115, 1.086, 1.158, 1.063, 1.060, 1.076, 1.063, 1.048, 1.032, falling with the rank but not monotone. The binding order is d = 2 (s = r - 3) at every rank 4..40 (P3); at rank 3 it is d = 1, J_hi* = 10.9, and 14.3 at rank 4, where no cell fails and the champion attains h from J = 7 and 13. (3) P2 holds. The exact real relaxation crosses last at J_hi^R = 34, 35, 60, 68, 88, 109, 130, 156, 170, 199, always in class s = r - 3, no near-tie: the lattice's share J_hi^R/J_hi = 1.133, 1.061, 1.034, 1.133, 1.000, 1.028, 1.066, 1.026, 1.012, 1.021, the limit's own J_hi*/J_hi^R = 1.034, 1.051, 1.049, 1.022, 1.063, 1.031, 1.010, 1.036, 1.036, 1.011. (4) K2 misses, P4 holds: J_hi*/r^1.5 = 4.3771, 4.5680, 4.6024, 4.6254, 4.6466 at r = 50, 100, 200, 400, 800 against the hand 4.663, the local exponent 400..800 1.5066, d = 2 binding among d = 1..4 at every r read; mu/kappa at d = 2 is 1.7415 at r = 800 and 1.7428 at 801 against the hand 1.7428, ragged at small r (1.40 at rank 5, 1.52 to 1.77 over ranks 6..14). RUN RECORD. 248.5 s, peak commit 139.5 MB under a 512 MB guard. Its exit code is red on C1..C6, on a beat not rechecked exactly, or on K1 or K2; P1..P4 are printed history. """ import math import os import sys import time from fractions import Fraction from math import comb os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import numpy as np # noqa: E402 from scipy.optimize import linprog # noqa: E402 import bandend # noqa: E402 JHI = {5: 30, 6: 33, 7: 58, 8: 60, 9: 88, 10: 106, 11: 122, 12: 152, 13: 168, 14: 195} # gap.py, bandgrowth.py JC = {5: 47, 6: 57, 7: 97, 8: 108, 9: 157, 10: 170, 11: 227, 12: 241, 13: 305, 14: 319} # bandend.py, bandgrowth.py YG = np.linspace(-15.0, 15.0, 30001) PHI = np.exp(-YG ** 2 / 2) W2 = PHI * PHI NMAX = 41 BIG = (50, 100, 200, 201, 400, 401, 800, 801) def htable(needed): """{n: h_n phi} on YG for n in needed, h_n = He_n / sqrt(n!).""" out, prev, cur = {}, None, np.ones_like(YG) for n in range(max(needed) + 1): if n in needed: out[n] = cur * PHI nxt = (YG * cur - (n ** 0.5 * prev if n else 0)) / (n + 1) ** 0.5 prev, cur = cur, nxt return out def cheb(F0, Fs): """min over x of max |F0 + sum x_i Fs[i]| on YG, by an exchange LP: solve on an active set, add the grid's violating local maxima, repeat. Returns (value on the full grid at x, x).""" k = len(Fs) if k == 0: return float(np.max(np.abs(F0))), np.zeros(0) act = np.unique(np.concatenate([np.arange(0, len(YG), 15), [int(np.argmax(np.abs(F0)))]])) # scaled so every column has max 1: HiGHS's tolerances are absolute s0 = float(np.max(np.abs(F0))) A = np.array(Fs).T cs = np.max(np.abs(A), axis=0) A, F0 = A / cs, F0 / s0 top, x = cheb_scaled(F0, A, act, k) return top * s0, x * s0 / cs def cheb_scaled(F0, A, act, k): for _ in range(40): Aa, f0 = A[act], F0[act] ones = np.ones((len(act), 1)) A_ub = np.vstack([np.hstack([Aa, -ones]), np.hstack([-Aa, -ones])]) b_ub = np.concatenate([-f0, f0]) c = np.zeros(k + 1) c[-1] = 1 res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=[(None, None)] * k + [(0, None)], method="highs") if res.status != 0: raise RuntimeError("exchange LP: " + res.message) x, t = res.x[:k], res.x[-1] g = np.abs(F0 + A @ x) top = float(np.max(g)) if top <= t * (1 + 1e-10): return top, x loc = np.where((g[1:-1] >= g[:-2]) & (g[1:-1] >= g[2:]) & (g[1:-1] > t * (1 + 1e-10)))[0] + 1 act = np.unique(np.concatenate([act, loc, [int(np.argmax(g))]])) return top, x def kappa_h(H, s, d): """kappa*_s in the h basis: the L2(dy) projection, int f^2 / int |f|.""" b = np.linalg.lstsq(np.array([H[s + m] for m in range(1, d + 1)]).T, -H[s], rcond=None)[0] f = H[s] + sum(bi * H[s + i] for i, bi in enumerate(b, 1)) return float(np.sum(f * f) / np.sum(np.abs(f))) def mu_h(H, s, d, odd=False): ms = [m for m in range(1, d + 1) if odd or m % 2 == 0] return cheb(H[s], [H[s + m] for m in ms])[0] def jstar_from(r, s, Mh, val): """4 (M_(r-1)/val_s)^(2/d) for He, from the h-basis values.""" d = r - 1 - s logfac = 0.5 * sum(math.log(j) for j in range(s + 1, r)) return 4 * math.exp((logfac + math.log(Mh / val)) * 2 / d) # ---------------------------------------------------------------- (R) def columns(J, r, s): """exact integer columns (1 + x)^J (1 - x)^(s+k), k = 0..d, width J + r.""" base = bandend.Q(J, s) cols, cur = [], base for k in range(r - s): cols.append(cur + [0] * (J + r - len(cur))) cur = bandend.pmul(cur, [1, -1]) return cols def real_least(J, r, s, want_x=False): """least height over real c' with c'(1) = 1 in class s, over its scale S_0 = max |column 0|; returns (t, S_0, cols, x scaled).""" cols = columns(J, r, s) S = [max(abs(v) for v in c) for c in cols] F = [np.array([v / S[k] for v in c]) for k, c in enumerate(cols)] k = len(F) - 1 if k == 0: return 1.0, S[0], cols, np.zeros(0) A = np.array(F[1:]).T n = A.shape[0] ones = np.ones((n, 1)) A_ub = np.vstack([np.hstack([A, -ones]), np.hstack([-A, -ones])]) b_ub = np.concatenate([-F[0], F[0]]) c = np.zeros(k + 1) c[-1] = 1 res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=[(None, None)] * k + [(0, None)], method="highs") return float(res.x[-1]), S[0], cols, res.x[:k] def exact_height(cols, S0, x): """height of col_0 + sum (x_k S_0 / S_k) col_k, in rationals.""" S = [max(abs(v) for v in c) for c in cols] coef = [Fraction(1)] + [Fraction(float(xk)) * S0 / S[k + 1] for k, xk in enumerate(x)] return max(abs(sum(cf * c[i] for cf, c in zip(coef, cols))) for i in range(len(cols[0]))) def golden(f, lo, hi, it=120): g = (5 ** 0.5 - 1) / 2 a, b = lo, hi c, d = b - g * (b - a), a + g * (b - a) fc, fd = f(c), f(d) for _ in range(it): if fc < fd: b, d, fd = d, c, fc c = b - g * (b - a) fc = f(c) else: a, c, fc = c, d, fd d = a + g * (b - a) fd = f(d) return min(fc, fd) def main(): t0 = time.time() H = htable(set(range(NMAX + 1))) # C1 and (G) print("(G) mu_s by exchange LP, ranks 3..40", flush=True) c1w, c6 = 0.0, True G = {} for r in range(3, 41): Mh = float(np.max(np.abs(H[r - 1]))) js, ratio2 = {}, None for s in range(r - 1): d = r - 1 - s mu = mu_h(H, s, d) if d <= 4: full = mu_h(H, s, d, odd=True) c1w = max(c1w, abs(full / mu - 1)) kap = kappa_h(H, s, d) c6 &= mu >= kap * (1 - 1e-9) if d == 2: ratio2 = mu / kap js[s] = jstar_from(r, s, Mh, mu) top = max(js, key=js.get) G[r] = (js[top], r - 1 - top, js, ratio2) row = " r = %2d J_hi* = %8.1f at d = %d; mu/kappa (d = 2) %.4f;" \ " J_hi*/r^1.5 = %.3f [%.1f s]" % ( r, js[top], r - 1 - top, ratio2, js[top] / r ** 1.5, time.time() - t0) if r in JHI: row += "; J_hi = %d, J_hi*/J_hi = %.3f" % (JHI[r], js[top] / JHI[r]) print(row, flush=True) c1 = c1w < 1e-6 print("C1 every-m LP against the even-m LP, d <= 4: worst relative " "%.2e (%s)" % (c1w, c1)) print("C6 mu_s >= kappa*_s at every s: %s [%.1f s]" % (c6, time.time() - t0), flush=True) # C2 worst = 0.0 r = 5 for J in range(30, 61, 6): for s in (3, 2): t, S0, cols, _ = real_least(J, r, s) Fl = [np.array([v / max(abs(u) for u in c) for v in c]) for c in cols] if s == 3: g = golden(lambda a: float(np.max(np.abs(Fl[0] + a * Fl[1]))), -50, 50) else: g = golden(lambda b: golden( lambda a: float(np.max(np.abs(Fl[0] + a * Fl[1] + b * Fl[2]))), -50, 50, 80), -50, 50, 80) worst = max(worst, abs(g / t - 1)) c2 = worst < 1e-6 print("C2 the (R) LP against golden search, rank 5, J = 30..60, d = 1, 2:" " worst relative %.2e (%s) [%.1f s]" % (worst, c2, time.time() - t0), flush=True) # C3 c3 = True row = "C3 (R) over its limit at J = 16 r^2, s = r - 3:" for r in range(5, 15): J, s = 16 * r * r, r - 3 t, S0, _, _ = real_least(J, r, s) ex = t * (S0 / bandend.champ_height(r, J)) Mh = float(np.max(np.abs(H[r - 1]))) mu = mu_h(H, s, 2) lim = (J ** 0.5 / 2) ** 2 * math.sqrt(math.factorial(s) / math.factorial(r - 1)) * mu / Mh c3 &= abs(ex / lim - 1) < 0.02 row += " %d: %.4f" % (r, ex / lim) print(row + " (%s) [%.1f s]" % (c3, time.time() - t0), flush=True) # (R) print("(R) the exact real relaxation, scanning down from J_c - 1") R, c5, exact_ok = {}, True, True for r in range(5, 15): H_at = {} found = None for J in range(JC[r] - 1, 1, -1): Hc = bandend.champ_height(r, J) best = None for s in range(r - 1): t, S0, cols, x = real_least(J, r, s) q = t * (S0 / Hc) if best is None or q < best[0]: best = (q, s, cols, S0, x) H_at[J] = best[0] if best[0] < 1 - 1e-6: found = J hx = exact_height(best[2], best[3], best[4]) exact_ok &= hx < Hc break near = [J for J, q in H_at.items() if abs(q - 1) < 1e-5] R[r] = found c5 &= found >= JHI[r] print(" r = %2d J_hi^R = %d (class s = r - %d, least/H_ch %.5f)," " J_hi = %d, J_hi^R/J_hi = %.3f, J_hi*/J_hi^R = %.3f;" " near-ties %s [%.1f s]" % (r, found, r - best[1], best[0], JHI[r], found / JHI[r], G[r][0] / found, near, time.time() - t0), flush=True) print("C5 J_hi^R >= J_hi at every rank: %s; every beat rechecked exactly:" " %s" % (c5, exact_ok)) # (A) need = set() for r in list(range(5, 41)) + list(BIG): need |= set(range(max(0, r - 5), r)) HA = htable(need) def mus(r, Hx): out = {} for d in range(1, 5): s = r - 1 - d F0 = Hx[s] if d == 1: out[d] = float(np.max(np.abs(F0))) elif d in (2, 3): F2 = Hx[s + 2] out[d] = golden(lambda b: float(np.max(np.abs(F0 + b * F2))), -5, 5) else: F2, F4 = Hx[s + 2], Hx[s + 4] out[d] = golden(lambda b4: golden( lambda b2: float(np.max(np.abs(F0 + b2 * F2 + b4 * F4))), -5, 5, 70), -5, 5, 70) return out worst = 0.0 for r in range(5, 41): m = mus(r, HA) Mh = float(np.max(np.abs(HA[r - 1]))) for d, v in m.items(): a = jstar_from(r, r - 1 - d, Mh, v) worst = max(worst, abs(a / G[r][2][r - 1 - d] - 1)) c4 = worst < 1e-4 print("(A) C4 against (G) at r = 5..40, d = 1..4: worst relative %.2e" " (%s) [%.1f s]" % (worst, c4, time.time() - t0), flush=True) top = {} for r in BIG: m = mus(r, HA) Mh = float(np.max(np.abs(HA[r - 1]))) js = {d: jstar_from(r, r - 1 - d, Mh, v) for d, v in m.items()} s = r - 3 kap = kappa_h(HA, s, 2) bind = max(js, key=js.get) top[r] = js[bind] print(" r = %3d J_hi*_d, d = 1..4: %s; binding d = %d;" " J_hi*/r^1.5 = %.4f; mu/kappa (d = 2) = %.4f [%.1f s]" % (r, ["%.1f" % js[d] for d in range(1, 5)], bind, top[r] / r ** 1.5, m[2] / kap, time.time() - t0), flush=True) # verdicts p1 = [r for r in JHI if not 0.9 <= G[r][0] / JHI[r] <= 1.2] k1 = [r for r in JHI if not 0.8 <= G[r][0] / JHI[r] <= 1.25] p2 = [r for r in JHI if not (1.0 <= R[r] / JHI[r] <= 1.15 and 0.9 <= G[r][0] / R[r] <= 1.1)] p3 = [r for r in JHI if G[r][1] != 2] expo = math.log2(top[800] / top[400]) k2 = not (4.2 <= top[800] / 800 ** 1.5 <= 5.1 and 1.45 <= expo <= 1.55) print("P1 ranks with J_hi*/J_hi outside [0.9, 1.2]: %s" % p1) print("P2 ranks off P2's bands: %s" % p2) print("P3 ranks whose binding order is not d = 2: %s" % p3) print("K1 ranks with J_hi*/J_hi outside [0.8, 1.25]: %s (K1 %s)" % (k1, "fires" if k1 else "misses")) print("K2 J_hi*(800)/800^1.5 = %.4f, local exponent 400..800 = %.4f" " (K2 %s)" % (top[800] / 800 ** 1.5, expo, "fires" if k2 else "misses")) print("%.1f s" % (time.time() - t0)) return (not (c1 and c2 and c3 and c4 and c5 and c6 and exact_ok) or bool(k1) or k2) if __name__ == "__main__": sys.exit(1 if main() else 0)