bandscale.py
Checks: FLATTEN: Every rank's band of failing depths ends theorem; rule; pattern.
Cited on: FLATTEN.
To run it: save bandscale.py, with bandend.py, flatten.py beside it; needs numpy, sympy.
"""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)