#!/usr/bin/env python3 """Theorem 9.24 (ii) and Theorem 9.27 of the paper "Analytic q-inequalities and a proof of Bergeron's conjecture": R(x, A) for A >= lam0 x^2. By Lemma 9.20 (closed-form evaluation bound), N_A(rho)/P_A(rho) <= (x/D)^2 (1 + W e^{1+1/lam}), W = lam x/(x-1), which is <= 0.187 for lam = A/x^2 >= 16 and x >= 101, and decreases in sigma = x s. By Lemma 9.19 it then suffices to bound pl <= sqrt(2 pi k2) pp and Lu >= sqrt(2 pi k2) Lambda for all Gaussian coefficients G(m, n) = [m+n choose m]_q with m >= mc = xc - 2, n >= (lam0 - 2/xc) m and every tilt that occurs, uniformly. By Lemma 9.19 (b) only degrees k >= D + x need the criterion, and there s < pi/sqrt(6 D) <= pi/sqrt(6 (lam0 x - 1)); with --xmax (only x <= XMAX) this bounds sigma. Regimes (tau = n s, sigma = m s, mu = m/n, E(t) = (1 - e^{-t})/t): N tau <= 8, n <= C0 m^2/E(tau): n-scaled (vartheta' = n theta); Taylor zone with the h-form bounds, zone A (theta <= pi/m) by closed-form pair factors, zone B (theta >= pi/m) by the Bernoulli bound of Lemma 9.21. C tau <= 8, U(rho) = [n+1]_rho >= C0 m^2: Lemma 9.23 (cycle index) and the i.i.d. sum U^m (family N with the single pair (1, n+1), no zone B). Sm SIG_SS <= sigma <= 6 (or the bound from --xmax): the m-scaled table of pn_certificate.py with n = oo allowed. Ss sigma <= SIG_SS, tau >= 8 (sums of nearly exponential variables): s-scaled; zone A (theta <= pi/m) by a Riemann integral in t = j/m, zone B (theta >= pi/m) by Lemma 9.22 (refined for each m <= XMAX if given). P sigma >= 6: regime P of pn_certificate.py (Section 9.4, value version). Every inequality is checked in ball arithmetic (python-flint) and counts only if it holds on the whole cell; each regime is checked at m = mc together with the monotonicity in m of its error terms. Runs used in the paper: --xc 101 (Theorem 9.24 (ii)); --xc 22 --xmax 40 --lam0 64 --sig-ss 0.25, --xc 41 --xmax 60 --lam0 32 and --xc 61 --xmax 100 --lam0 16 (Theorem 9.27). """ import sys, time, argparse, math from flint import arb, ctx import cww_certificate as W from cww_certificate import A, up, lo, require, CertError, PI, ZETA2, AR, H, sinhc, sinc, floor3, Psibar from rq_certificate import gauss_tail, lam_bounds import pn_certificate as PN ctx.prec = 96 SQ2PI = (2 * PI).sqrt() SIN8 = (PI / 8).sin() H40 = H(4, 0) # = 2/15 = sup_{y real} |h''''(y)| C0 = 20 # family C needs U(rho) >= C0 m^2 def dn(v, d=4): """lower bound of the ball v, rounded down to d decimals (for the log)""" return "%.*f" % (d, math.floor(float(A(v).lower()) * 10**d) / 10**d) def upd(v, d=4): """upper bound of the ball v, rounded up to d decimals (for the log)""" return "%.*f" % (d, math.ceil(float(A(v).upper()) * 10**d) / 10**d) def ups(v, g=3): """upper bound of the ball v, rounded up to g significant digits (for the log)""" u = float(A(v).upper()) if u <= 0: return "%.*g" % (g, u) e = math.floor(math.log10(u)) - (g - 1) return "%.*g" % (g, math.ceil(u / 10**e) * 10**e) def amin(*xs): b = None for x in xs: if x is None: continue x = A(x) if b is None or x.upper() < b.upper(): b = x return b def local_generic(e4, y3, y1, uZ, zc, ucell=0.02, Ucap=12, umin=1.5): """(Ilow, Iup, U) with u = theta sqrt(kappa_2) (true variance). Hypotheses on the Taylor zone |theta| <= theta0: log phi~ = -u^2/2 - i kappa_3 theta^3/6 + R4, |R4| <= K4 theta^4/24; zc >= K4 theta0^2/(24 kappa_2) (so that Re log phi~ <= -u^2/4 on the zone); for |u| <= U: K4 theta^4/24 <= e4 u^4, |kappa_3 theta^3/6| <= y3 u^3, and the offset phase |delta theta| <= y1 u; U <= uZ <= theta0 sqrt(kappa_2).""" require(A(zc).upper() < 0.25, "Taylor-zone condition") U = min(Ucap, floor3(float(A(uZ).lower()))) require(U >= umin and U > 0 and arb(U) <= uZ, "U") # umin < 1.5 only where just the upper bound Iu is used ncell = int(U / ucell) + 1 du = arb(U) / ncell Il = arb(0); Iu = arb(0) for i in range(ncell): u = (du * i).union(du * (i + 1)) E4 = e4 * u**4 Y = y3 * u**3 + y1 * u + E4 ex = (-u**2 / 2 - E4).union(-u**2 / 2 + E4) cy = (1 - Y**2 / 2).union(arb(1)) Il += du * lo(ex.exp() * cy) Iu += du * up((-u**2 / 2 + E4).exp()) tail = gauss_tail(arb(U), 0) return lo(2 * (Il - tail) / SQ2PI), up(2 * (Iu + tail) / SQ2PI), U # ============================================================================================= # Lemma 9.20 (closed-form evaluation bound) # ============================================================================================= def np_crude(lam, xc, sig_lo=None): """sup over x >= xc, psi, and s with x s >= sig_lo (if given) of N/P, for A = lam x^2 (lam >= 2): N/P <= Q S / B, B >= 1; Q <= (sinh(x s/2)/sinh(D s/2))^2 <= (x/D)^2; S <= 1 + t Z^{t+1}, t <= W = lam x/(x-1), (t+1) log Z <= (t+1)(x-2)/(D+2) <= 1 + 1/lam.""" lam = A(lam); x = arb(xc) Dx = lam - 1 / x # D/x >= lam - 1/x Wt = lam * x / (x - 1) S = 1 + Wt * (1 + 1 / lam).exp() Q = 1 / Dx**2 if sig_lo is not None: # sinh(sigma/2)/sinh(Dx sigma/2) decreases in sigma (Dx > 1) sg = A(sig_lo) Q = min(up(Q), up(((sg / 2).sinh() / (Dx * sg / 2).sinh())**2)) return up(Q * S) # ============================================================================================= # helper functions # ============================================================================================= def gfun(y): """g(y) = h''(y) = 1/y^2 - 1/sinh^2(y) = Psi(y)/y^2: decreasing on y >= 0, g(0) = 1/3 (point evaluation).""" y = A(y) if y.upper() < 0.5: s = arb(0) for k in range(1, 30): s += W._q(W.a_k[k] * 2 * k * (2 * k - 1)) * y**(2 * k - 2) q = (up(abs(y)) / PI)**2 tail = arb('8.8') / PI**2 * 30 * q**29 / (1 - q)**2 return s + arb(0, 1) * tail return 1 / y**2 - 1 / y.sinh()**2 def Efun(t): """E(t) = (1 - e^{-t})/t (E(0) = 1), decreasing; point evaluation.""" return PN._Epoint(A(t)) def Bfun(y): """B(y) = y/(e^y - 1) (B(0) = 1), decreasing, |B'| <= 1/2.""" y = A(y) if y.upper() == 0: return arb(1) return y / y.expm1() def psi2(y): """(y/2)^2/sinh^2(y/2) = 1/sinhc(y/2)^2, decreasing from 1.""" return 1 / sinhc(A(y) / 2)**2 def N4(ta): """sup_{y >= ta} y^4 Li_{-3}(e^{-y}) for ta >= 4 (y^4 e^{-y} and (1+4e^{-y}+e^{-2y})/(1-e^{-y})^4 decrease).""" y = A(ta) require(y.lower() >= 4, "N4: ta >= 4") e = (-y).exp() return up(y**4 * e * (1 + 4 * e + e * e) / (1 - e)**4) # ============================================================================================= # Family N (and the i.i.d. part of family C): n-scaled. kappa_2 = m n^2 b2, theta = vartheta'/n. # ============================================================================================= def b2_bounds_N(ta, tb, mub): """kappa_2/(m n^2) in [b2m, b2p] for tau in [ta, tb], every j/n <= mub: kappa_2 = s^-2 sum_j [Psi((n+j)s/2) - Psi(js/2)] with Psi(y) = y^2 g(y), g decreasing, g <= 1/3.""" b2m = (gfun(A(tb) / 2) - A(mub)**2 / 3) / 4 b2p = (1 + A(mub))**2 / 4 * gfun(A(ta) / 2) require(lo(b2m) > 0, "b2 lower bound") return lo(b2m), up(b2p) def deriv_consts_N(ta, tb, mub, th0): """(k4n, c3n): sup_{|theta| <= th0/n} |L''''(s - i theta)| <= m n^4 k4n and |kappa_3| <= m n^3 c3n. h-form L'''' = sum_j [((n+j)/2)^4 h''''((n+j)w/2) - (j/2)^4 h''''(jw/2)], |h''''| <= H(4, v) on |Im| <= v < pi; numerators alternatively <= 6(1 + A_3)/s^4 each (decayed form). |h'''(y)| <= min(H(4,0) y, 2/y^3), y >= 0.""" ta = A(ta); tb = A(tb); mub = A(mub); th0 = A(th0) v = (1 + mub) * th0 / 2 vd = mub * th0 / 2 require(vd.upper() < float(PI.lower()) - 0.05, "N: denominator strip") den = mub**4 * H(4, vd) / 16 numc = [] if v.upper() < float(PI.lower()) - 0.05: numc.append(up((1 + mub)**4 * H(4, v) / 16)) if ta.lower() > 0: numc.append(up(6 * (1 + AR[3]) / ta**4)) require(numc, "N: no numerator bound") k4n = up(amin(*numc) + den) y1 = ta / 2; y2 = (1 + mub) * tb / 2 hn = up(H40 * y2) if y1.lower() > 0: hn = amin(hn, up(2 / y1**3)) c3n = up((1 + mub)**3 / 8 * hn + mub**3 / 8 * H40 * mub * tb / 2) return k4n, c3n def FN_upper(ta, tb, tha, thb, mub): """sup of F(u) = log(tau^2 H1^2 + th^2 S1^2) - log(tau^2 H0^2 + th^2 S0^2) + 2 log(H0/H1) (= log of the pair factor |1-rho^{n+j}e^{i(n+j)theta}|^2/(1-rho^{n+j})^2 * (1-rho^j)^2/|1-rho^j e^{ij theta}|^2, u = j/n) over tau in [ta, tb], th = n theta in [tha, thb], u in [0, mub], in zone A (u th <= pi): S1 = |sinc((1+u)th/2)|, H1 = sinhc((1+u)tau/2), S0 = sinc(u th/2) >= sinc(min(mub thb/2, pi/2)), H0 = sinhc(u tau/2).""" ta, tb, tha, thb, mub = map(A, (ta, tb, tha, thb, mub)) xa = tha / 2; xb = (1 + mub) * thb / 2 S1M = up(abs((xa.union(xb)).sinc())) S1M = amin(S1M, arb(1), up(1 / xa)) H1m = lo(sinhc(ta / 2)); H1M = up(sinhc((1 + mub) * tb / 2)) z0 = mub * thb / 2 if z0.upper() > float((PI / 2).lower()): z0 = PI / 2 S0m = lo(sinc(z0)) H0M = up(sinhc(mub * tb / 2)) # with rho = tau/th in [ra, rb]: F = log(rho^2 H1^2 + S1^2) - log(rho^2 H0^2 + S0^2) + 2 log(H0/H1) # = log(1 + S1^2/(rho^2 H1^2)) - log(1 + S0^2/(rho^2 H0^2)) ra = ta / thb; rb = tb / tha c1 = (rb**2 * H1M**2 + S1M**2).log() - (ra**2 + S0m**2).log() + 2 * (H0M / H1m).log() best = up(c1) if ra > 0: c2 = up((1 + S1M**2 / (ra**2 * H1m**2)).log() - (1 + S0m**2 / (rb**2 * H0M**2)).log()) if c2 < best: best = c2 require(best.is_finite(), "FN not finite") return best def zoneA_N(ta, tb, mub, th0, mc, b2p, ratio='1.06', thmax=None): """F >= sqrt(kappa_2/(2pi)) * 2 int_{zone A} |phi| for every m >= mc; zone A = [th0, pi/mu] in vartheta' (integrated to thmax or to oo with the tail bound). log|phi|^2 <= m sup F on each cell.""" th = A(th0); tot = arb(0); worst = None r = arb(ratio) ta = A(ta); tb = A(tb) H1m = lo(sinhc(ta / 2)); H1M = up(sinhc((1 + A(mub)) * tb / 2)); H0M = up(sinhc(A(mub) * tb / 2)) K = up((tb**2 * H1M**2 + 4) * H0M**2 / H1m**2) # tail: exp(F) <= K (pi/(2 th))^2 thT = max(60.0, 3.0 * float(K.sqrt().upper()) * math.pi / 2) if thmax is not None: thT = min(thT, float(A(thmax).upper())) while th.upper() < thT: tb_ = th * r if tb_.upper() > thT: tb_ = arb(thT) F = FN_upper(ta, tb, th, tb_, mub) require(F < 0 and (-F * mc).lower() > 1.1, "N zone A: rate (%s at th=%.3f)" % (F.str(3), float(th.lower()))) worst = F if worst is None or F > worst else worst tot += (tb_ - th) * (mc * F / 2).exp() th = tb_ if thmax is None or thT < float(A(thmax).lower()): # tail th >= thT: exp(m F/2) <= K^{m/2} (pi/(2 th))^m; int = thT/(m-1) (sqrt(K) pi/(2 thT))^m rr = K.sqrt() * PI / (2 * thT) require(rr.upper() < 0.5, "N zone A tail ratio") tot += arb(thT) / (mc - 1) * rr**mc return up((mc * b2p / (2 * PI)).sqrt() * 2 * tot), worst def zoneB_N(ta, tb, mub, mc, b2p, nmax_kind, nmax_coef): """F >= sqrt(kappa_2/(2pi)) * 2 int_{pi/m}^{pi} |phi| for every m >= mc (Bernoulli-Abel): |phi| <= (4/3)(4 sqrt(m) + 1) exp(-m(1 - P0) + pimax/sin(pi/(2m))), p_0 <= mu/E(tau) <= P0 = mub/E(tb), V >= m(1 - P0); kappa_2 <= n^2 m b2p with n <= nmax = coef m^p (p = 1: n = m/mu <= m/mua; p = 2: n <= C1 m^2).""" P0 = up(A(mub) / Efun(tb)) require(P0 < 1, "N zone B: P0 < 1") pim = up(P0 * (1 - P0)) if P0 <= 0.5 else arb(1) / 4 x = arb(mc) e = -(1 - P0) + pim / (x * (PI / (2 * x)).sin()) # exponent per m, decreasing in m require(e < 0, "N zone B exponent") p = nmax_kind require((-e * mc).lower() > p + 1.5, "N zone B: m-monotonicity") F = arb(nmax_coef) * x**p * (x * b2p / (2 * PI)).sqrt() * 2 * arb(4) / 3 * (4 * x.sqrt() + 1) * PI * (x * e).exp() return up(F), e def table_cell_N(ta, tb, mua, mub, mc, side, iid=False, y1_override=None, cands=None, umin=1.5): """(pl, Lu, info) for G(m, n), tau in [ta, tb], mu = m/n in [mua, mub] (iid: the i.i.d. sum of m truncated geometrics on {0..n}, 1/n <= mub), every m >= mc; side 'P' includes the offset of G(x, D) (Delta <= m/2 and <= B(tau)/s), side 'N' has no offset.""" ta = A(ta); tb = A(tb); mua = A(mua); mub = A(mub) b2m, b2p = b2_bounds_N(ta, tb, mub) if y1_override is not None: y1 = A(y1_override) elif side == 'P': y1 = up(mub / (2 * (mc * b2m).sqrt())) if ta.lower() > 0: y1 = amin(y1, up(Bfun(ta) / (ta * (mc * b2m).sqrt()))) else: y1 = arb(0) if not iid: if mua.lower() > 0: nk, nc = 1, up(1 / mua) else: nk, nc = 2, up(C0 / Efun(tb)) FB, eB = zoneB_N(ta, tb, mub, mc, b2p, nk, nc) else: FB, eB = arb(0), None best = None; why = {} if cands is None: cands = ['0.8', '1.0', '1.2', '1.4', '1.6', '1.8', '2.0', '2.2', '2.4', '2.7', '3.0', '3.5', '4.0', '5.0'] for c in cands: th0 = arb(c) try: k4n, c3n = deriv_consts_N(ta, tb, mub, th0) zc = k4n * th0**2 / (24 * b2m) e4 = k4n / (24 * mc * b2m**2) y3 = c3n / (6 * arb(mc).sqrt() * b2m * b2m.sqrt()) uZ = th0 * (mc * b2m).sqrt() Il, Iu, U = local_generic(e4, y3, y1, uZ, zc, umin=umin) FA, worstA = zoneA_N(ta, tb, mub, th0, mc, b2p) except CertError as ex: why[c] = str(ex) continue pl = lo(Il - FA - FB); Lu = up(Iu + FA + FB) sc = float(pl.lower()) / float(Lu.upper()) if best is None or sc > best[0]: best = (sc, c, pl, Lu, dict(Il=float(Il.lower()), Iu=float(Iu.upper()), FA=float(FA.upper()), FB=float(FB.upper()), U=U, th0=c, worstA=float(worstA.upper()))) require(best is not None, "N: no admissible Taylor zone: %s" % why) return best[2], best[3], best[4] # ============================================================================================= # Family C: cycle index. G(m,n)(q) = (1/m!) sum_{sigma in S_m} prod_{cycles c} U(q^{|c|}), U = [n+1]_q. # ============================================================================================= def cell_C(ta, tb, mc, side, c0=C0): """(pl, Lu, info) for G(m, n) with tau in [ta, tb], U(rho) >= c0 m^2 (guaranteed by n >= c0 m^2/E(tb)), every m >= mc. pp_G >= pp_{U^m} / prod(1 + i/U); Lambda_G <= sum_sigma U^{cyc - m} Lambda_{f(sigma)}.""" ta = A(ta); tb = A(tb) x = arb(mc) nmin = arb(c0) * x**2 / Efun(tb) eta = up(1 / nmin) # 1/n <= eta b2m, b2p = b2_bounds_N(ta, tb, eta) # offset of j = mean(G(x, D-1)) from mean(U^m) (and from mean(U_D^m) for G(x, D)): <= m(m-1)/4 + m/2 if side == 'P': y1 = up(Efun(tb) * (1 + 1 / x) / (4 * c0 * (x * b2m).sqrt())) else: y1 = arb(0) plI, LuI, infI = table_cell_N(ta, tb, 0, eta, mc, side, iid=True, y1_override=y1) mh = (mc + 1) // 2 _, LuH, _ = table_cell_N(ta, tb, 0, eta, mh, 'N', iid=True, umin=0.25) # only the upper bound Lu is used # variance ratio kappa_2(G)/kappa_2(U^m): mu = m/n <= E(tb)/(c0 m) mu = up(Efun(tb) / (c0 * x)) gb = gfun(tb / 2) corr = 1 / (1 - 1 / (3 * nmin**2 * gb)) rlo = lo(1 - mu**2 / (3 * gb) * corr) rhi = up(1 + mu * (2 + mu) * corr) etaC = up((1 / (2 * arb(c0))).exp() - 1) # prod(1 + i/U) - 1 <= e^{m^2/(2U)} - 1 # k > m/2 non-fixed points: weight <= sum_{k > m/2} (m/sqrt U)^k, times sqrt(2 pi kappa_2(U^m)) <= ... (uniform in n): # <= sqrt(2 pi m b2p) m^2/E c0^{1 - k0/2}/(1 - c0^{-1/2}), k0 = floor(m/2) + 1 >= (m+1)/2; with (m+1)/4 in place of # k0/2 the bound m^{5/2} c0^{-(m+1)/4} decreases in m once m > 10/log(c0), so its value at mc covers every m >= mc require(arb(10) / arb(c0).log() < mc, "C tail monotonicity") ek = arb(mc + 1) / 4 tailC = (2 * PI * x * b2p).sqrt() * x**2 / Efun(tb) * arb(c0)**(1 - ek) / (1 - 1 / arb(c0).sqrt()) pl = lo(rlo.sqrt() * plI / (1 + etaC)) Lu = up(rhi.sqrt() * (LuI + etaC * arb(2).sqrt() * LuH + tailC)) return pl, Lu, dict(plI=float(plI.lower()), LuI=float(LuI.upper()), LuH=float(LuH.upper()), tailC=float(tailC.upper()), etaC=float(etaC.upper()), th0=infI['th0']) # ============================================================================================= # Family Sm: m-scaled, nu in [na, oo]. # ============================================================================================= def a2_bounds_inf(sa, sb, na, eps): sa = lo(sa); sb = up(sb); na = lo(na); e1 = up(eps) lo_ = 2 * (Psibar(sa * (na + 1) / 2) - Psibar(sa * na / 2) - (Psibar(sb * (1 + e1) / 2) - Psibar(sb * e1 / 2))) / sb**3 hi_ = (sb - 2 * Psibar(sa / 2)) / sa**3 # sum_j Psi((n+j)s/2) <= m require(lo_.is_finite() and hi_.is_finite() and lo(lo_) > 0, "a2 inf bracket") return lo(lo_), up(hi_) def deriv_consts_inf(sa, sb, na, th0, eps): """(k4z, c3): as deriv_consts_G with n = oo allowed (tau >= na sa).""" sa = lo(sa); th0 = A(th0); eps = A(eps) k4c = [] vd = th0 / 2 if vd.upper() < float(PI.lower()) - 0.05: k4c.append(up(H(4, vd) * (1 + up(eps))**5 / 80 + 6 * (1 + AR[3]) / sa**4 * (1 + up(eps)))) ta = na * sa if ta.lower() >= 4: k4c.append(up((6 * AR[3] + N4(ta)) / sa**4)) # l-form: j^4 Li_{-3}(e^{-js}) <= 6 A_3/s^4 if sa > 0.5: k4c.append(up(AR[4] * 24 * ZETA2 / sa**5)) require(k4c, "inf: no fourth-derivative bound") c3c = [up(arb(2) / sa**4)] if sa > 0.5: c3c.append(up(AR[3] * 6 * ZETA2 / sa**5)) return amin(*k4c), amin(*c3c) def F_upper_inf(ta, tb, sa, sb, tha, thb, na): """F(t) <= 2 log coth((na + ta) sa/2) - log(1 + X0), X0 = (vt/sig)^2 S0^2/H0^2.""" ta, tb, sa, sb, tha, thb, na = map(A, (ta, tb, sa, sb, tha, thb, na)) S0m, S0M = PN._absrange(lambda z: z.sinc(), ta * tha / 2, tb * thb / 2) H0m, H0M = PN._absrange(sinhc, ta * sa / 2, tb * sb / 2) require(S0m > 0, "sinc vanishes") y1 = (na + ta) * sa / 2 X0m = (tha / sb)**2 * S0m**2 / H0M**2 return up(2 * (1 / y1.tanh()).log() - (1 + X0m).log()) def pair_cell_inf(sa, sb, na, tha, thb): NX = PN.NX ts = [arb(i) / NX for i in range(NX + 1)] M = [F_upper_inf(ts[c], ts[c + 1], sa, sb, tha, thb, na) for c in range(NX)] require(all(m.is_finite() for m in M), "cell bounds not finite") Aval = sum(M, arb(0)) / NX PV = sum(((M[c] - M[c - 1]).nonnegative_part() for c in range(1, NX)), arb(0)) return up(Aval), up(PV) def profile_inf(sa, sb, na, mc): """(v, pm): V >= v m, max p(1-p) <= pm, for nu >= na (Fermi level lam in [log(e^sa - 1), lam+(sb, na)]).""" sa = lo(sa); sb = up(sb); na = lo(na) lam_lo = lo(sa.expm1().log()) lam_hi = up(lam_bounds(sb, na)) lam = lam_lo.union(lam_hi) eps_all = arb(0).union(arb(1) / mc) sig = sa.union(sb) a1 = lam - sig * eps_all a2 = sa * (1 + na) - sb / mc - lam_hi # <= s(M'-1) - s c for every nu >= na require(a2.lower() > 0, "profile a2") t1 = lo((a1 / 2).tanh()); t2 = lo((a2 / 2).tanh()) v = lo((t1 + t2) / (2 * sb) - arb(1) / (4 * mc)) require(v.is_finite() and v > 0, "profile v") pm = arb(1) / 4 if lam_hi >= 0 else up((1 - (lam_hi / 2).tanh()**2) / 4) return v, pm def far_zone_inf(sa, sb, na, th0, mc, a2p): sa = A(sa); sb = A(sb); na = A(na) x = arb(mc) v, pm = profile_inf(sa, sb, na, mc) pref = arb(4) / 3 * (4 * x.sqrt() + 1) # V <= m tot = arb(0); rateA = None; rateB = None th = A(th0) e_sw = max(0.05, 4.5 / mc) # Bernoulli blocks then have rate * mc > 4 (checked) while True: e_here = PN.abel_expo(v, pm, th / x, x) / x if e_here.upper() < -e_sw or th.upper() >= 8: break thb = th * arb('1.04') Aa, Bb = pair_cell_inf(sa, sb, na, th, thb) require(Aa.is_finite() and Bb.is_finite(), "pair cell not finite") rateA = Aa if rateA is None or Aa > rateA else rateA tot += (thb - th) / x * ((x * Aa + Bb) / 2).exp() th = thb edges = [] t = th while t.upper() < 60: tb = t * arb('1.1') edges.append((t / x, tb / x)); t = tb tt = t / x # for m > mc the scaled blocks end at theta = t/m < t/mc: the gap [t/m, t/mc] has length <= t/mc and there # -V + |Z| <= m (-v + pm/(m sin(t/(2m)))) <= m e(mc), e(mc) = abel_expo(v, pm, t/mc, mc)/mc; its share of F, # ~ m^2 e^{m e(mc)}, decreases in m (rateB mc > 4), so its value at mc is added gap_expo = PN.abel_expo(v, pm, tt, x) require(gap_expo.is_finite(), "gap exponent") tot += tt * pref * gap_expo.exp() while tt.upper() < float(PI.lower()): tb = tt * arb('1.2') if tb.upper() > float(PI.lower()): tb = PI edges.append((tt, tb)); tt = tb for (a_, b_) in edges: expo = PN.abel_expo(v, pm, a_, x) require(expo.is_finite(), "Abel exponent not finite") tot += (b_ - a_) * pref * expo.exp() r = lo(-expo / x) rateB = r if rateB is None or r < rateB else rateB F = (a2p * x**3 / (2 * PI)).sqrt() * 2 * tot return up(F), rateA, rateB def table_cell_inf(sa, sb, na, mc, side): sa = A(sa); sb = A(sb); na = A(na) if side == 'P': doff = amin(arb(1) / 2, up(Bfun(na * sa) / sa)) else: doff = arb(0) sig = sa.union(sb) eps = arb(0).union(arb(1) / mc).nonnegative_part() a2m, a2p = a2_bounds_inf(sa, sb, na, eps) best = None; why = {} for t0 in PN.taylor_candidates(float(sa.lower()), float(sb.upper()), 1e9)[4:]: try: k4z, c3 = deriv_consts_inf(sa, sb, na, t0, eps) require(mc == PN.MC, "local_value uses MC") Il, Iu, U = PN.local_value(sig, a2m, k4z, c3, t0, eps, doff) F, rateA, rateB = far_zone_inf(sa, sb, na, t0, mc, a2p) require(rateA is None or rateA < -arb(2) / mc, "Riemann rate") require((rateB * mc).lower() > 4, "Bernoulli rate") except CertError as e: why[float(t0.lower())] = str(e) continue plc = lo(Il - F); Luc = up(Iu + F) score = float(plc.lower()) / float(Luc.upper()) if best is None or score > best[0]: best = (score, t0, plc, Luc, dict(Il=float(Il.lower()), Iu=float(Iu.upper()), far=float(F.upper()), U=U, th0=float(t0.lower()))) require(best is not None, "Sm: no admissible Taylor zone: %s" % why) return best[2], best[3], best[4] # ============================================================================================= # Family Ss: s-scaled, sigma in (0, sb], tau >= ta >= 4 (Gamma regime). # ============================================================================================= def good_fraction(m): """#{1 <= j <= m : |sin(j theta/2)| >= sin(pi/8)} >= 7m/16 - 5/8 for pi/m <= theta <= pi (Lemma 9.22).""" return arb(7) * m / 16 - arb(5) / 8 def zoneA_S(sb, ta, gam, mc, nt=16, ratio='1.1', cT=40): """sqrt(kappa_2/(2pi)) 2 int_{gam s}^{pi/m} |phi| <= sqrt(m/(2pi)) 2 int_gam^{pi/sigma} |phi(c s)| dc, every m >= mc, sigma in (0, sb]: log|phi|^2 <= 2m log coth(ta/2) - m I(c) + log(1 + c^2), I(c) >= int_0^1 log(1 + ca^2 sinc^2(vb t/2)/sinhc^2(sb t/2)) dt, vb = min(cb sb, pi).""" sb = up(sb); x = arb(mc) lc = (1 / (A(ta) / 2).tanh()).log() tot = arb(0); worst = None c = A(gam); r = arb(ratio) while c.upper() < cT: cb = c * r if cb.upper() > cT: cb = arb(cT) vb = cb * sb if vb.upper() > float(PI.lower()): vb = PI f = lambda t: (1 + c**2 * sinc(vb * t / 2)**2 / sinhc(sb * t / 2)**2).log() I = W.lower_riemann(f, arb(0), arb(1), nt) rate = lc - I / 2 # log|phi| <= m rate + log(1+cb^2)/2 require(rate < 0 and (-rate * mc).lower() > 0.6, "S zone A rate") worst = rate if worst is None or rate > worst else worst tot += (cb - c) * (x * rate + (1 + cb**2).log() / 2).exp() c = cb # tail c >= cT: I >= log(1 + (2c/pi)^2/sinhc^2(sb/2)) >= 2 log(2c/(pi sh)); 1 + c^2 <= 2c^2 sh = sinhc(sb / 2) K = (1 / (A(ta) / 2).tanh()) * PI * sh / 2 rr = K / cT require(rr.upper() < 0.5, "S zone A tail") tot += arb(2).sqrt() * arb(cT)**2 / (mc - 2) * rr**mc return up((x / (2 * PI)).sqrt() * 2 * tot), worst def zoneB_S(sb, ta, mc, mmax=None): """sqrt(kappa_2/(2pi)) 2 int_{pi/m}^{pi} |phi| <= sqrt(2 pi m) (m/sigma) coth(ta/2)^m (sinh(sigma/2)/sin(pi/8))^g, g >= 7m/16 - 5/8; increasing in sigma, decreasing in m (checked). With mmax (only mc <= m <= mmax): for a good j the factor is <= sinh(js/2)/sin(pi/8) <= (j/m) r (sinh is convex), so the product over g0 = ceil(7m/16 - 5/8) good j is <= r^g0 prod_{j=m-g0+1}^m (j/m); evaluated for every m.""" sb = up(sb) r = (sb / 2).sinh() / SIN8 require(r < 1, "S zone B ratio") lc = (1 / (A(ta) / 2).tanh()).log() if mmax is None: x = arb(mc) require((lc + arb(7) / 16 * r.log()).upper() < -1.6 / mc, "S zone B m-monotonicity") F = (2 * PI * x).sqrt() * (x / sb) * (x * lc).exp() * r**good_fraction(mc) return up(F) worst = arb(0) for m in range(mc, mmax + 1): g0 = -((-7 * m + 10) // 16) # ceil((7m - 10)/16) = ceil(7m/16 - 5/8) require(arb(g0) >= good_fraction(m) and g0 >= 2, "g0") prod = arb(1) for j in range(m - g0 + 1, m + 1): prod *= arb(j) / m x = arb(m) F = (2 * PI * x).sqrt() * (x / sb) * (x * lc).exp() * r**g0 * prod if F > worst: worst = up(F) return worst def cell_Ss(sb, ta, mc, side, mmax=None, cands=('0.5', '0.6', '0.7', '0.8', '0.9')): sb = up(sb); ta = lo(ta) c2m = lo(psi2(sb) - psi2(ta)) require(c2m > 0, "c2m") k4 = up(6 * AR[3] + N4(ta)) # K4 <= m k4 / s^4 (l-form, valid for all theta) x = arb(mc) y3 = up(1 / (3 * c2m * c2m.sqrt() * x.sqrt())) # |kappa_3| <= 2 m/s^3 y1 = up(Bfun(ta) / (c2m * x).sqrt()) if side == 'P' else arb(0) e4 = up(k4 / (24 * x * c2m**2)) FB = zoneB_S(sb, ta, mc, mmax) best = None; why = {} for g in cands: gam = arb(g) try: zc = k4 * gam**2 / (24 * c2m) uZ = gam * (c2m * x).sqrt() Il, Iu, U = local_generic(e4, y3, y1, uZ, zc) FA, worst = zoneA_S(sb, ta, gam, mc) except CertError as ex: why[g] = str(ex); continue pl = lo(Il - FA - FB); Lu = up(Iu + FA + FB) sc = float(pl.lower()) / float(Lu.upper()) if best is None or sc > best[0]: best = (sc, pl, Lu, dict(Il=float(Il.lower()), Iu=float(Iu.upper()), FA=float(FA.upper()), FB=float(FB.upper()), U=U, gam=g)) require(best is not None, "Ss: no admissible zone: %s" % why) return best[1], best[2], best[3] # ============================================================================================= # Driver # ============================================================================================= def _run_cells(name, fn, cells, split, max_depth=3, log=sys.stdout): """cells: list of parameter tuples; fn(*cell) -> (pl, Lu, info); failing cells are split by split(cell).""" ok = []; failed = [] todo = list(cells) for depth in range(max_depth + 1): nxt = [] for c in todo: try: pl, Lu, info = fn(*c) ok.append((c, pl, Lu, info)) except CertError as e: if depth < max_depth and split is not None: nxt += split(c) else: failed.append((c, str(e)[:200])) todo = nxt if not todo: break return ok, failed def main(): ap = argparse.ArgumentParser() ap.add_argument('--xc', type=int, default=101) ap.add_argument('--lam0', type=float, default=16.0) ap.add_argument('--sig-ss', type=float, default=0.5, help="boundary between the regimes Ss and Sm") ap.add_argument('--xmax', type=int, default=None, help="only xc <= x <= XMAX; then sigma = x s < pi XMAX / sqrt(6 (lam0 XMAX - 1)) (low-degree lemma)") args = ap.parse_args() xc = args.xc; mc = xc - 2 PN.MC = mc # pn_certificate.local_value uses MC t0 = time.time() # degrees D < k < D + x are covered by the low-degree lemma; for k >= D + x the tilted mean k - x of G(x, D-1) # is >= D >= lam0 x - 1, so s < pi/sqrt(6 (lam0 x - 1)) and sigma = x s < pi x/sqrt(6 (lam0 x - 1)) sig_top = None if args.xmax is not None: require(args.xmax >= xc, "xmax >= xc") sig_top = up(PI * args.xmax / (6 * (arb(args.lam0) * args.xmax - 1)).sqrt()) print("sources: " + PN.source_digest(['pn_families.py', 'pn_certificate.py', 'rq_certificate.py', 'cww_certificate.py'])) print("Theorem 9.24 (ii): R(x, A) for lam = A/x^2 >= %g, %s (regimes N, C, Sm, Ss, P); arb prec = %d" % (args.lam0, ("x >= %d" % xc) if sig_top is None else ("%d <= x <= %d, sigma < %.4f" % (xc, args.xmax, float(sig_top.upper()))), ctx.prec)) nu0 = args.lam0 - 2.0 / xc # every nu = n/m >= nu0 mu_top = 1 / arb(nu0) res = {} allok = True # ---- family N t = time.time() tau_edges = [0.25 * i for i in range(33)] # [0, 8] mu_cells = [(arb(0), mu_top / 8), (mu_top / 8, mu_top / 4), (mu_top / 4, mu_top / 2), (mu_top / 2, mu_top)] cellsN = [(arb(tau_edges[i]), arb(tau_edges[i + 1]), ma, mb, mc, side) for i in range(32) for (ma, mb) in mu_cells for side in ('P', 'N')] def splitN(c): ta, tb, ma, mb, m_, side = c tm = (ta + tb) / 2 return [(ta, tm, ma, mb, m_, side), (tm, tb, ma, mb, m_, side)] okN, failN = _run_cells("N", table_cell_N, cellsN, splitN) res['N'] = okN print("[N] tau <= 8, mu <= 1/%.3f, n <= %d m^2/E(tau): %d cells, %d failed %.0fs" % (nu0, C0, len(okN), len(failN), time.time() - t), flush=True) for f in failN[:5]: print(" failed:", f) allok &= not failN # ---- family C t = time.time() cellsC = [(arb(tau_edges[i]), arb(tau_edges[i + 1]), mc, side) for i in range(32) for side in ('P', 'N')] def splitC(c): ta, tb, m_, side = c tm = (ta + tb) / 2 return [(ta, tm, m_, side), (tm, tb, m_, side)] okC, failC = _run_cells("C", cell_C, cellsC, splitC) res['C'] = okC print("[C] tau <= 8, U(rho) >= %d m^2: %d cells, %d failed %.0fs" % (C0, len(okC), len(failC), time.time() - t), flush=True) # (N) with mu_a = 0 covers n <= C0 m^2/E(tb_N), (C) covers n >= C0 m^2/E(tb_C): for every tau both cells contain, # the C cell must end no later (tb_C <= tb_N), so that the two ranges of n meet for side in ('P', 'N'): nz = [(float(c[0].lower()), float(c[1].upper())) for (c, pl, Lu, info) in okN if c[-1] == side and c[2] == 0] cz = [(float(c[0].lower()), float(c[1].upper())) for (c, pl, Lu, info) in okC if c[-1] == side] for (a1, b1) in cz: for (a2, b2) in nz: if max(a1, a2) < min(b1, b2): require(b1 <= b2, "N/C threshold mismatch on [%g, %g] vs [%g, %g]" % (a1, b1, a2, b2)) print(" N/C split: every C cell ends no later than the N cells it meets", flush=True) for f in failC[:5]: print(" failed:", f) allok &= not failC # ---- family Sm t = time.time() sm_top = 6.0 if sig_top is None else min(6.0, float(sig_top.upper()) * (1 + 1e-12)) sig_ss = args.sig_ss sig_edges = [sig_ss] while sig_edges[-1] < sm_top: sig_edges.append(min(round(sig_edges[-1] * 1.1, 6), sm_top)) cellsS = [(arb(sig_edges[i]), arb(sig_edges[i + 1]), arb(nu0), mc, side) for i in range(len(sig_edges) - 1) for side in ('P', 'N')] def splitS(c): sa, sb, na, m_, side = c sm = (sa + sb) / 2 return [(sa, sm, na, m_, side), (sm, sb, na, m_, side)] okS, failS = _run_cells("Sm", table_cell_inf, cellsS, splitS) res['Sm'] = okS print("[Sm] %g <= sigma <= %.4g, nu >= %.3f: %d cells, %d failed %.0fs" % (sig_ss, sig_edges[-1], nu0, len(okS), len(failS), time.time() - t), flush=True) for f in failS[:5]: print(" failed:", f) allok &= not failS # ---- family Ss t = time.time() ss_top = arb(sig_ss) if sig_top is None or sig_top > sig_ss else sig_top okG, failG = _run_cells("Ss", lambda sb, ta, m_, side: cell_Ss(sb, ta, m_, side, args.xmax), [(ss_top, arb(8), mc, side) for side in ('P', 'N')], None) res['Ss'] = okG print("[Ss] sigma <= %.4g, tau >= 8: %d cells, %d failed %.0fs" % (float(ss_top.upper()), len(okG), len(failG), time.time() - t), flush=True) for f in failG[:5]: print(" failed:", f) allok &= not failG # ---- constants for sigma <= 6: all families, and the families N, C, Ss that cover sigma <= 1/2 plmin = None; Lumax = None; plmin_s = None; Lumax_s = None fams = ('N', 'C', 'Sm', 'Ss') if sm_top > sig_ss else ('N', 'C', 'Ss') for fam in fams: pls = [pl for (c, pl, Lu, info) in res[fam] if c[-1] == 'P'] Lus = [Lu for (c, pl, Lu, info) in res[fam] if c[-1] == 'N'] require(pls and Lus, "family %s has no verified cells" % fam) fpl = min(pls, key=lambda v: float(v.lower())); fLu = max(Lus, key=lambda v: float(v.upper())) print(" %-2s: min pl = %s, max Lu = %s" % (fam, dn(fpl), upd(fLu))) plmin = fpl if plmin is None or fpl < plmin else plmin Lumax = fLu if Lumax is None or fLu > Lumax else Lumax if fam != 'Sm': plmin_s = fpl if plmin_s is None or fpl < plmin_s else plmin_s Lumax_s = fLu if Lumax_s is None or fLu > Lumax_s else Lumax_s # variance loss kappa_2(G(x-2,e)) >= kappa_2(G(x,D-1)) (1 - loss), loss <= min(2/s^2, (D+x)^2/6)/kappa_2(G(x, D-1)) lossN = arb(0) for (c, pl, Lu, info) in res['N']: ta, tb, ma, mb = c[0], c[1], c[2], c[3] b2m, _ = b2_bounds_N(ta, tb, mb) l = (1 + mb * (1 + arb(1) / mc))**2 / (6 * b2m) if ta > 0: l = amin(l, 2 / (ta**2 * b2m)) l = up(l / mc) if l > lossN: lossN = l for (c, pl, Lu, info) in res['C']: # n >= C0 m^2/E(tb) >= C0 m^2: mu = m/n <= 1/(C0 mc) and (m+1)/n <= (mc+1)/(C0 mc^2) for m >= mc ta, tb = c[0], c[1] muC = up(1 / (arb(C0) * mc)) b2m, _ = b2_bounds_N(ta, tb, muC) l = up((1 + arb(mc + 1) / (arb(C0) * mc * mc))**2 / (6 * b2m) / mc) if ta > 0: l = amin(l, up(2 / (ta**2 * b2m) / mc)) if l > lossN: lossN = l for (c, pl, Lu, info) in res['Sm']: sa, sb, na = c[0], c[1], c[2] a2m, _ = a2_bounds_inf(sa, sb, na, arb(0).union(arb(1) / mc)) l = up(2 / (sa**2 * a2m) / mc) if l > lossN: lossN = l c2m = lo(psi2(ss_top) - psi2(arb(8))) l = up(2 / (c2m * mc)) if l > lossN: lossN = l require(lossN < 1, "loss") # sigma = x s <= sig_ss: the Gaussian coefficients of both sides have sigma <= sig_ss, so they lie in N, C or Ss rhs_s = lo(plmin_s * (1 - lossN).sqrt() / Lumax_s) npc = np_crude(args.lam0, xc) okA = bool(npc < rhs_s) allok &= okA print("[E] sigma <= %g: N/P <= %s (Lemma 9.20 at lam = %g) < (min pl / max Lu) sqrt(1 - loss) over N, C, Ss" " = %s / %s * %s = %s: %s (loss <= %s)" % (sig_ss, upd(npc), args.lam0, dn(plmin_s), upd(Lumax_s), dn((1 - lossN).sqrt(), 5), dn(rhs_s), "verified" if okA else "FAILED", upd(lossN, 5)), flush=True) # sig_ss <= sigma <= 6: all families; Lemma 9.20 in its sinh form decreases in sigma if sm_top > sig_ss: rhs = lo(plmin * (1 - lossN).sqrt() / Lumax) npc2 = np_crude(args.lam0, xc, sig_lo=arb(sig_ss)) okA2 = bool(npc2 < rhs) allok &= okA2 print("[E] %g <= sigma <= %.4g: N/P <= %s < (min pl / max Lu) sqrt(1 - loss) = %s / %s * %s = %s: %s" % (sig_ss, sm_top, ups(npc2), dn(plmin), upd(Lumax), dn((1 - lossN).sqrt(), 5), dn(rhs), "verified" if okA2 else "FAILED"), flush=True) if sig_top is not None and sig_top < 6: print("ALL VERIFIED" if allok else "NOT VERIFIED", "(%.0fs)" % (time.time() - t0)) return 0 if allok else 1 # ---- regime P (sigma >= 6): pl/Lu from pn_certificate's regime P, N/P <= Lemma 9.20 with sigma >= 6 t = time.time() # s < pi/sqrt(6 (lam0 x - 1)) <= pi/sqrt(6 (lam0 xc - 1)) (see sig_top above); regime P needs m s >= SIGP0 for # m = x - 2 when x s >= 6, and n s >= NUMIN SIGP0 (here nu >= nu0 >= 2) require(6 * (1 - arb(2) / xc) >= PN.SIGP0 and arb(nu0) >= PN.NUMIN, "regime P constants need a larger xc") smax = PI / (6 * (arb(args.lam0) * xc - 1)).sqrt() cells = [(arb(0), arb('0.0128'))] v = arb('0.0128') while v < smax: w = v * arb('1.15') if w > smax: w = up(smax) cells.append((v, w)); v = w worstP = None for (a, b) in cells: okc = False for g in ('0.5', '0.55', '0.6', '0.45'): try: pl, Lu, info = PN.check_P_cell(a, b, arb(g)) okc = True break except CertError: continue require(okc, "regime P cell failed") r = pl / Lu worstP = r if worstP is None or r < worstP else worstP lossP = up(2 * smax / (PN.K2P - smax - PN.Gtail(PN.SIGP0) - PN.Gtail(PN.NUMIN * PN.SIGP0))) rhsP = lo(worstP) * (1 - lossP).sqrt() npP = np_crude(args.lam0, xc, sig_lo=6) okP = bool(npP < rhsP) allok &= okP print("[P] sigma >= 6: N/P <= %s < %s: %s %.0fs" % (ups(npP), dn(rhsP), "verified" if okP else "FAILED", time.time() - t), flush=True) print("ALL VERIFIED" if allok else "NOT VERIFIED", "(%.0fs)" % (time.time() - t0)) return 0 if allok else 1 if __name__ == "__main__": sys.exit(main())