"""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 (|| = |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)