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