"""bandscale.py -- does the Gaussian limit of the dual bound predict the depth J_c(r) from which bandend.py's certificate holds, and its growth like the square of the rank? QUESTION. bandgrowth.py reads J_c / r^2 in 1.58..1.98 over ranks 5..14, the binding order s = r - 3. Does the fixed-s Gaussian limit of DUAL_s / H_ch predict J_c(r) and its parity split? THE LIMIT (by hand). sigma = sqrt(J)/2, phi(y) = e^(-y^2/2). The coefficients of P = (1 + x)^J (1 - x)^s c'(x), c'(1) = 1, tend to 2^J sigma^(-1-s) f(y) / sqrt(2 pi), f = (He_s + sum_(m=1..d) b_m He_(s+m)) phi, the b_m free reals, d = r - 1 - s; the champion's to 2^J sigma^-r He_(r-1) phi / sqrt(2 pi). So DUAL_s / H_ch -> G_s(J) = sigma^d kappa*_s / M_(r-1), kappa*_s = int f*^2 / int |f*| at the f* least in L2(dy), M_n = max |He_n phi|, and the order s certifies from J*_s = 4 (M_(r-1) / kappa*_s)^(2/d) on. The Gaussian J* is the largest J*_s over s <= r - 2. On the circle each class's factor against its Gaussian is about exp(-s^2/(4J)), mostly cancelling between the class and the champion: (X) reads what is left. THE ARMS. (G) kappa*_s, M_(r-1) and J*_s at every s <= r - 2, ranks 3..40: the L2 projection by Gauss-Hermite quadrature (exact on polynomials), int |f*| and the max on a grid. (X) the exact E = DUAL_(r-3) / H_ch (bandend.bounds, champ_height, in integers) over its limit G_(r-3), at J = J_c and at J = a r^2, a = 1, 2, 4, 8, 16, ranks 5..14. (A) J*_s at d = 1..4 to rank 801, from He_n/sqrt(n!) by its recurrence and grid integrals, the factorials carried in logs. THE GROWTH (by hand, after (G)). Near 0, He_n e^(-y^2/4) ~ A_n cos(w_n y - n pi/2), w_n = sqrt(n + 1/2). At d = 2 the minimiser matches amplitudes, f* ~ A_s delta y sin(w_s y - s pi/2) e^(-y^2/4), delta = w_(s+2) - w_s ~ 1/sqrt(s), so kappa* ~ A_s delta pi sqrt(2 pi)/16 and J* ~ 64 r^(3/2)/(pi sqrt(2 pi)) = 8.128 r^(3/2). At d = 2k, delta^k He_k(y) is left, about 4.97 r^(3/2) at k = 2; odd d = 2k + 1 grows as r^((3k+1)/(2k+1)). H1 J*(r)/r^(3/2) within 3% of 8.128 at r = 400 and 800. H2 the binding order among d = 1..4 is s = r - 3 at every r computed. H3 J*_(r-2)/r settles (linear). H4 the parity gap J*(r)/J*(r + 1) ((r + 1)/r)^(3/2) - 1 below 1% at r >= 200. KILL K2 J*(800)/800^(3/2) outside [7.7, 8.6], or log2(J*(800)/J*(400)) outside [1.45, 1.55]. CONTROL C4 (A) reproduces (G)'s J*_s at r = 5..40, d = 1..4, to relative 1e-6. PREDICTIONS (frozen before the engine). G1 J*_(r-3)(r) / r settles (within +-25% of its mean over r = 10..40); the Gaussian limit does not give r^2, and J*_(r-3) / J_c falls with r, below 0.5 by rank 14. G2 E / G at J = J_c is below 1 and falls with r over ranks 5..14. G3 E / G -> 1 as J grows at fixed r. KILL. K1 J*_(r-3)(r) / r^2 in [1.4, 2.2] at every rank 5..14 (the limit would explain the r^2 itself). G2, G3 read, not gated. CONTROLS, read first. C1 G_(r-2) / sqrt(J) at r = 3, 4 against 0.2218, 0.1986 (bandend.py FINDINGS (2), exact at J = 1000), within 2%. C2 the exact E / G at s = r - 2, r = 3, 4, J = 1000, within 2% of 1. C3 the quadrature's int He_m He_n e^(-y^2) against the grid sum, worst relative below 1e-8 at m, n <= 41. C5 (added after the run) the projection's residual against each free term on the nodes, over |q| |He_(s+m)|, below 1e-6 at every s of ranks 3..40: the normal equations the run solved lose the projection at large d; the projection is now by least squares on the nodes, columns normalised, and every print but C4's is unchanged. FINDINGS (copied from the printed output). (1) Controls: C3 worst relative 4.78e-12; C5 5.83e-08 (its later run, 21.8 s, 269.1 MB, C4 there 9.55e-15); C1 Gaussian 0.2216, 0.1985 sqrt(J) against 0.2218, 0.1986; C2 E/G at J = 1000 1.0009, 1.0006. (2) K1 FIRES, G1 fails: the Gaussian limit predicts J_c. Its J* = 49.2, 60.4, 100.4, 112.5, 161.8, 175.0, 232.3, 246.5, 310.8, 325.9 at ranks 5..14 against J_c = 47, 57, 97, 108, 157, 170, 227, 241, 305, 319: J*/J_c = 1.019..1.059. The binding order is s = r - 2 at ranks 3, 4 (J* = 20.4, 25.4 against J_c = 19, 26) and s = r - 3 from rank 5 on, as the exact scans found. The parity split is the limit's: J*/r^2 = 1.969, 1.677, 2.048, 1.758, ... at ranks 5, 6, 7, 8. (3) The hand growth law (J* linear in r) is wrong: J*_(r-3)/r rises from 9.84 (r = 5) to 47.01 (r = 40), and J*/r^2 falls past rank 7 within each parity (2.048, 1.997, 1.920, 1.839 at ranks 7, 9, 11, 13), to 1.175 at rank 40. The r^2 is a window reading. (4) G3 holds: E/G at J = 16 r^2 is 1.002..1.007 at ranks 5..14. G2 fails: E/G at J_c is 1.031..1.111, above 1, falling but not monotone (1.044 at r = 10, 1.053 at r = 11); at J = r^2, odd at odd r, it is 1.045..1.222. The exact bound runs a little above its limit, so J_c sits a little below J*. E zigzags with the parity of J (by hand, bandend.bounds: 1.025, 0.9901, 1.0609, 1.036 at r = 5, J = 45..48; 1.0044, 0.9955, 1.0142, 1.008 at r = 10, J = 168..171), so E/G at J_c and at r^2 mixes the branches and G3 reads even J only. At r = 5 the zigzag is H_ch's alone (its steps 1.8, 2.0, 1.8077, 2.0, 1.8148 over J = 44..49, the dual bound's without a period); at r = 10 both factors alternate. (5) (A): C4 worst relative 7.99e-15. K2 misses: J*/r^1.5 = 7.5703, 7.8477, 7.9876, 8.0577, 8.0926 at r = 50, 100, 200, 400, 800, rising to 8.128 (H1 holds), the local exponent 400..800 1.5062, J* the largest over d = 1..4. Among them the binding order is d = 2 (s = r - 3) at every r read (H2), no d > 4 read past r = 40; J*_(r-2)/r = 12.276 .. 12.945, settling (H3); the parity gaps -0.0054, -0.0027, -0.0013 at r = 200, 400, 800 (H4). At r = 800 d = 4 gives 111900.6 = 4.95 r^1.5 against the hand 4.97, and d = 3 grows by 2^1.34 from 400 to 800 against the hand r^(4/3). RUN RECORD. 18.7 s, peak commit 268.9 MB under a 512 MB guard. Its exit code is red on C1..C5, on G3 read as E/G within 1% of 1 at J = 16 r^2 (its gated check of the limit against the exact bound) or on K2; G1, G2 and K1 are printed history. """ import os import sys import time os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import numpy as np # noqa: E402 import bandend # noqa: E402 JC = {5: 47, 6: 57, 7: 97, 8: 108, 9: 157, 10: 170, 11: 227, 12: 241, 13: 305, 14: 319} # bandend.py (3), bandgrowth.py (2) NMAX = 42 XQ, WQ = np.polynomial.hermite.hermgauss(60) # weight e^(-y^2) SW = np.sqrt(WQ) YG = np.linspace(-30.0, 30.0, 240001) DY = YG[1] - YG[0] def he_table(y): H = [np.ones_like(y), y.copy()] for n in range(1, NMAX): H.append(y * H[n] - n * H[n - 1]) return H HQ = he_table(XQ) HG = he_table(YG) PHI = np.exp(-YG ** 2 / 2) def inner(m, n): return float(np.sum(WQ * HQ[m] * HQ[n])) def M(n): return float(np.max(np.abs(HG[n] * PHI))) def orth(s, d, b): """the projection's defining property: its residual against each free term, relative, worst over the terms.""" q = HQ[s] + sum(bi * HQ[s + i] for i, bi in enumerate(b, 1)) nq = inner_vec(q, q) ** 0.5 return max(abs(inner_vec(q, HQ[s + m])) / (nq * inner(s + m, s + m) ** 0.5) for m in range(1, d + 1)) def inner_vec(u, v): return float(np.sum(WQ * u * v)) def kappa(s, d): """kappa*_s with d free higher terms: the L2(phi^2) projection, by least squares on the quadrature's nodes (the normal equations lose the projection at large d).""" A = SW[:, None] * np.array(HQ[s + 1:s + d + 1]).T cn = np.linalg.norm(A, axis=0) # the He_n span 20 orders of size b = np.linalg.lstsq(A / cn, -SW * HQ[s], rcond=None)[0] / cn q = HG[s].copy() for i, bi in enumerate(b, 1): q = q + bi * HG[s + i] f = q * PHI return float(np.sum(f * f) * DY / (np.sum(np.abs(f)) * DY)), b def jstar(r): out = {} for s in range(r - 1): d = r - 1 - s k, _ = kappa(s, d) out[s] = 4 * (M(r - 1) / k) ** (2 / d) return out def gauss_ratio(r, s, J): d = r - 1 - s return (J ** 0.5 / 2) ** d * kappa(s, d)[0] / M(r - 1) def exact_ratio(r, s, J): return float(bandend.bounds(r, J, s)[1] / bandend.champ_height(r, J)) def main(): t0 = time.time() # C3 worst = 0.0 for m in range(NMAX): for n in range(m, NMAX): g = float(np.sum(HG[m] * HG[n] * PHI * PHI) * DY) q = inner(m, n) scale = (inner(m, m) * inner(n, n)) ** 0.5 worst = max(worst, abs(g - q) / scale) c3 = worst < 1e-8 print("C3 quadrature against the grid, worst relative %.2e (%s)" % (worst, c3)) # C1, C2 c1 = c2 = True for r, lit in ((3, 0.2218), (4, 0.1986)): g = gauss_ratio(r, r - 2, 1000) / 1000 ** 0.5 e = exact_ratio(r, r - 2, 1000) / 1000 ** 0.5 c1 &= abs(g / lit - 1) < 0.02 c2 &= abs(e / g - 1) < 0.02 print("C1/C2 r = %d, s = r - 2: Gaussian %.4f sqrt(J), literal %.4f," " exact at J = 1000 %.4f, E/G %.4f" % (r, g, lit, e, e / g)) print("C1 %s C2 %s [%.1f s]" % (c1, c2, time.time() - t0), flush=True) # (G) print("(G) the Gaussian crossings J*_s, ranks 3..40") ratios = {} c5w = 0.0 for r in range(3, 41): for s in range(r - 1): c5w = max(c5w, orth(s, r - 1 - s, kappa(s, r - 1 - s)[1])) js = jstar(r) top = max(js, key=js.get) j3 = js.get(r - 3) ratios[r] = j3 row = " r = %2d J* = %8.1f at s = r - %d;" % (r, js[top], r - top) row += " J*_(r-2) = %7.1f" % js[r - 2] if j3 is not None: row += " J*_(r-3) = %7.1f /r = %.2f /r^2 = %.3f" % ( j3, j3 / r, j3 / r ** 2) if r in JC: row += " J_c = %d, J*/J_c = %.3f" % (JC[r], js[top] / JC[r]) print(row) c5 = c5w < 1e-6 print("C5 the projection's residual against its free terms, every s," " ranks 3..40: worst relative %.2e (%s)" % (c5w, c5)) per_r = [ratios[r] / r for r in range(10, 41)] mean = sum(per_r) / len(per_r) spread = max(abs(v / mean - 1) for v in per_r) g1 = spread <= 0.25 and ratios[14] / JC[14] < 0.5 k1 = all(1.4 <= ratios[r] / r ** 2 <= 2.2 for r in range(5, 15)) print("G1 J*_(r-3)/r over r = 10..40: mean %.2f, worst deviation %.1f%%;" " J*_(r-3)/J_c at rank 14 %.3f (G1 %s)" % (mean, 100 * spread, ratios[14] / JC[14], "holds" if g1 else "fails")) print("K1 J*_(r-3)/r^2 in [1.4, 2.2] at every rank 5..14: %s (K1 %s) " "[%.1f s]" % (k1, "fires" if k1 else "misses", time.time() - t0), flush=True) # (X) print("(X) exact E = DUAL_(r-3)/H_ch over its limit G") g2, g3v = [], [] for r in range(5, 15): s = r - 3 cells = [("J_c", JC[r])] + [("%dr^2" % a, a * r * r) for a in (1, 2, 4, 8, 16)] row = " r = %2d" % r for name, J in cells: e = exact_ratio(r, s, J) g = gauss_ratio(r, s, J) row += " %s=%d: E %.3g, E/G %.3f" % (name, J, e, e / g) if name == "J_c": g2.append(e / g) if name == "16r^2": g3v.append(e / g) print(row + " [%.1f s]" % (time.time() - t0), flush=True) falls = all(b < a for a, b in zip(g2, g2[1:])) print("G2 E/G at J_c, ranks 5..14: %s (below 1 %s, falling %s)" % (["%.3f" % v for v in g2], all(v < 1 for v in g2), falls)) g3 = all(abs(v - 1) < 0.01 for v in g3v) print("G3 E/G at J = 16 r^2, ranks 5..14: %s, all within 1%% of 1: %s" % (["%.4f" % v for v in g3v], g3)) c4, k2 = arm_a(t0) print("%.1f s" % (time.time() - t0)) return not (c1 and c2 and c3 and c4 and c5 and g3) or k2 BIG = (50, 100, 200, 201, 400, 401, 800, 801) def jstar_h(r, h): """{d: J*_(r-1-d)} at d = 1..4 from the normalised h_n = He_n/sqrt(n!).""" w = PHI * PHI Mh = float(np.max(np.abs(h[r - 1] * PHI))) out = {} for d in range(1, 5): s = r - 1 - d A = np.array([[float(np.sum(h[s + i] * h[s + k] * w)) for k in range(1, d + 1)] for i in range(1, d + 1)]) rhs = -np.array([float(np.sum(h[s] * h[s + i] * w)) for i in range(1, d + 1)]) b = np.linalg.solve(A, rhs) q = h[s].copy() for i, bi in enumerate(b, 1): q = q + bi * h[s + i] f = q * PHI kh = float(np.sum(f * f) / np.sum(np.abs(f))) logfac = 0.5 * sum(np.log(j) for j in range(s + 1, r)) out[d] = 4 * np.exp((logfac + np.log(Mh / kh)) * 2 / d) return out def arm_a(t0): ranks = list(range(5, 41)) + list(BIG) need = set() for r in ranks: need |= set(range(max(0, r - 5), r)) h, prev, cur = {}, None, np.ones_like(YG) for n in range(max(need) + 1): if n in need: h[n] = cur nxt = (YG * cur - (n ** 0.5 * prev if n else 0)) / (n + 1) ** 0.5 prev, cur = cur, nxt worst = 0.0 for r in range(5, 41): a = jstar_h(r, h) g = jstar(r) for d, v in a.items(): worst = max(worst, abs(v / g[r - 1 - d] - 1)) c4 = worst < 1e-6 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) rows = {} for r in BIG: a = jstar_h(r, h) rows[r] = a bind = max(a, key=a.get) print(" r = %3d J*_d, d = 1..4: %s; binding d = %d; J*/r^1.5 = %.4f;" " J*_(r-2)/r = %.3f; J*/r^2 = %.3f [%.1f s]" % (r, ["%.1f" % a[d] for d in range(1, 5)], bind, a[bind] / r ** 1.5, a[1] / r, a[bind] / r ** 2, time.time() - t0), flush=True) top = {r: max(rows[r].values()) for r in BIG} h1 = all(abs(top[r] / r ** 1.5 / 8.128 - 1) < 0.03 for r in (400, 800)) h2 = all(max(rows[r], key=rows[r].get) == 2 for r in BIG) gaps = {r: top[r] / top[r + 1] * ((r + 1) / r) ** 1.5 - 1 for r in (200, 400, 800)} expo = np.log2(top[800] / top[400]) k2 = not (7.7 <= top[800] / 800 ** 1.5 <= 8.6 and 1.45 <= expo <= 1.55) print("H1 J*/r^1.5 at 400, 800: %.4f, %.4f against 8.128 (%s)" % (top[400] / 400 ** 1.5, top[800] / 800 ** 1.5, h1)) print("H2 binding d = 2 among d = 1..4 at every r of %s: %s" % (BIG, h2)) print("H3 J*_(r-2)/r at %s: %s" % (BIG, ["%.3f" % (rows[r][1] / r) for r in BIG])) print("H4 parity gaps at r = 200, 400, 800: %s" % ["%+.4f" % gaps[r] for r in (200, 400, 800)]) print("K2 J*(800)/800^1.5 = %.4f, local exponent 400..800 = %.4f " "(K2 %s)" % (top[800] / 800 ** 1.5, expo, "fires" if k2 else "misses")) return c4, k2 if __name__ == "__main__": sys.exit(1 if main() else 0)