"""vartail.py -- the dual tail with ||Q_s c*||_1 bounded whole: does every order close at ranks 11..14, and from what depth? QUESTION. dualtail.py's tail bounds ||Q_s c*||_1 term by term (Taylor at 1, the triangle inequality, the variation lemma per term). At ranks 11..14 that tail starts past depth 500 at several orders (821 at rank 11, s = r - 5), where the certificate stops. Bound the norm whole. THE BOUND (by hand, before the engine). (a) c* is palindromic: G_s is symmetric Toeplitz, so y = G_s^-1 1 is fixed by reversal. (b) P = Q_s c* = (1 + x)^J p, p = (1 - x)^s c* of degree r - 1, and multiplying by 1 + x adds no sign change (each nonzero coefficient f_k + f_(k-1) takes the sign of f_k or f_(k-1), and the sign block it takes it from never moves back), so V = var(P) <= r - 1. (c) At s >= 1, P = (1 - x) R, R = (1 + x)^J (1 - x)^(s-1) c*: P's coefficients are the differences of 0, R_0, .., R_n, 0, which run monotone in at most V + 1 runs, the end runs from or to 0, so ||P||_1 <= 2 V max|R_k| <= 2 V mean|R(e^(i theta))|. (d) d = r - 1 - s = 2e + eps: c* = (1 + x)^eps ct, ct palindromic of degree 2e, ct = sum_i gamma_i x^(e-i) (1 - x)^(2i). On the circle (1 - x)^2 = -4xU, U = sin^2(theta/2), so |c*| = |1 + x|^eps |g(U)|, g(U) = sum_i gamma_i (-4U)^i, and mean|R| = I_s(J + eps) E|g(U)|, U ~ Beta(s/2, (J + 1 + eps)/2), I_n(J) = 2^(J+n-1) B((J+1)/2, n/2)/pi. (e) E|g| <= sqrt(E g^2) (exact at d = 1, g = 1/2), E g^2 = sum gamma_i gamma_j (-4)^(i+j) E U^(i+j), E U^n = prod_(k K^2 E g^2, K = 2V F_r 2^(s-1+eps) Gamma(s/2) pi^(-3/2) Gamma(y)/Gamma(y + s/2), y = (J + 1 + eps)/2, both sides positive (m_s = 1/(1' G_s^-1 1) > 0), so squaring is exact; each Gamma ratio bounded above by Wendel at the half step, sqrt(y + 1/2) <= beta (t + (1 + eps/2)/t), beta > 1/sqrt2, J = t^2; times S^2 it is one polynomial in t, its positive real roots counted exactly. s = 0 keeps the triangle tail. The margin is dualtail.py's, J^((r-s-1)/2): the loss moves the tail at r - 3 by itself, at r - 2 by its square. PREDICTIONS (frozen before the engine). PV1 at r - 3, ranks 5..10, at dualtail's J_t: 2V mean|R|/||P||_1, V = r - 1, in [2.5, 5], below l2 l3 at ranks 6..10; the Cauchy-Schwarz factor sqrt(E g^2)/E|g| in [1, 1.3]. PV2 the new r - 3 tail starts at most 0.5 dualtail's at ranks 7..10. PV3 at ranks 11..14 every order's tail closes, the largest start in [700, 2500], at order r - 3 or r - 2. PV4 at ranks 11..14 the new start is below the triangle's at every order 1 <= s <= r - 3. KILL (printed): the new r - 3 tail starts at or past dualtail's at a rank in 7..10, or a rank in 11..14 with an order open or a start past 5000. CONTROLS, read before any verdict. V1 (b)-(d) exactly: r = 5..8, s = 1..r - 2, J in {r + 7, 60, 150}: var(P) <= var(p) <= r - 1, ||P||_1 <= 2 var(P) max|R|, max|R| <= the Beta form, the Beta form against mpmath quadrature of |R| on the circle at J <= 60, and the peeled c* re-expanded equal to c*. PLANT: V for 2V, off at some cell. V2 the rational chain against lgamma: K^2 E g^2 at J = t^2 never below the same from lgamma with E g^2 from the exact c* peeled, the assembled E g^2 (vtail's polynomials) equal to that exactly, and the polys' rho_s equal to m_s/C(2J, J); ranks 5..14, every s >= 1, J in {r + 7, 100, 400, 1000}. PLANT: beta = 0.7, off at some cell. V3 (positive) ranks 5..14: the tails taken claim a depth only where the exact DUAL_s > H_ch, at the 41 depths from each order's start (the smaller of the two bounds; the triangle's alone at s = 0). FINDINGS (copied from the printed output). (1) Controls: V1 54 cells, 0 off (the plant V for 2V: 46 off); V2 300 cells, 0 off, true/bound at most 0.999744 (the plant beta = 0.7: 165 off), the assembled E g^2 equal to the exact c*'s at every cell; V3 3485 depths, every order's taken start, 0 off. (2) PV1 holds. At s = r - 3 and dualtail's J_t, ranks 5..10: the loss 2(r - 1) mean|R|/||P||_1 = 2.5718, 2.8443, 2.6739, 2.9843, 2.9017, 3.1833; the Cauchy-Schwarz factor 1.1384 .. 1.2049; their product 2.9277 .. 3.8357 against l2 l3 = 4.9187 .. 13.7008. There var(P) is r - 1 at even r and r - 3 at odd r (2, 5, 4, 7, 6, 9). (3) PV2, PV3, PV4 hold; the KILL misses. The whole-norm r - 3 tail starts at 313, 427, 559, 709 against dualtail's 817, 1267, 1839, 2542 (ranks 7..10). Every order closes at ranks 5..14 and the binding order moves to r - 2, where the whole-norm starts above the triangle at ranks 5..12 (by 29% of it at rank 5, 0.5% at 12) and below it at 13 and 14: J_t(r) = 163, 261, 383, 526, 693, 880, 1091, 1322, 1571, 1833 at r = 5..14. (4) The carry: the certificate DUAL_s > H_ch at every (order, depth) cell from 501 to J_t(r, s) - 1, 966, 1383, 1836 and 2320 cells at ranks 11..14, none failing. With the certificate from J_c through 500 (bandgrowth.py), the champion is the only least-height vector up to sign at every depth J >= J_c at ranks 11..14. Not used, by hand after the run: at odd d, p = (1 + x)(1 - x)^s ct, so V <= r - 2; at d = 1 that would cut the binding r - 2 start by about ((r - 2)/(r - 1))^2. RUN RECORD. Tails at ranks 5..14 in 170.6 s; controls and reads to 237.8 s; the carry to 383.4 s; peak 169.5 MB under a 512 MB guard, against an estimate of 20 minutes (the carry's cells 0.02 to 0.06 s, timed first). V1 first read 28 of 54 cells off: the quadrature of |R| on the circle was not split at g's zeros (off by 1e-5 to 1e-7 at d >= 2, exact at d = 1); split there, 0 off, the tolerance 1e-12 kept. A first V2 read E g^2 off vtail's own gamma_i on both sides and V3 only the whole-norm starts; both widened and rerun, every other print the same. Planted, in a scratch copy: the i + j = 1 term of vtail's E g^2 doubled puts V2 off at 40 checks over ranks 5, 6 (0 unplanted). Exit red on a control off, an order open, a carried cell failing or the KILL. """ import os import sys import time from fractions import Fraction from math import comb, lgamma, log, ceil, sqrt, pi os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import mpmath as mp # noqa: E402 import sympy as sp # noqa: E402 from bandend import Q, pmul, champ_height, bounds # noqa: E402 from dualtail import (polys, tail as tri_tail, gamma_half, pi_pow_up, # noqa: E402 ratio_up, exact_cstar, BETA, t, Jv) mp.mp.dps = 30 JT = {5: 240, 6: 478, 7: 817, 8: 1267, 9: 1839, 10: 2542} # dualtail.py def peel(c, d): """(eps, gamma): c = (1 + x)^eps sum_i gamma_i x^(e-i) (1 - x)^(2i); c a palindromic coefficient list (numbers or Polys), d = 2e + eps.""" eps, e = d % 2, d // 2 ct = list(c) if eps: q = [ct[0]] for k in range(1, d): q.append(ct[k] - q[-1]) ct = q rem = list(ct) gam = [None] * (e + 1) for i in range(e, -1, -1): g = rem[e - i] gam[i] = g for j in range(2 * i + 1): rem[e - i + j] = rem[e - i + j] - g * comb(2 * i, j) * (-1) ** j return eps, gam def moments(a, b, n): """E U^k, k = 0..n, U ~ Beta(a, b).""" out = [1] for k in range(n): out.append(out[-1] * (a + k) / (a + b + k)) return out def ratio_y(n, eps, beta=BETA): """Upper bound on Gamma(y)/Gamma(y + n/2), y = (t^2 + 1 + eps)/2.""" y = (t ** 2 + 1 + eps) / 2 k = n // 2 if n % 2 == 0: return 1 / sp.prod([y + i for i in range(k)]) b = beta * (t + (1 + sp.Rational(eps, 2)) / t) return b / y / sp.prod([y + sp.Rational(1, 2) + i for i in range(k)]) def K_up(r, s, eps, V, beta=BETA): """Upper bound, rational in t, on K = 2V F_r 2^(s-1+eps) Gamma(s/2) pi^(-3/2) Gamma(y)/Gamma(y + s/2).""" gr, er = gamma_half(r) gs, es = gamma_half(s) F = 2 ** (r - 1) * gr * ratio_up(r, beta) * (t + 1 / (4 * t)) return (2 * V * F * 2 ** (s - 1 + eps) * gs * ratio_y(s, eps, beta) * pi_pow_up(-3 + er + es)) def tpoly(P): """A Poly in J as a Poly in t, J = t^2.""" cs = P.all_coeffs() out = [] for c in cs[:-1]: out += [c, 0] out.append(cs[-1]) return sp.Poly(out, t, domain="QQ") def ev2(P, J): """A Poly in t of even powers only, at t^2 = J, exactly.""" out = Fraction(0) for (k,), c in P.terms(): assert k % 2 == 0 out += Fraction(int(c.p), int(c.q)) * J ** (k // 2) return out def var(seq): sg = [v > 0 for v in seq if v] return sum(sg[i] != sg[i + 1] for i in range(len(sg) - 1)) def vtail(r, s, beta=BETA, V=None): """J_t of the whole-norm bound at order s >= 1 (None if open), info.""" d = r - 1 - s V = r - 1 if V is None else V detN, A, S, D = polys(r, s) eps, G = peel(A, d) e = d // 2 a = sp.Rational(s, 2) bb = (t ** 2 + 1 + eps) / 2 # common denominator of E U^n, n <= 2e: prod_(k < 2e) (a + bb + k) Dm = sp.prod([a + bb + k for k in range(2 * e)]) Mn = [sp.Poly(sp.expand(sp.prod([a + k for k in range(n)]) * sp.prod([a + bb + k for k in range(n, 2 * e)])), t, domain="QQ") for n in range(2 * e + 1)] Gt = [tpoly(g) for g in G] Eg = sp.Poly(0, t, domain="QQ") for i in range(e + 1): for j in range(e + 1): Eg += Gt[i] * Gt[j] * Mn[i + j] * (-4) ** (i + j) Kn, Kd = sp.fraction(sp.together(K_up(r, s, eps, V, beta) ** 2)) Kn, Kd = sp.Poly(Kn, t, domain="QQ"), sp.Poly(Kd, t, domain="QQ") Dmp = sp.Poly(sp.expand(Dm), t, domain="QQ") dN, dD = tpoly(detN), tpoly(D) num = dN ** 2 * Kd * Dmp - dD ** 2 * Kn * Eg # denominators D^2 S^2 Kd Dm: Kd and Dm are products of t and # positive shifts of t^2, read below for no root past t0 if Kd.LC() < 0: num, Kd = -num, -Kd rts = [float(hi) for (lo, hi), _ in num.intervals(eps=sp.Rational(1, 1000)) if hi > 0] dr = [float(hi) for P in (Kd, Dmp) for (lo, hi), _ in P.intervals(eps=sp.Rational(1, 100)) if hi > 0] t0 = sp.Rational(ceil(max(rts + dr + [1.0]) * 1000) + 1, 1000) Jt = ceil(t0 ** 2) ok = (num.LC() > 0 and num.count_roots(t0, None) == 0 and Kd.count_roots(t0, None) == 0 and Dmp.count_roots(t0, None) == 0 and num.eval(t0) > 0) return (Jt if ok else None), dict(detN=detN, A=A, S=S, D=D, G=G, eps=eps, Eg=Eg, Dm=Dmp, degT=num.degree()) # ------------------------------------------------------------ exact reads def g_roots(gam): """g's coefficients (highest first, in U) and its real roots in (0, 1).""" g = [mp.mpf(v.numerator) / v.denominator if isinstance(v, Fraction) else mp.mpf(v) for v in gam] poly = [g[i] * (-4) ** i for i in range(len(g))][::-1] roots = [rt.real for rt in (mp.polyroots(poly, maxsteps=200, extraprec=60) if len(poly) > 1 else []) if abs(rt.imag) < 1e-20 and 0 < rt.real < 1] return poly, sorted(roots) def Eabs_g(gam, a, b): """E|g(U)|, U ~ Beta(a, b), by quadrature split at g's roots in (0, 1).""" poly, roots = g_roots(gam) pts = [mp.mpf(0)] + roots + [mp.mpf(1)] lB = mp.log(mp.beta(a, b)) def f(u): return abs(mp.polyval(poly, u)) * mp.exp( (a - 1) * mp.log(u) + (b - 1) * mp.log(1 - u) - lB) return mp.quad(f, pts) def I(n, J): return (mp.mpf(2) ** (J + n - 1) * mp.beta(mp.mpf(J + 1) / 2, mp.mpf(n) / 2) / mp.pi) def exact_read(r, s, J): """At one depth: V, ||P||_1, max|R|, the Beta form mean|R|, E|g|, sqrt(E g^2), and l2 l3 (dualtail.py's two losses).""" d = r - 1 - s c = exact_cstar(r, J, s) P = pmul(Q(J, s), c) p = pmul([comb(s, i) * (-1) ** i for i in range(s + 1)], c) R = pmul(Q(J, s - 1), c) eps, gam = peel(c, d) a, b = mp.mpf(s) / 2, mp.mpf(J + 1 + eps) / 2 Eg = Eabs_g(gam, a, b) M = moments(Fraction(s, 2), Fraction(J + 1 + eps, 2), 2 * len(gam)) Eg2 = sum(gam[i] * gam[j] * (-4) ** (i + j) * M[i + j] for i in range(len(gam)) for j in range(len(gam))) fr = lambda q: mp.mpf(q.numerator) / q.denominator # noqa: E731 mean = I(s, J + eps) * Eg l1P = sum(abs(v) for v in P) b_ = [(-1) ** m * sum(comb(k, m) * c[k] for k in range(m, d + 1)) for m in range(d + 1)] Lh = [2 * (s + m) * I(s + m, J) for m in range(d + 1)] l2l3 = sum(abs(fr(b_[m])) * Lh[m] for m in range(d + 1)) / fr(l1P) return dict(V=var(P), varp=var(p), l1P=l1P, maxR=max(abs(v) for v in R), mean=mean, Eg=Eg, Eg2=fr(Eg2), gam=gam, eps=eps, R=R, l2l3=l2l3, c=c) # ------------------------------------------------------------ controls def v1(two=2): off = n = 0 for r in range(5, 9): for s in range(1, r - 1): for J in (r + 7, 60, 150): x = exact_read(r, s, J) d = r - 1 - s eps, e = x["eps"], d // 2 back = [Fraction(0)] * (2 * e + 1) for i, g in enumerate(x["gam"]): for j in range(2 * i + 1): back[e - i + j] += g * comb(2 * i, j) * (-1) ** j if eps: back = pmul(back, [1, 1]) bad = back != x["c"] bad |= not x["V"] <= x["varp"] <= r - 1 bad |= x["l1P"] > two * x["V"] * x["maxR"] mR = x["maxR"] bad |= mp.mpf(mR.numerator) / mR.denominator > x["mean"] * ( 1 + mp.mpf(10) ** -25) if J <= 60: R = [mp.mpf(v.numerator) / v.denominator for v in x["R"]] def mod(th): z = mp.expj(th) return abs(mp.polyval(R[::-1], z)) # split at the kinks: g(U) = 0, U = sin^2(theta/2) th = [2 * mp.asin(mp.sqrt(u)) for u in g_roots(x["gam"])[1]] pts = sorted(set(mp.linspace(0, 2 * mp.pi, 9)) | set(th) | {2 * mp.pi - v for v in th}) q = mp.quad(mod, pts) / (2 * mp.pi) bad |= abs(q / x["mean"] - 1) > mp.mpf(10) ** -12 n += 1 off += bool(bad) return n, off def v2(cache, beta=BETA): off = n = 0 worst = 0.0 for r in range(5, 15): for s in range(1, r - 1): info = cache[(r, s)] d, eps = r - 1 - s, info["eps"] K = K_up(r, s, eps, r - 1, beta) for J in (r + 7, 100, 400, 1000): m_s, _ = bounds(r, J, s) Sv = info["S"].eval(J) rho = Fraction(str(info["detN"].eval(J) / (info["D"].eval(J) * Sv))) off += rho != m_s / comb(2 * J, J) # the true side from the exact c*, peeled; the bound's E g^2 # from vtail's assembled polynomials, equal exactly _, gam = peel(exact_cstar(r, J, s), d) M = moments(Fraction(s, 2), Fraction(J + 1 + eps, 2), 2 * len(gam)) Eg2 = sum(gam[i] * gam[j] * (-4) ** (i + j) * M[i + j] for i in range(len(gam)) for j in range(len(gam))) Sq = Fraction(int(Sv.p), int(Sv.q)) if hasattr(Sv, "p") else Fraction(Sv) asm = ev2(info["Eg"], J) / (Sq ** 2 * ev2(info["Dm"], J)) off += asm != Eg2 Eg2 = float(Eg2) ub = float(K.subs(t, sp.sqrt(J))) ** 2 * float(asm) lC = (2 * J * log(2) + lgamma(J + 0.5) - lgamma(J + 1) - 0.5 * log(pi)) lIr = ((J + r - 1) * log(2) - log(pi) + lgamma((J + 1) / 2) + lgamma(r / 2) - lgamma((J + r + 1) / 2)) Je = J + eps lIs = ((Je + s - 1) * log(2) - log(pi) + lgamma((Je + 1) / 2) + lgamma(s / 2) - lgamma((Je + s + 1) / 2)) from math import exp tr = (2 * (r - 1) * exp(lIr + lIs - lC)) ** 2 * Eg2 n += 1 off += ub < tr * (1 - 1e-12) worst = max(worst, tr / ub) return n, off, worst def v3(Jt): bad = n = 0 for (r, s), jt in Jt.items(): for J in range(jt, jt + 41): n += 1 bad += not bounds(r, J, s)[1] > champ_height(r, J) return n, bad def carry(Jt): """The per-depth certificate DUAL_s > H_ch from 501 to J_t(r, s) - 1 at each order of ranks 11..14 whose tail starts past 500.""" out = {} for r in range(11, 15): fails, n, H = [], 0, {} for s in range(r - 1): for J in range(501, Jt[(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 from 501: %d fail %s" % (r, n, len(fails), fails[:5]), flush=True) return out def main(): t0 = time.time() new, old, cache = {}, {}, {} for r in range(5, 15): for s in range(r - 1): if s >= 1: new[(r, s)], cache[(r, s)] = vtail(r, s) if r >= 11 or s >= r - 3 or s == 0: old[(r, s)] = tri_tail(r, s)[0] print("tails built (%.1f s)" % (time.time() - t0), flush=True) n1, o1 = v1() n1p, o1p = v1(two=1) print("V1 sign changes, variation, Beta form, quadrature and peel, r = 5..8," " s = 1..r - 2, J in {r + 7, 60, 150}: %d cells, %d off; PLANT V for" " 2V: %d off" % (n1, o1, o1p), flush=True) n2, o2, w2 = v2(cache) n2p, o2p, _ = v2(cache, beta=sp.Rational(7, 10)) print("V2 the rational chain against lgamma and rho exactly, r = 5..14," " s >= 1, J in {r + 7, 100, 400, 1000}: %d cells, %d off, true/bound" " at most %.6f; PLANT beta = 0.7: %d off" % (n2, o2, w2, o2p), flush=True) Jt = {} for r in range(5, 15): for s in range(r - 1): v = [x for x in (new.get((r, s)), old.get((r, s))) if x is not None] Jt[(r, s)] = min(v) if v else None opened = [k for k, v in list(new.items()) + list(old.items()) if v is None] if any(v is None for v in Jt.values()) or opened: print("an order is open: %s" % opened) return True n3, o3 = v3(Jt) print("V3 DUAL_s > H_ch at the 41 depths from each order's taken start" " (the smaller of the two bounds), r = 5..14, every s: %d checked," " %d off (%.1f s)" % (n3, o3, time.time() - t0), flush=True) ctl = o1 == 0 and o1p > 0 and o2 == 0 and o2p > 0 and o3 == 0 print("controls: %s" % ("all ok" if ctl else "OFF"), flush=True) print("PV1 at s = r - 3, dualtail's J_t: V = var(P), the loss 2(r - 1)" " mean|R|/||P||_1, the Cauchy-Schwarz factor, l2 l3") pv1 = True for r in range(5, 11): x = exact_read(r, r - 3, JT[r]) loss = 2 * (r - 1) * x["mean"] / x["l1P"] cs = mp.sqrt(x["Eg2"]) / x["Eg"] pv1 &= 2.5 <= loss <= 5 and 1 <= cs <= 1.3 and ( r == 5 or loss < x["l2l3"]) print(" r = %2d J %4d var(P) %d loss %.4f CS %.4f loss x CS %.4f" " l2 l3 %.4f" % (r, JT[r], x["V"], loss, cs, loss * cs, x["l2l3"]), flush=True) print("TAIL: J_t(r, s), the smaller of the whole-norm and triangle starts") Jr = {} for r in range(5, 15): row = [Jt[(r, s)] for s in range(r - 1)] Jr[r] = max(x if x is not None else 10 ** 9 for x in row) bind = max(range(r - 1), key=lambda s: row[s] or 10 ** 9) print(" r = %2d J_t = %5d binding s = r - %d by s %s whole-norm" " %s" % (r, Jr[r], r - bind, row, [new.get((r, s)) for s in range(r - 1)]), flush=True) pv2 = all(new[(r, r - 3)] <= 0.5 * old[(r, r - 3)] for r in range(7, 11)) pv3 = all(Jt[(r, s)] is not None for r in range(11, 15) for s in range(r - 1)) and all( 700 <= Jr[r] <= 2500 for r in range(11, 15)) and all( max(range(r - 1), key=lambda s: Jt[(r, s)]) >= r - 3 for r in range(11, 15)) pv4 = all(new[(r, s)] < old[(r, s)] for r in range(11, 15) for s in range(1, r - 2)) kill = [r for r in range(7, 11) if new[(r, r - 3)] >= old[(r, r - 3)]] + [ r for r in range(11, 15) if Jr[r] > 5000] print("PV1 %s; PV2 r - 3 at most half dualtail's, ranks 7..10: %s; PV3" " ranks 11..14 close in [700, 2500] at r - 3 or r - 2: %s; PV4 below" " the triangle at 1 <= s <= r - 3, ranks 11..14: %s" % tuple("holds" if v else "fails" for v in (pv1, pv2, pv3, pv4))) print("KILL %s %s" % (kill, "fires" if kill else "misses"), flush=True) print("CARRY: the certificate from 501 to J_t(r, s) - 1, ranks 11..14 " "(%.1f s)" % (time.time() - t0), flush=True) cr = carry(Jt) held = all(not cr[r][1] for r in cr) print("carried: no (order, depth) cell fails: %s [%.1f s]" % (held, time.time() - t0)) return not (ctl and held) or bool(kill) if __name__ == "__main__": sys.exit(1 if main() else 0)