bandend.py
Checks: FLATTEN: Every rank's band of failing depths ends theorem; FLATTEN: Past the last failure seen the champion attains the least height rule.
Cited on: FLATTEN.
To run it: save bandend.py, with flatten.py beside it; needs numpy, sympy.
"""bandend.py -- every rank's band of failing depths ends: the champion
(1 + x)^(r-1) (x - 1)^J is the unique least-height vector of the rank-r
flattening lattice past a depth J0(r), explicit at ranks 3 and 4.
QUESTION. band.py reads the champion attaining h(J + r, J) past the last
failure at ranks 2..10 as far as its scans reach. Is there a proof, and
from which depth does it hold at ranks 3 and 4?
THE OBJECTS. x -> -x keeps heights and sends the lattice to
{c(x) (1 + x)^J : deg c <= r - 1}. Every nonzero c is (1 - x)^s c'
with a = c'(1) a nonzero integer, deg c' <= d = r - 1 - s; Q_s =
(1 + x)^J (1 - x)^s. The class s = r - 1 is a times the champion.
G_s, the (d + 1)-square Gram matrix of the shifts x^k Q_s, entries
g_m = (-1)^s [x^(J+s+m)] (1 + x)^(2J) (1 - x)^(2s).
y = G_s^-1 1, m_s = 1 / sum(y), the least ||Q_s c'||^2 over real c'
with c'(1) = 1, at c* = m_s y.
DUAL_s = 1 / ||Q_s y||_1 = m_s / ||Q_s c*||_1, a lower bound on the
height of every integer P of exact s (|<P, Q_s c*>| = |a| m_s).
L2_s = sqrt(m_s / (J + r)), the weaker lower bound.
CERTIFIED (r, J): DUAL_s > H_ch(J) at every s <= r - 2, exactly, in
integers. Then the champion is the unique minimiser up to sign.
I_r(J) = 2^(J+r-1) B((J + 1)/2, r/2) / pi >= H_ch(J).
J0(r), the scan's: the least depth from which every depth through
JMAX is CERTIFIED. A sufficient condition's start, not the depth
where the champion starts to attain h (band.py's J_ch).
THE PROOF AT EVERY RANK (by hand, before the engine). Drop the circle
outside |theta| <= 2K/sqrt(J) (the integrand of ||P||^2 is >= 0), put
theta = 2u/sqrt(J) and write c'(x) = g(w), w = sqrt(J)(x - 1), g(0) = a.
The window's Gram matrix in the basis w^m, divided by
4^J (4/J)^s (2/sqrt(J)) / (2 pi), tends entrywise to that of
(2iu)^m against u^(2s) e^(-u^2) on |u| <= K, positive definite; the
least form at g(0) = a is |a|^2 / (G^-1)_00, continuous there, so no
cofactor escapes by huge coefficients. Hence ||P||^2 >= (1 - o(1))
lambda 4^J 4^s J^(-s-1/2) / pi, the height >= c 2^J J^(-s/2-3/4), while
H_ch <= I_r(J) ~ c' 2^J J^(-r/2): a margin J^(1/4) at s = r - 2. For
every r there is J0(r) past which the champion is the unique minimiser
up to sign, h(J + r, J) = H_ch(J), and the cell cannot fail.
PREDICTIONS (Gaussian limits, TRANSPLANTS). The binding class is
s = r - 2, DUAL/H_ch about (sqrt(J)/2) kappa_(r-2)/M_(r-1), kappa_s =
int He_s^2 e^(-y^2) / int |He_s| e^(-y^2/2) and M_n = max |He_n(y)
e^(-y^2/2)| (He the probabilists' Hermite polynomials): 0.22 sqrt(J)
at r = 3 (J near 20), 0.20 sqrt(J) at r = 4 (near 25), about 0.12 sqrt(J)
at r = 5 and 0.10 sqrt(J) at r = 6. The l2 bound
against H_ch near J = 82 (r = 3) and 130 (r = 4); against I_4, 234.
PR1 J0(3) in [12, 40], J0(4) in [15, 50].
PR2 J0(r) > J_hi(r) + 1 at every r = 5..10, rising with r.
PR3 the tail's J_t <= 400 at r = 3, 4.
KILLS. K1 a control off. K2 J0(4) > 140. K3 a J_t past JMAX at r = 3, 4.
CONTROLS, read before any verdict.
C1 m_s against a float least-squares solve, r = 3, 4, J = 4..40.
C2 least height and least ||P||^2 per exact s by ball enumeration,
r = 3, 4, J = 4..40: DUAL and L2 never above the least height,
m_s never above the least norm, the minimum over s equal to
flatten.least_height's h.
C3 at r = 3..5, J = 2..40, no certification where h < H_ch.
C4 DUAL at s = r - 1 never above H_ch.
C5 the tail's arm: its rho_s(J) equal to bounds()'s m_s / C(2J, J)
exactly, and its lower bound F_low never above F = m_s / ((J + r)
I_r(J)^2) by lgamma, at J = r + 7, 100, J_t, 1000. A rho with one
factor too many in D, and pi bounded from above (3.15), each make
it print off (planted, scratch copies).
C2 and C5 carry the soundness; C3 and C4 catch only gross faults (C3's
cells sit far below every crossing, C4 is Holder's at d = 0).
TAIL (r = 3, 4) the l2 inequality m_s / (J + r) > I_r(J)^2 at every
integer J >= J_t, by Wendel's bounds on the gamma ratios and an
exact real-root count of one polynomial per s.
FINDINGS.
(1) THEOREM at r = 3 and 4: the champion is the unique least-height
vector up to sign at every depth J >= 19 (r = 3) and J >= 26
(r = 4). Exact certificates at every J = J0..1000; the l2 tail
from J_t = 85 (r = 3, s = 1) and 239 (r = 4, s = 2), the other
classes' tails from 9, 9 and 37.
(2) The dual bound crosses where its Gaussian limit says: DUAL/H_ch
/sqrt(J) at J = 1000 is 0.2218 (r = 3) and 0.1986 (r = 4)
against the frozen 0.22 and 0.20. The l2 bound crosses at 83
and 135.
(3) Ranks 5..10: J0 = 47, 57, 97, 108, 157, 170 through J = 400
(r = 5..8) and 500 (r = 9, 10), each past band.py's J_hi + 1 =
31, 34, 59, 61, 89, 107 (PR2 holds), and past the last depth
band.py reads, J_hi + 12 = 42, 45, 70, 72, 100, 118, so a gap of depths
stands between each scan and its certificate, and the l2 tail
does not begin by JMAX: the l2 bound at s = r - 2 never crosses
by then at r = 5..10.
(4) TRANSPLANTS missed: at r = 5..10 the binding class is s = r - 3,
not r - 2 (r = 5: s = 2 crosses at 47, s = 3 at 43). The ratio
at JMAX, 0.1466 (r = 5) and 0.1441 (r = 6), against the frozen
0.12 and 0.10: a slip in evaluating kappa and M by hand, not a
miss of the limit.
PR1, PR2, PR3 hold; K1, K2, K3 miss.
RUN RECORD. The first full run died at C4 on C(2J, k) with k < 0 at
J < s (a zero binomial, now guarded). The second carried the dual sum
in Fractions; replaced by one integer sum over Y = den y (0.03 s a cell
at r = 10, J = 500). The third ran the tail in sympy's rational
functions and was refused at the 512 MB ceiling at r = 4, s = 0; the
tail now clears the common denominator (J + 1)...(J + r - 1) and works
in polynomials alone (50 MB). Controls: C1 worst relative 1.47e-15;
C2 259 classes, 0 bounds above their least, 0 cells off h; C3 42 cells
with h < H_ch, 0 certified; C4 0 of 952; C5 20
depths, 0 off, F_low/F at most 0.999874. Wall 42.2 s, peak 89 MB.
"""
import os
import sys
import time
from fractions import Fraction
from math import comb, exp, factorial, gcd, lgamma, log, pi
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
from flatten import least_height, height, lll # noqa: E402
JMAX = {3: 1000, 4: 1000, 5: 400, 6: 400, 7: 400, 8: 400, 9: 500, 10: 500}
# band.py's FINDINGS F5: the last failing depth its scan reads, ranks 5..10
JHI = {5: 30, 6: 33, 7: 58, 8: 60, 9: 88, 10: 106}
def pmul(a, b):
out = [0] * (len(a) + len(b) - 1)
for i, x in enumerate(a):
if x:
for j, y in enumerate(b):
out[i + j] += x * y
return out
def Q(J, s):
return pmul([comb(J, i) for i in range(J + 1)],
[comb(s, i) * (-1) ** i for i in range(s + 1)])
def champ_height(r, J):
p = [comb(J, i) * (-1) ** (J - i) for i in range(J + 1)]
for _ in range(r - 1):
p = pmul(p, [1, 1])
return height(p)
def C(n, k):
return comb(n, k) if 0 <= k <= n else 0
def gram(J, s, d):
g = [(-1) ** s * sum((-1) ** j * comb(2 * s, j) * C(2 * J, J + s + m - j)
for j in range(2 * s + 1)) for m in range(d + 1)]
return [[g[abs(i - k)] for k in range(d + 1)] for i in range(d + 1)]
def solve(G, b):
n = len(G)
A = [[Fraction(G[i][k]) for k in range(n)] + [Fraction(b[i])]
for i in range(n)]
for c in range(n):
p = next(i for i in range(c, n) if A[i][c] != 0)
A[c], A[p] = A[p], A[c]
for i in range(n):
if i != c and A[i][c] != 0:
f = A[i][c] / A[c][c]
A[i] = [u - f * v for u, v in zip(A[i], A[c])]
return [A[i][n] / A[i][i] for i in range(n)]
def bounds(r, J, s):
"""(m_s, DUAL_s) as Fractions; y is scaled to integers Y = den y
so the l1 norm is one integer sum."""
d = r - 1 - s
y = solve(gram(J, s, d), [1] * (d + 1))
den = 1
for v in y:
den = den * v.denominator // gcd(den, v.denominator)
Y = [int(v * den) for v in y]
m = 1 / sum(y)
l1 = sum(abs(v) for v in pmul(Q(J, s), Y))
return m, Fraction(den, l1)
def certified(r, J, H=None):
H = champ_height(r, J) if H is None else H
return all(bounds(r, J, s)[1] > H for s in range(r - 1))
# ------------------------------------------------------------ controls
def c1():
import numpy as np
worst = 0.0
for r in (3, 4):
for J in range(4, 41):
for s in range(r):
d = r - 1 - s
q = np.array(Q(J, s), dtype=float)
scale = np.abs(q).max()
B = np.zeros((d + 1, len(q) + d))
for k in range(d + 1):
B[k, k:k + len(q)] = q / scale
# c' = e0 + N z with sum(c') = 1: N spans sum zero
N = np.zeros((d + 1, d))
for k in range(d):
N[0, k], N[k + 1, k] = -1.0, 1.0
e0 = np.zeros(d + 1)
e0[0] = 1.0
if d:
z = np.linalg.lstsq((B.T @ N), -(B.T @ e0), rcond=None)[0]
c = e0 + N @ z
else:
c = e0
val = float(np.sum((B.T @ c) ** 2)) * scale * scale
m = float(bounds(r, J, s)[0])
worst = max(worst, abs(val - m) / m)
print("C1 m_s against float least squares, r = 3, 4, J = 4..40, "
"every s: worst relative %.2e" % worst)
return worst < 1e-9
def exact_s_minima(r, J, s):
"""(least height, least ||P||^2) over integer c' with c'(1) != 0,
by Schnorr-Euchner enumeration of the shifts of Q_s. The height
search walks the ball of radius (best - 1) sqrt(width) shrinking on
each new best, so it misses no vector of lower height; the norm
search the ball of radius^2 = best norm."""
d = r - 1 - s
q = Q(J, s)
n, M = d + 1, len(q) + d
rows = []
for k in range(n):
v = [0] * M
v[k:k + len(q)] = q
rows.append(v)
# reduce, then read each reduced vector's c' coordinates exactly
b, dd, lam = lll(rows)
T = []
for v in b:
coef = solve([[sum(rows[i][t] * rows[k][t] for t in range(M))
for k in range(n)] for i in range(n)],
[sum(rows[i][t] * v[t] for t in range(M))
for i in range(n)])
assert all(c.denominator == 1 for c in coef)
T.append([int(c) for c in coef])
Bs = [dd[i + 1] / dd[i] for i in range(n)]
mu = [[lam[k][j] / dd[j + 1] for j in range(k)] for k in range(n)]
out = []
for mode in ("h", "n"):
best = [height(q) if mode == "h" else sum(t * t for t in q)]
R2 = [((best[0] - 1) ** 2 * M if mode == "h" else best[0])
* (1 + 1e-9) + 1e-9]
x = [0] * n
def leaf():
cp = [sum(x[i] * T[i][k] for i in range(n)) for k in range(n)]
if sum(cp) == 0:
return
v = [0] * M
for i in range(n):
if x[i]:
for t in range(M):
v[t] += x[i] * b[i][t]
val = height(v) if mode == "h" else sum(t * t for t in v)
if 0 < val < best[0]:
best[0] = val
R2[0] = (((val - 1) ** 2 * M if mode == "h" else val)
* (1 + 1e-9) + 1e-9)
def rec(i, partial):
c = -sum(mu[j][i] * x[j] for j in range(i + 1, n))
x0 = round(c)
sg = 1 if c >= x0 else -1
t = 0
while True:
xi = (x0 if t == 0 else x0 + sg * ((t + 1) // 2)
if t % 2 else x0 - sg * (t // 2))
t += 1
p = partial + (xi - c) ** 2 * Bs[i]
if p > R2[0]:
break
x[i] = xi
if i == 0:
if any(x):
leaf()
else:
rec(i - 1, p)
x[i] = 0
rec(n - 1, 0.0)
out.append(best[0])
return out[0], out[1]
def c2():
bad = hoff = 0
cells = 0
for r in (3, 4):
for J in range(4, 41):
hs = []
for s in range(r):
lh, ln = exact_s_minima(r, J, s)
m, D = bounds(r, J, s)
L2sq = m / (J + r)
bad += (D > lh) + (L2sq > lh * lh) + (m > ln)
hs.append(lh)
cells += 1
hoff += min(hs) != least_height(J + r, J)[0]
print("C2 exact-s minima, r = 3, 4, J = 4..40: %d (r, J, s) classes, "
"%d bounds above their least, %d cells where min over s != h"
% (cells, bad, hoff))
return bad == 0 and hoff == 0
def c3():
n = wrong = 0
for r in (3, 4, 5):
for J in range(2, 41):
H = champ_height(r, J)
if least_height(J + r, J)[0] < H:
n += 1
wrong += certified(r, J, H)
print("C3 cells r = 3..5, J = 2..40 with h < H_ch: %d, certified "
"among them: %d" % (n, wrong))
return n > 0 and wrong == 0
def c4():
above = n = 0
for r in range(3, 11):
for J in range(2, 121):
n += 1
above += bounds(r, J, r - 1)[1] > champ_height(r, J)
print("C4 DUAL at s = r - 1 above H_ch: %d of %d (r = 3..10, "
"J = 2..120)" % (above, n))
return above == 0
# ---------------------------------------------------------- the scan
def scan(r):
"""J0: least J with every depth J0..JMAX certified; first crossings
per s of DUAL and of L2 against H_ch, as the last failing depth + 1."""
jm = JMAX[r]
last_fail = r - 2
last_s = {s: r - 2 for s in range(r - 1)}
last_l2 = {s: r - 2 for s in range(r - 1)}
ratio = None
for J in range(r - 1, jm + 1):
H = champ_height(r, J)
ok = True
for s in range(r - 1):
m, D = bounds(r, J, s)
if not D > H:
ok = False
last_s[s] = J
if not m > (J + r) * H * H:
last_l2[s] = J
if s == r - 2 and J == jm:
ratio = float(D / H) / jm ** 0.5
if not ok:
last_fail = J
return last_fail + 1, {s: v + 1 for s, v in last_s.items()}, \
{s: v + 1 for s, v in last_l2.items()}, ratio
# ------------------------------------------------------------ the tail
def tail(r, s, start):
"""Least integer J_t >= start such that m_s/(J + r) > I_r(J)^2 at
every integer J >= J_t, from T(J) > 0 at every real J >= J_t by an
exact real-root count; None when T keeps a real root past J_t.
D = (J + 1) ... (J + r - 1) clears every C(2J, J + k)/C(2J, J),
|k| <= r - 1, so G / C(2J, J) = N / D with N a polynomial matrix and
rho_s = m_s / C(2J, J) = det N / (D S), S = 1^T adj(N) 1."""
import sympy as sp
J = sp.symbols("J")
d, K = r - 1 - s, r - 1
def n_k(k): # D * C(2J, J + k) / C(2J, J)
k = abs(k)
return sp.prod([J - i + 1 for i in range(1, k + 1)]
+ [J + i for i in range(k + 1, K + 1)])
g = [sp.expand((-1) ** s * sum((-1) ** j * comb(2 * s, j)
* n_k(s + m - j)
for j in range(2 * s + 1)))
for m in range(d + 1)]
N = sp.Matrix(d + 1, d + 1, lambda i, k: g[abs(i - k)])
detN = sp.Poly(N.det(method="berkowitz"), J, domain="QQ")
S = sp.Poly(sum(N.adjugate(method="berkowitz")), J, domain="QQ")
D = sp.Poly(sp.prod([J + i for i in range(1, K + 1)]), J, domain="QQ")
# F = rho/(J + r) * pi^(3/2)/(4^(r-1) Gamma(r/2)^2)
# * Gamma(J + 1/2)/Gamma(J + 1) * [Gamma(x + r/2)/Gamma(x)]^2,
# x = (J + 1)/2. Wendel: Gamma(J+1/2)/Gamma(J+1) >= (J + 1/2)^(-1/2);
# r = 2p: Gamma(x + p)/Gamma(x) = x (x + 1) ... (x + p - 1); r = 2p + 1:
# Gamma(x + p + 1/2)/Gamma(x) >= (x + 1/2) ... (x + p - 1/2) x / sqrt(x + 1/2).
# F > 1 follows from F_low^2 > 1, F_low^2 = rho^2/(J + r)^2 * C2 *
# gam2 / (J + 1/2), gam2 = gn / gd; cleared of its positive
# denominators: T = detN^2 C2 gn - (J + r)^2 (J + 1/2) D^2 S^2 gd > 0.
x = (J + 1) / 2
piL = sp.Rational(314159, 100000)
half = sp.Rational(1, 2)
p = r // 2
if r % 2 == 0: # Gamma(r/2) = (p - 1)!
C2 = piL ** 3 / (sp.Integer(4) ** (2 * (r - 1)) * factorial(p - 1) ** 4)
gn, gd = sp.prod([x + i for i in range(p)]) ** 4, sp.Integer(1)
else: # Gamma(r/2) = sqrt(pi) qq, Wendel at the half step
qq = sp.Rational(factorial(2 * p), 4 ** p * factorial(p))
C2 = piL / (sp.Integer(4) ** (2 * (r - 1)) * qq ** 4)
gn = (sp.prod([x + half + i for i in range(p)]) ** 2 * x ** 2) ** 2
gd = (x + half) ** 2
gn, gd = sp.Poly(gn, J, domain="QQ"), sp.Poly(gd, J, domain="QQ")
T = (detN ** 2 * gn * C2
- sp.Poly((J + r) ** 2 * (J + half), J, domain="QQ")
* D ** 2 * S ** 2 * gd)
top = max([float(hi) for (lo, hi), _ in T.intervals(eps=sp.Rational(1, 1000))]
+ [0.0])
Jt = max(start, int(top) + 1)
def rho(j):
return Fraction(str(detN.eval(j) / (D.eval(j) * S.eval(j))))
def f_low(j):
return float(sp.sqrt(detN.eval(j) ** 2 * gn.eval(j) * C2
/ ((j + r) ** 2 * (j + half) * D.eval(j) ** 2
* S.eval(j) ** 2 * gd.eval(j))))
if T.count_roots(Jt, None) == 0 and T.eval(Jt) > 0:
return Jt, top, rho, f_low
return None, top, rho, f_low
def c5(tails):
"""The tail's arm against the exact per-depth values: rho_s(J)
equal to m_s / C(2J, J) from bounds() exactly, and F_low <= F, F =
m_s / ((J + r) I_r(J)^2) by lgamma, at J = r + 7, 100, J_t, 1000."""
off = n = 0
worst = 0.0
for (r, s), (jt, _, rho, f_low) in tails.items():
for J in sorted({r + 7, 100, jt or 100, 1000}):
m = bounds(r, J, s)[0]
n += 1
off += rho(J) != m / comb(2 * J, J)
lI = ((J + r - 1) * log(2) - log(pi) + lgamma((J + 1) / 2)
+ lgamma(r / 2) - lgamma((J + r + 1) / 2))
F = exp(log(m.numerator) - log(m.denominator) - log(J + r)
- 2 * lI)
fl = f_low(J)
off += fl > F * (1 + 1e-12)
worst = max(worst, fl / F)
print("C5 the tail's rho against bounds() exactly and F_low <= F, "
"r = 3, 4, every s, %d depths: %d off, F_low/F at most %.6f"
% (n, off, worst))
return off == 0
def main():
t0 = time.time()
oks = [c1(), c2(), c3(), c4()]
arms = {(r, s): tail(r, s, r - 1) for r in (3, 4) for s in range(r - 1)}
oks.append(c5(arms))
print("controls: %s (%.1f s)" % ("all ok" if all(oks) else "OFF",
time.time() - t0))
print("SCAN: J0 = least depth from which DUAL certifies every depth "
"through JMAX")
J0 = {}
for r in range(3, 11):
j0, ds, l2, rt = scan(r)
J0[r] = j0
print(" r = %2d J0 = %4d (JMAX %d) DUAL crossings by s %s "
"L2 crossings by s %s DUAL/H_ch/sqrt(J) at JMAX %.4f"
% (r, j0, JMAX[r], ds, l2, rt))
print("TAIL: J_t per s with m_s/(J+r) > I_r(J)^2 at every integer "
"J >= J_t")
tails = {}
for (r, s), (jt, top, _, _) in arms.items():
tails[(r, s)] = jt
print(" r = %d s = %d J_t = %s (largest real root %.2f)"
% (r, s, jt, top))
for r in (3, 4):
jt = max(tails[(r, s)] or 10 ** 9 for s in range(r - 1))
print(" r = %d: tail from %d, per-depth certificates through %d: "
"%s" % (r, jt, JMAX[r], "CLOSED" if jt <= JMAX[r] else "OPEN"))
pr1 = 12 <= J0[3] <= 40 and 15 <= J0[4] <= 50
pr2 = all(J0[r] > JHI[r] + 1 for r in JHI) and \
all(J0[r] < J0[r + 1] for r in range(5, 10))
pr3 = all((tails[(r, s)] or 10 ** 9) <= 400
for r in (3, 4) for s in range(r - 1))
print("PR1 J0(3) in [12, 40], J0(4) in [15, 50]: %s" % ("holds" if pr1
else "fails"))
print("PR2 J0(r) > J_hi(r) + 1 (band.py's J_hi %s), rising, r = 5..10: "
"%s" % (JHI, "holds" if pr2 else "fails"))
print("PR3 every J_t <= 400 at r = 3, 4: %s" % ("holds" if pr3
else "fails"))
k1 = not all(oks)
print("K1 %s" % ("fires" if k1 else "misses"))
k2 = J0[4] > 140
print("K2 J0(4) = %d: %s" % (J0[4], "fires" if k2 else "misses"))
k3 = any((tails[(r, s)] or 10 ** 9) > JMAX[r]
for r in (3, 4) for s in range(r - 1))
print("K3 %s" % ("fires" if k3 else "misses"))
print("%.1f s" % (time.time() - t0))
return k1 or k2 or k3 or not (pr1 and pr2 and pr3)
if __name__ == "__main__":
sys.exit(1 if main() else 0)