dualtail.py
Checks: FLATTEN: Every rank's band of failing depths ends theorem; rule; FLATTEN: Past the last failure seen the champion attains the least height rule.
Cited on: FLATTEN.
To run it: save dualtail.py, with bandend.py, flatten.py beside it; needs numpy, sympy.
"""dualtail.py -- the dual certificate's tail: past which depth does the
dual bound beat the champion at every depth, at ranks 3..10?
QUESTION. bandend.py certifies the champion (1 + x)^(r-1)(x - 1)^J as
the unique least-height vector at each depth from J_c(r) through JMAX
(1000 at ranks 3, 4; 400 at 5..8; 500 at 9, 10), and closes the depths
past JMAX at ranks 3, 4 by the l2 bound's tail. At ranks 5..10 the l2
bound at the order s = r - 2 never passes H_ch by JMAX, so no l2 tail
closes there by JMAX. Does the dual bound, made explicit in J, close one?
THE OBJECTS (bandend.py's). Q_s = (1 + x)^J (1 - x)^s, c* the least
||Q_s c'||^2 cofactor with c*(1) = 1, m_s its value, DUAL_s =
m_s / ||Q_s c*||_1, I_r(J) = 2^(J+r-1) B((J + 1)/2, r/2) / pi >= H_ch.
THE BOUND (by hand, before the engine).
(a) Taylor at 1: c*(x) = sum_m b_m (1 - x)^m, m = 0..d, b_0 = 1,
b_m = (-1)^m sum_k C(k, m) c*_k. By the triangle inequality
||Q_s c*||_1 <= sum_m |b_m| L_(s+m), L_n = ||(1 + x)^J (1 - x)^n||_1.
(b) L_0 = 2^J. At n >= 1 the coefficients of (1 + x)^J (1 - x)^n are
the differences of R = (1 + x)^J (1 - x)^(n-1), zero-padded at both
ends, so L_n is R's total variation. By Krawtchouk reciprocity,
C(N, n) K_k(n) = C(N, k) K_n(k), N = J + n, their signs in k are
those of a polynomial of degree n, so they change sign at most n
times, R has at most n turning points and L_n <= 2n max|R| <=
2n I_n(J), the mean modulus of R on the circle.
(c) Divided by C(2J, J) = 4^J Gamma(J + 1/2) / (sqrt(pi) Gamma(J + 1)),
x = (J + 1)/2: I_r L_n / C(2J, J) = F_r w_n with
F_r = 2^(r-1) Gamma(r/2) [Gamma(x)/Gamma(x + r/2)]
[Gamma(J + 1)/Gamma(J + 1/2)],
w_0 = pi^(-1/2),
w_n = 2n 2^(n-1) Gamma(n/2) pi^(-3/2) Gamma(x)/Gamma(x + n/2).
SUFFICIENT: rho_s = m_s / C(2J, J) > F_r sum_m |b_m| w_(s+m).
(d) Upper bounds, every one rational once J = t^2: Wendel,
Gamma(J + 1)/Gamma(J + 1/2) <= a = sqrt(J + 1/2) <= t + 1/(4t);
Gamma(x)/Gamma(x + 1/2) <= b/x, b = sqrt(x + 1/2) =
sqrt((t^2 + 2)/2) <= beta (t + 1/t), beta = 0.70711 > 1/sqrt 2;
the integer steps of the gamma ratios exact; pi^(-e) from
pi >= 3.14159 and sqrt(pi) >= 1.77245.
(e) rho_s = det N / (D S) and b_m = B_m / S as in bandend.py's tail,
N, S and B_m polynomials in J got by exact evaluation at
integer depths and interpolation (degree <= (d + 1)(r - 1)).
S > 0 at every integer depth (G_s positive definite), so the
condition is det N / D > F_r sum_m |B_m| w_(s+m). With
sigma_m the sign of B_m past its last real root, |B_m| =
sigma_m B_m there, and the sufficient condition is one polynomial T_s(t) > 0, decided by an exact
real-root count past t_s; J_t(r, s) = ceil(t_s^2).
Asymptotically rho_s ~ J^(-s) and each |b_m| w_(s+m) F_r ~
J^(-(r+s-1)/2), so the margin is J^((r-s-1)/2): sqrt(J) at s = r - 2,
and every order's tail closes; the question is where.
PREDICTIONS (frozen before the engine; the losses are the triangle
inequality over d + 1 terms, 2n max|R| against the true variation and
H_ch against I_r, guessed a factor 3 to 6 together).
PD1 the tail closes at every rank 3..10 and every order s <= r - 2.
PD2 the binding order of the tail is s = r - 2 at every rank.
PD3 J_t(r) = max_s J_t(r, s) in [300, 2000] at rank 5 and in
[600, 5000] at rank 10, rising with r.
KILLS (printed).
KD1 (frozen at the last close) a rank in 5..10 whose tail starts past
its JMAX. Firing means the per-depth certificate must be carried
from JMAX to J_t, at the orders whose tail has not begun.
KD2 a rank in 5..10 with an order whose tail does not close, or
starts past 5000 (the per-depth carry out of reach).
CONTROLS, read before any verdict.
D1 the lemma (b): L_n <= 2n max|R| and at most n sign changes, exact,
n = 1..11, J = 0..150. A bound with 2 for 2n prints off at some
cell (planted, scratch copy).
D2 (a) exactly: ||Q_s c*||_1 <= sum_m |b_m| L_(s+m), r = 3..10, every
s <= r - 2, J in {r + 7, 60, 150}.
D3 the rational chain against lgamma: the bound's F_r sum |b_m|
w_(s+m) at J = t^2 never below the same sum computed in floats
from lgamma with the true |b_m|, and rho_s from the polynomials
equal to bandend.bounds()'s m_s / C(2J, J) exactly, at
J in {r + 7, 100, 400, 1000}. A beta of 0.7 (below 1/sqrt 2)
prints off at J = 1000 (planted, scratch copy).
D4 (positive) at ranks 3 and 4 the dual tail closes, and every depth
it covers that bandend.py also certifies is consistent: the bound
never claims a depth where the exact DUAL_s <= H_ch.
FINDINGS.
(1) Every order's dual tail closes at ranks 3..10 (PD1 holds):
J_t(r) = 35, 91, 240, 478, 817, 1267, 1839, 2542 at r = 3..10.
At ranks 3 and 4 it starts far below the l2 tail's 85 and 239.
(2) The tail's binding order is s = r - 3 from rank 4 on (r = 3:
s = 1), not r - 2: PD2 fails. The per-depth certificate's is
r - 3 from rank 5 on. J_t(5) = 240 is below the frozen
[300, 2000]: PD3 fails, the guessed losses too large.
(3) KD1 fires at ranks 6..10 (J_t past JMAX 400, 400, 400, 500,
500); KD2 misses. carry() takes the per-depth certificate from
JMAX + 1 to J_t(r, s) - 1 at each open order: 0, 77, 416, 991,
1530 and 2570 (order, depth) cells at ranks 5..10, none failing.
With bandend.py's scan, the champion is the unique least-height
vector up to sign at every depth J >= J_c = 47, 57, 97, 108, 157,
170 at ranks 5..10 (per-depth to J_t, the tail past it).
RUN RECORD. Tails, D1..D4: 602.1 s against an estimate of 3 to 6
minutes (the tails were not timed per rank first; the rank-10 cells
dominate), peak 94 MB under a 512 MB guard; D1 1661 cells, D2 132, D3
176 (true/bound at most 0.999858), D4 205, all 0 off. The carry, run
alone on the printed J_t table, 245.1 s, peak 52 MB. Planted, in
scratch copies: 2 for 2n puts D1 off at 1510 of 1661 cells; beta = 0.7
puts D3 off at 26 of the 36 cells of ranks 3..5 (their tails 2.7 s in
all, so the wall sits at ranks 6..10). Timed per (r, s) since: sympy's
interpolate, expanding symbolically, was the whole of it (36 s of
polys(7, 0) alone); exact divided differences in fractions give the
same polynomials (checked against it at (r, s) = (4, 0), (5, 1), (6, 0))
and build all 44 tails in 16.6 s. main() on that build printed D1 1661,
D2 132, D3 176 (true/bound at most 0.999858) and D4 205 cells, all 0
off, D2 and D3 reading every cell's polynomials off their nodes, and
every J_t and carried count as above; under the gate, tails, D1..D4 and
the carry, 265.7 s, peak 74.3 MB. Its exit code is red on a control
off, a tail open (PD1), KD2 or a carried cell failing.
"""
import os
import sys
import time
from fractions import Fraction
from math import comb, factorial, lgamma, log, exp, ceil, sqrt
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
from bandend import (JMAX, pmul, Q, champ_height, bounds, # noqa: E402
gram, solve)
import sympy as sp # noqa: E402
PI_LO = sp.Rational(314159, 100000) # < pi
SQPI_LO = sp.Rational(177245, 100000) # < sqrt(pi)
BETA = sp.Rational(70711, 100000) # > 1/sqrt(2)
t = sp.symbols("t", positive=True)
Jv = sp.symbols("J")
def n_k(r, k, J):
"""D * C(2J, J + k) / C(2J, J), D = (J + 1)...(J + r - 1)."""
k, K = abs(k), r - 1
out = 1
for i in range(1, k + 1):
out *= J - i + 1
for i in range(k + 1, K + 1):
out *= J + i
return out
def N_at(r, s, J):
d = r - 1 - s
g = [(-1) ** s * sum((-1) ** j * comb(2 * s, j) * n_k(r, s + m - j, 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 det_adj1(M):
"""(det M, adj(M) 1) exactly."""
n = len(M)
A = [[Fraction(v) for v in row] for row in M]
det = Fraction(1)
for c in range(n):
p = next(i for i in range(c, n) if A[i][c] != 0)
if p != c:
A[c], A[p] = A[p], A[c]
det = -det
det *= A[c][c]
for i in range(c + 1, n):
f = A[i][c] / A[c][c]
if f:
A[i] = [u - f * v for u, v in zip(A[i], A[c])]
z = solve(M, [1] * n)
return det, [det * v for v in z]
def newton(xs, ys):
"""The interpolating polynomial's coefficients, highest first, by
exact divided differences (sympy's interpolate expands symbolically,
minutes at rank 10)."""
c = [Fraction(y) for y in ys]
n = len(xs)
for j in range(1, n):
for i in range(n - 1, j - 1, -1):
c[i] = (c[i] - c[i - 1]) / (xs[i] - xs[i - j])
out = [c[-1]]
for i in range(n - 2, -1, -1):
out = [a - xs[i] * b for a, b in zip(out + [0], [0] + out)]
out[-1] += c[i]
return [sp.Rational(v.numerator, v.denominator) for v in out]
def polys(r, s):
"""detN, A = adj(N) 1 (list), S = sum A, D, as Polys in J over QQ,
by exact values at integer depths and interpolation."""
d = r - 1 - s
deg = (d + 1) * (r - 1)
pts = list(range(r + 2, r + 2 + deg + 2))
vals = [det_adj1(N_at(r, s, J)) for J in pts]
def interp(ys):
return sp.Poly(newton(pts, ys), Jv, domain="QQ")
detN = interp([v[0] for v in vals])
A = [interp([v[1][k] for v in vals]) for k in range(d + 1)]
S = sum(A[1:], A[0])
D = sp.Poly(sp.prod([Jv + i for i in range(1, r)]), Jv, domain="QQ")
return detN, A, S, D
def taylor(A, d):
"""B_m = (-1)^m sum_k C(k, m) A_k: c*(x) = sum_m (B_m/S)(1 - x)^m."""
return [sum((A[k] * comb(k, m) for k in range(m, d + 1)),
sp.Poly(0, Jv, domain="QQ")) * (-1) ** m
for m in range(d + 1)]
def gamma_half(n):
"""Gamma(n/2) = val * sqrt(pi)^e, e in {0, 1}: (val, e)."""
if n % 2 == 0:
return sp.Integer(factorial(n // 2 - 1)), 0
p = (n - 1) // 2
return sp.Rational(factorial(2 * p), 4 ** p * factorial(p)), 1
def ratio_up(n, beta=BETA):
"""Upper bound on Gamma(x)/Gamma(x + n/2), x = (t^2 + 1)/2, rational
in t (Wendel at the half step, b <= beta (t + 1/t))."""
x = (t ** 2 + 1) / 2
k = n // 2
if n % 2 == 0:
return 1 / sp.prod([x + i for i in range(k)])
b = beta * (t + 1 / t)
return b / x / sp.prod([x + sp.Rational(1, 2) + i for i in range(k)])
def pi_pow_up(e2):
"""Upper bound on pi^(e2/2), e2 in {0, -1, -2, -3}."""
return {0: sp.Integer(1), -1: 1 / SQPI_LO, -2: 1 / PI_LO,
-3: 1 / (PI_LO * SQPI_LO)}[e2]
def weights(r, s, beta=BETA):
"""Upper bounds, rational in t, on F_r w_n for n = s..r - 1."""
gr, er = gamma_half(r)
a = t + 1 / (4 * t)
F = 2 ** (r - 1) * gr * ratio_up(r, beta) * a
out = []
for n in range(s, r):
if n == 0:
out.append(F * pi_pow_up(-1 + er))
else:
gn, en = gamma_half(n)
out.append(F * 2 * n * 2 ** (n - 1) * gn * ratio_up(n, beta)
* pi_pow_up(-3 + er + en))
return out
def tail(r, s, beta=BETA):
"""J_t: from it on, rho_s > F_r sum |b_m| w_(s+m) at every real
J >= J_t, by an exact real-root count in t = sqrt(J); None if the
polynomial keeps a root past every bound. Returns (J_t, info)."""
d = r - 1 - s
detN, A, S, D = polys(r, s)
B = taylor(A, d)
# signs of B_m past their last real root
jb = 1
sig = []
for Bm in B:
if Bm.is_zero:
sig.append(0)
continue
rts = [float(hi) for (lo, hi), _ in Bm.intervals(eps=sp.Rational(1, 100))]
jb = max([jb] + [int(x) + 1 for x in rts])
sig.append(1 if Bm.LC() > 0 else -1)
W = weights(r, s, beta)
U = sum(sig[m] * B[m].as_expr().subs(Jv, t ** 2) * W[m]
for m in range(d + 1))
lhs = (detN.as_expr() / D.as_expr()).subs(Jv, t ** 2)
num, den = sp.fraction(sp.together(lhs - U))
num, den = sp.Poly(num, t, domain="QQ"), sp.Poly(den, t, domain="QQ")
# den's sign at large t is fixed; make T = num * sign(den) > 0 the test
if den.LC() < 0:
num = -num
dr = [float(hi) for (lo, hi), _ in den.intervals(eps=sp.Rational(1, 100))
if hi > 0]
rts = [float(hi) for (lo, hi), _ in num.intervals(eps=sp.Rational(1, 1000))
if hi > 0]
t0 = sp.Rational(ceil(max(rts + dr + [1.0, sqrt(jb)]) * 1000) + 1, 1000)
Jt = ceil(t0 ** 2)
ok = (num.LC() > 0 and num.count_roots(t0, None) == 0
and den.count_roots(t0, None) == 0 and num.eval(t0) > 0
and all(B[m].count_roots(Jt, None) == 0 for m in range(d + 1)
if sig[m]))
info = dict(detN=detN, A=A, S=S, D=D, B=B, sig=sig, W=W, top=float(t0),
degT=num.degree())
return (Jt if ok else None), info
# ------------------------------------------------------------ controls
def d1(two_n=True):
off = n = 0
for k in range(1, 12):
for J in range(0, 151):
R = pmul([comb(J, i) for i in range(J + 1)],
[comb(k - 1, i) * (-1) ** i for i in range(k)])
P = pmul(R, [1, -1])
sg = [v > 0 for v in P if v]
changes = sum(sg[i] != sg[i + 1] for i in range(len(sg) - 1))
bound = (2 * k if two_n else 2) * max(abs(v) for v in R)
n += 1
off += changes > k or sum(abs(v) for v in P) > bound
print("D1 the l1 lemma, n = 1..11, J = 0..150: %d cells, %d off"
% (n, off))
return off == 0
def exact_cstar(r, J, s):
d = r - 1 - s
y = solve(gram(J, s, d), [1] * (d + 1))
tot = sum(y)
return [v / tot for v in y]
def d2(cache):
off = n = 0
for r in range(3, 11):
for s in range(r - 1):
d = r - 1 - s
B, S = cache[(r, s)]["B"], cache[(r, s)]["S"]
for J in (r + 7, 60, 150):
c = exact_cstar(r, J, s)
b = [(-1) ** m * sum(comb(k, m) * c[k]
for k in range(m, d + 1))
for m in range(d + 1)]
bp = [Fraction(str(B[m].eval(J) / S.eval(J)))
for m in range(d + 1)]
l1 = sum(abs(v) for v in pmul(Q(J, s), c))
L = []
for m in range(d + 1):
k = s + m
P = pmul([comb(J, i) for i in range(J + 1)],
[comb(k, i) * (-1) ** i for i in range(k + 1)])
L.append(sum(abs(v) for v in P))
n += 1
off += (b != bp) + (l1 > sum(abs(b[m]) * L[m]
for m in range(d + 1)))
print("D2 Taylor and triangle exactly, r = 3..10, s <= r - 2, J in "
"{r + 7, 60, 150}: %d cells, %d off" % (n, off))
return off == 0
def lI(r, J):
"""ln I_r(J)."""
return ((J + r - 1) * log(2) - log(3.141592653589793)
+ lgamma((J + 1) / 2) + lgamma(r / 2) - lgamma((J + r + 1) / 2))
def d3(cache, beta=BETA):
"""The rational chain against lgamma: U(sqrt J) >= the true
F_r sum |b_m| w_(s+m), and rho from the polynomials = m_s/C(2J, J)."""
off = n = 0
worst = 0.0
for r in range(3, 11):
for s in range(r - 1):
d = r - 1 - s
info = cache[(r, s)]
W = weights(r, s, beta)
for J in (r + 7, 100, 400, 1000):
m_s, _ = bounds(r, J, s)
rho = Fraction(str(info["detN"].eval(J)
/ (info["D"].eval(J) * info["S"].eval(J))))
off += rho != m_s / comb(2 * J, J)
tv = sqrt(J)
bm = [abs(float(info["B"][m].eval(J) / info["S"].eval(J)))
for m in range(d + 1)]
ub = sum(bm[m] * float(W[m].subs(t, tv)) for m in range(d + 1))
lC = (2 * J * log(2) + lgamma(J + 0.5) - lgamma(J + 1)
- 0.5 * log(3.141592653589793))
tr = 0.0
for m in range(d + 1):
k = s + m
lL = J * log(2) if k == 0 else log(2 * k) + lI(k, J)
tr += bm[m] * exp(lI(r, J) + lL - lC)
n += 1
off += ub < tr * (1 - 1e-12)
worst = max(worst, tr / ub)
print("D3 the rational chain against lgamma and rho exactly, r = 3..10,"
" s <= r - 2, J in {r + 7, 100, 400, 1000}: %d cells, %d off, "
"true/bound at most %.6f" % (n, off, worst))
return off == 0
def d4(tails):
"""At ranks 3, 4 the bound claims a depth only where DUAL_s > H_ch."""
bad = n = 0
for r in (3, 4):
for s in range(r - 1):
jt = tails[(r, s)]
if jt is None:
bad += 1
continue
for J in range(jt, jt + 41):
n += 1
bad += not bounds(r, J, s)[1] > champ_height(r, J)
print("D4 ranks 3, 4: every order's tail closes, and DUAL_s > H_ch at "
"the 41 depths from each J_t: %d checked, %d off" % (n, bad))
return bad == 0
def carry(Jt_by_s):
"""The per-depth certificate DUAL_s > H_ch carried from JMAX + 1 to
J_t(r, s) - 1 at each order whose tail starts past JMAX: the depths
that fail it, per rank."""
out = {}
for r in range(5, 11):
fails, n = [], 0
H = {}
for s in range(r - 1):
for J in range(JMAX[r] + 1, Jt_by_s[(r, s)]):
if J not in H:
H[J] = champ_height(r, J)
n += 1
if not bounds(r, J, s)[1] > H[J]:
fails.append((s, J))
out[r] = (n, fails)
print(" r = %2d carried %5d (order, depth) cells past JMAX %d: "
"%d fail %s" % (r, n, JMAX[r], len(fails), fails[:5]),
flush=True)
return out
def main():
t0 = time.time()
tails, cache = {}, {}
for r in range(3, 11):
for s in range(r - 1):
jt, info = tail(r, s)
tails[(r, s)], cache[(r, s)] = jt, info
print("tails built (%.1f s)" % (time.time() - t0))
oks = [d1(), d2(cache), d3(cache), d4(tails)]
print("controls: %s (%.1f s)" % ("all ok" if all(oks) else "OFF",
time.time() - t0))
print("TAIL: J_t(r, s), from which the dual bound passes I_r(J) at "
"every depth")
Jt = {}
for r in range(3, 11):
row = {s: tails[(r, s)] for s in range(r - 1)}
Jt[r] = max(v if v is not None else 10 ** 9 for v in row.values())
bind = max(row, key=lambda s: row[s] if row[s] is not None else 10 ** 9)
print(" r = %2d J_t = %5d binding s = %d JMAX %d by s %s"
% (r, Jt[r], bind, JMAX[r], row))
pd1 = all(v is not None for v in tails.values())
pd2 = all(max(range(r - 1), key=lambda s: tails[(r, s)] or 10 ** 9)
== r - 2 for r in range(3, 11))
pd3 = 300 <= Jt[5] <= 2000 and 600 <= Jt[10] <= 5000 and \
all(Jt[r] <= Jt[r + 1] for r in range(5, 10))
print("PD1 every order's tail closes: %s" % ("holds" if pd1 else "fails"))
print("PD2 the binding order is r - 2 at every rank: %s"
% ("holds" if pd2 else "fails"))
print("PD3 J_t(5) in [300, 2000], J_t(10) in [600, 5000], rising: %s"
% ("holds" if pd3 else "fails"))
kd1 = [r for r in range(5, 11) if Jt[r] > JMAX[r]]
kd2 = [r for r in range(5, 11) if Jt[r] > 5000]
print("KD1 ranks whose tail starts past JMAX: %s: %s"
% (kd1, "fires" if kd1 else "misses"))
print("KD2 ranks with a tail open or past 5000: %s: %s"
% (kd2, "fires" if kd2 else "misses"))
print("CARRY: the per-depth certificate from JMAX + 1 to J_t(r, s) - 1 "
"(%.1f s)" % (time.time() - t0))
if not pd1:
print("carry skipped: a tail is open")
return True
cr = carry(tails)
held = all(not cr[r][1] for r in cr)
print("carried: no (order, depth) cell fails: %s" % held)
print("%.1f s" % (time.time() - t0))
return not (all(oks) and pd1 and held) or bool(kd2)
if __name__ == "__main__":
sys.exit(1 if main() else 0)