#!/usr/bin/env python3 """Certificate for Theorem 9.16 (ii) of the paper "Analytic q-inequalities and a proof of Bergeron's conjecture" (Section 9.5). Gamma_{x,y}(q) = prod_{j=2}^{x} (1 - q^{y+j}) / (1 - q^j) = sum_k gamma(k) q^k (a power series), mu = (x-1) y / 2. This is the rational q-Catalan polynomial C_{x,y+1} when gcd(x, y+1) = 1. Claim, for every x >= X0 = 201 and x <= y <= 2x: gamma(k-1) <= gamma(k) for 3 <= k <= mu - 1/2, strictly for KL < k <= mu - 1/2 (KL = 5000). [L] Exact integers, degrees k <= KL, for all such (x, y). Modulo q^(KL+1) only finitely many series occur. [P] sigma = x s >= 6 and 0 < s <= S_TOP = 0.0182: the regime-P check of cww_certificate.py. Its analysis (Lemmas 9.3(a), 9.10, 9.11(b)) uses only that the numerator exponents y+j are at least x+2, so it applies verbatim; it is evaluated at S_TOP. [F] sigma <= 6: cells in (sigma, nu = y/x), each valid for every x >= max(X0, sigma_min/S_TOP). A degree k > KL has saddle point s < pi/sqrt(6(KL+1)) < S_TOP, and a degree k <= mu - 1/2 has s > 0 with L'(s) = mu - k >= 1/2. So [L], [P] and [F] together cover 3 <= k <= mu - 1/2. [F] and [P] use ball arithmetic (python-flint); an inequality counts only if it holds on the whole cell. The inequalities checked are those of Lemmas 9.9-9.11 and 9.17 of the paper, adapted to Gamma_{x,y} (Section 9.5): pairs (j, y+j) in the pair decomposition and the product formula for |psi|, and, beyond the Riemann-sum range, Bernoulli pairs for G = [x+y choose x]_q together with the factor [y+1]_rho / |[y+1]_{rho e^{i theta}}|, integrated over theta instead of bounded pointwise. Every hypothesis is checked with require(); nothing depends on Python's assert. Needs cww_certificate.py (Theorem C') in the same directory for the shared special functions and part [P]. Usage: python3 rq_certificate.py [--jobs J] [--parts LPF] [--lsplit i/n] (--lsplit runs the i-th of n interleaved pieces of [L], so that [L] can be spread over processes) """ import sys, time, argparse import flint from flint import arb, ctx, fmpz_poly, fmpz_series import cww_certificate as W from cww_certificate import (A, up, lo, require, CertError, PI, ZETA2, AR, NS, b_k, _q, Psibar, sinhc, sinc, H, floor3, _Gparts) ctx.prec = 96 X0 = 201 # every x >= X0 KL = 5000 # [L] degree bound S_TOP = arb('0.0182') # [P] and the lower end of the x-range of an [F] cell NX = 200 # t-cells in the Riemann sum (needs X0 > NX) THETA_M = 50 # Bernoulli blocks: scaled cells up to vartheta = THETA_M, then theta in [THETA_M/x, 1], [1, pi] # ============================================================================================= # [F] scaled variables: sigma = x s, nu = y/x in [1, 2], eps = 1/x, vartheta = x theta, u = theta sqrt(K) # ============================================================================================= def a2_bounds(sig, nu, eps): """(a2m, a2p): a2m <= kappa_2/x^3 <= a2p on the cell (Lemma 9.10 with the pairs (j, y+j)): kappa_2 = s^-2 [sum_{m=y+2}^{y+x} Psi(ms/2) - sum_{j=2}^{x} Psi(js/2)], sums compared with integrals.""" sig = A(sig); nu = A(nu); eps = A(eps) if sig.upper() < 1.2: lo_ = arb(0); hi_ = arb(0) for k in range(1, NS + 1): b = _q(b_k[k]); p = 2 * k + 1 lo_ += b * sig**(2 * k - 2) * (((1 + nu) / 2)**p - ((nu + eps) / 2)**p - ((1 + eps) / 2)**p + eps**p) hi_ += b * sig**(2 * k - 2) * (((1 + nu + eps) / 2)**p - ((nu + 2 * eps) / 2)**p - arb(2)**(-p) + (eps / 2)**p) # tail: |b_k| <= 4.4 pi^{-2k}, |bracket| <= 2 * 1.51^{2k+1} q = (up(sig) * arb('1.51') / PI)**2 require(q.upper() < 1, "series ratio") tail = arb('13.3') * (arb('1.51') / PI)**2 * q**NS / (1 - q) return lo(2 * (lo_ - arb(0, 1) * tail)), up(2 * (hi_ + arb(0, 1) * tail)) sa = lo(sig); sb = up(sig); na = lo(nu); nb = up(nu); e1 = up(eps) # lower: Psibar(s(1+nu)/2) - Psibar(s(nu+eps)/2) increases in sigma and nu and decreases in eps lo_ = 2 * (Psibar(sa * (1 + na) / 2) - Psibar(sa * (na + e1) / 2) - Psibar(sb * (1 + e1) / 2)) / sb**3 # upper: the widest interval for the numerator sum hi_ = 2 * (Psibar(sb * (1 + nb + e1) / 2) - Psibar(sa * na / 2) - Psibar(sa / 2) + Psibar(sb * e1 / 2)) / sa**3 require(lo_.lower() > 0, "a2 lower bracket must be positive") return lo(lo_), up(hi_) def S_r(r, nu, eps): """sum_{m in M} (m/2)^r <= x^{r+1} S_r for M = {2..x} u {y+2..y+x}.""" return ((1 + eps)**(r + 1) + (1 + nu + eps)**(r + 1) - nu**(r + 1)) / (arb(2)**r * (r + 1)) def deriv_consts(sig, nu, th0, eps): """k4 >= |kappa_4|/x^5; k5 >= sup_{|t| <= th0/x} |L5(s-it)|/x^6; kappa_3 <= x^4 sigma c3; sup |Re L5(s-it)| <= x^6 sigma c5. h-form (Lemma 9.3(b); needs v = (1+nu) th0/2 < pi; parity, Lemma 9.4) and, for sigma > 0.5, l-form (Lemma 9.3(a), all t); the smaller bound is used.""" sig = A(sig); nu = A(nu); eps = A(eps) v = (1 + up(nu)) * A(th0) / 2 cands = [] if v.upper() < float(PI.lower()) - 0.05: k4 = H(4, 0) * S_r(4, nu, eps) cands.append((up(k4), up(H(5, v) * S_r(5, nu, eps)), up(k4), up(H(6, v) * S_r(6, nu, eps)))) if sig.lower() > 0.5: k4l = up(AR[4] * 24 * ZETA2 / sig**5); k5l = up(AR[5] * 120 * ZETA2 / sig**6) cands.append((k4l, k5l, up(AR[3] * 6 * ZETA2 / sig**5), up(k5l / sig))) require(cands, "no derivative bound available") return (min(c[0] for c in cands), min(c[1] for c in cands), min(c[2] for c in cands), min(c[3] for c in cands)) def gauss_tail(U, j): """int_U^oo u^j e^{-u^2/4} du = 2^j Gamma((j+1)/2, U^2/4).""" U = A(U) return arb(2)**j * (U * U / 4).gamma_upper(arb(j + 1) / 2) def local_F(sig, a2m, a2p, k4, k5, c3, c5, th0, eps, ucell=0.02, Ucap=18): """Lower bounds for sqrt(K) G0 and K^{3/2} G2 of Lemma 9.9 on the Taylor zone vartheta <= th0, all x >= X0.""" sig = A(sig); eps = A(eps) cond = k4 * th0**2 / (24 * a2m) + k5 * th0**3 / (120 * a2m) require(cond.upper() < 0.25, "Taylor-zone condition") # gives Re f <= -u^2/4 on the zone uZ = th0 * (a2m * X0).sqrt() # u_Z at x = X0; it grows with x U = min(Ucap, floor3(float(uZ.lower()))) require(U >= 2 and arb(U) <= uZ, "U") ncell = int(U / ucell) + 1 du = arb(U) / ncell e15 = eps * eps.sqrt(); a15 = a2m * a2m.sqrt(); a25 = a2m**2 * a2m.sqrt() ratio = a2p / a2m G0 = arb(0); G2 = arb(0) for i in range(ncell): u = (du * i).union(du * (i + 1)) E = k4 * eps * u**4 / (24 * a2m**2) + k5 * e15 * u**5 / (120 * a25) lo_ = -ratio * u**2 / 2 - E hi_ = -u**2 / 2 + E Phi1 = c3 * u**3 / (6 * a15) + c5 * eps * u**5 / (120 * a25) m1 = min(up(sig**2 * eps * Phi1**2 / 2), arb(2)) b = m1 + u * eps * Phi1 / a2m.sqrt() + sig * u**2 * eps**3 * Phi1**2 / (4 * a2m) G0 += du * lo(lo_.exp() - hi_.exp() * b) G2 += du * lo((u**2 / 2) * (1 - u**2 * eps**3 / (12 * a2m)) * lo_.exp()) C3 = c3 / (6 * a15); C5 = c5 * eps / (120 * a25) I = lambda j: gauss_tail(arb(U), j) neg = (2 * I(0) + eps / a2m.sqrt() * (C3 * I(4) + C5 * I(6)) + sig * eps**3 / (4 * a2m) * (C3**2 * I(8) + 2 * C3 * C5 * I(10) + C5**2 * I(12))) return lo(G0 - up(neg)), lo(G2), U def Xw(w, sig, th): """X(w) = sin^2(w th/2)/sinh^2(w sig/2) = (th/sig)^2 (sinc(w th/2)/sinhc(w sig/2))^2, sigma > 0.""" return (th / sig)**2 * (sinc(w * th / 2) / sinhc(w * sig / 2))**2 def F_upper(t, sig, th, nu): """upper bound of F(t) = log(1 + X(nu+t)) - log(1 + X(t)) on the ball (Lemma 9.11: log|psi|^2 = sum_j F(j/x)).""" if A(sig).lower() >= 0.3: X1 = Xw(nu + t, sig, th); X0_ = Xw(t, sig, th) return up((1 + up(X1)).log() - (1 + lo(X0_).nonnegative_part()).log()) # small sigma: 1 + X(w) = (sig^2 S(w) + th^2 T(w)) / (sig^2 S(w)), S = sinhc^2, T = sinc^2 (sigma = 0 allowed) S1 = sinhc((nu + t) * sig / 2)**2; S0 = sinhc(t * sig / 2)**2 s1 = sinc((nu + t) * th / 2)**2; s0 = sinc(t * th / 2)**2 t1 = sig**2 * S1 + th**2 * s1 t0 = sig**2 * S0 + th**2 * s0 return up(up(t1).log() - lo(t0).log() + 2 * (up(sinhc(t * sig / 2)).log() - lo(sinhc((nu + t) * sig / 2)).log())) def F_box_upper(ta, tb, sa, sb, tha, thb, na, nb): """upper bound of F on a box: the smaller of the direct ball bound and a mean-value bound.""" ta, tb, sa, sb, tha, thb, na, nb = map(A, (ta, tb, sa, sb, tha, thb, na, nb)) t = ta.union(tb); sig = sa.union(sb); th = tha.union(thb); nu = na.union(nb) direct = F_upper(t, sig, th, nu) if ta.lower() <= 0 or sa.lower() < 0.3: return direct fc = F_upper((ta + tb) / 2, (sa + sb) / 2, (tha + thb) / 2, (na + nb) / 2) dy1, ds1, dt1 = _Gparts(nu + t, sig, th) dy0, ds0, dt0 = _Gparts(t, sig, th) cen = (fc + abs(dy1 - dy0) * (tb - ta) / 2 + abs(ds1 - ds0) * (sb - sa) / 2 + abs(dt1 - dt0) * (thb - tha) / 2 + abs(dy1) * (nb - na) / 2) return min(direct, up(cen)) def pair_cell(sa, sb, na, nb, tha, thb): """(A, B) with log|psi(theta)|^2 <= x A + B for x theta in [tha, thb], sigma in [sa, sb], nu in [na, nb] and every x >= X0 (Lemma 9.11(a) with the pairs (j, y+j); j = 1 lies in the first t-cell because x > NX).""" ts = [arb(i) / NX for i in range(NX + 1)] M = [F_box_upper(ts[c], ts[c + 1], sa, sb, tha, thb, na, nb) for c in range(NX)] require(all(m.is_finite() for m in M), "cell bounds not finite") Aval = sum(M, arb(0)) / NX TV = sum((abs(M[c] - M[c - 1]) for c in range(1, NX)), arb(0)) return up(Aval), up(TV - M[0]) def lam_bounds(sig, nu, iters=60): """Ball containing lam+, the root of h(lam) = log(1+e^lam) - log(1+e^{lam-a}) - sigma, a = sigma(1+nu). Always lam+ in [-log nu, -log nu + a] (h = sigma[(1+nu) logistic(xi) - 1], xi in [lam - a, lam]); for sigma > 0.05 refined by bisection from both sides (h increases in lam).""" sig = A(sig); nu = A(nu) l0 = -nu.log() lo_a = lo(l0); hi_b = up(l0 + sig * (1 + nu)) if sig.lower() > 0.05: h = lambda l: (1 + l.exp()).log() - (1 + (l - sig * (1 + nu)).exp()).log() - sig a0, a1 = lo_a, hi_b for _ in range(iters): m = (a0 + a1) / 2 if h(m).upper() < 0: a0 = lo(m) else: a1 = up(m) b0, b1 = lo_a, hi_b for _ in range(iters): m = (b0 + b1) / 2 if h(m).lower() > 0: b1 = up(m) else: b0 = lo(m) lo_a, hi_b = a0, b1 return lo_a.union(hi_b) def V_unit(sig, nu, xc, lam, eps_all): """lower bound for V/x, V = sum_{i= xc (eps_all = [0, 1/xc]).""" sig = A(sig); nu = A(nu) lamr = lam - sig * eps_all a1 = lamr; a2 = sig * (1 + nu) - sig * eps_all - lamr w = lambda a: (1 - (A(a) / 2).tanh()**2) / 4 # p(1-p) = 1/(4 cosh^2(a/2)), no wide cosh^2 balls w1, w2 = w(a1), w(a2) require(w1.is_finite() and w2.is_finite(), "V_unit: profile weights not finite") v = (1 + nu) * min(lo(w1), lo(w2)) # (x+y) min_i p_i(1-p_i) if sig.lower() > 1e-3: # sum >= integral - max for unimodal terms v = max(lo(v), lo(((a1 / 2).tanh() + (a2 / 2).tanh()) / (2 * sig) - arb(1) / (4 * xc))) return lo(v) def bernoulli_far(sa, sb, na, nb, th1, xc, a2p0): """At x = xc: I0 >= int |psi|, I1 >= int theta |psi| over th1/x <= theta <= pi (Lemma 9.17), and the smallest exponential rate among the blocks. a2p0 >= kappa_2(s')/x^3 for 0 <= s' <= s, so that k <= mu - 1/2 forces sigma >= 1/(2 a2p0 x^2).""" sig = A(sa).union(A(sb)); nu = A(na).union(A(nb)); x = arb(xc); eps = 1 / x eps_all = arb(0).union(eps) s = sig * eps V = V_unit(sig, nu, xc, lam_bounds(sig, nu), eps_all) * x pref = arb(4) / 3 * (4 * ((1 + up(nu)) * x / 4).sqrt() + 1) # >= 1/Pr(S = x) (Darroch, Chebyshev) sig_eff = max(lo(A(sa)), lo(1 / (2 * A(a2p0) * x**2))) require(sig_eff.lower() > 0, "sigma_eff > 0") Rp_up = up((-(sig_eff * lo(nu))).exp()) # R' = rho^{y+1} Rp_lo = lo((-(up(sig) * (up(nu) + eps))).exp()) require(Rp_up.upper() < 1 and Rp_lo.lower() > 0, "0 < R' < 1") J = up(PI / Rp_lo.sqrt() * (2 * Rp_up.sqrt() / (1 - Rp_up)).asinh()) # >= int_{-pi}^{pi} dphi/|1 - R' e^{i phi}| blocks = [] t = A(th1) while t.upper() < THETA_M: tb = t * arb('1.25') if tb.upper() > THETA_M: tb = arb(THETA_M) blocks.append((t / x, tb / x, t, None)); t = tb blocks.append((arb(THETA_M) / x, arb(1), arb(THETA_M), None)) tops = [arb('1.0'), arb('1.25'), arb('1.5'), arb('1.75'), arb('2.0'), arb('2.25'), arb('2.5'), arb('2.75'), PI] for a_, b_ in zip(tops[:-1], tops[1:]): best = None for d in (1, 2, 3): la = d * a_ / 2; hb = d * b_ / 2 if hb.upper() >= float(PI.lower()) and la.lower() <= float(PI.upper()): continue # [d a/2, d b/2] meets pi m = min(lo(la.sin()**2), lo(hb.sin()**2)) if best is None or m > best[1]: best = (d, m) require(best is not None, "top block without admissible d") blocks.append((a_, b_, None, best)) I0 = arb(0); I1 = arb(0); rate = None for (tha, thb, vta, fixed) in blocks: if fixed is not None: dmax = arb(fixed[0]); sin2 = fixed[1] else: # d = nearest integer to pi/theta dmax = PI * x / vta + arb(1) / 2; sin2 = (thb / 4).cos()**2 require(dmax.upper() <= ((1 + lo(nu)) * x - 1).lower(), "d <= x + y - 1") u = sin2 / (dmax * s / 2).cosh()**2 om = 1 - (1 - lo(u)).sqrt() # <= 1 - R for every split pair require((V - dmax / 4).lower() > 0, "V > d/4") expo = -om * (V - dmax / 4) nper = up((up(nu) + eps) * (thb - tha) * x / (2 * PI)) Q = 2 * (nper + 3) * J # >= int over the block of [y+1]_rho/|[y+1]_{rho e^{it}}| g = pref * expo.exp() * Q I0 += g I1 += g * thb r = lo(-expo / x) rate = r if rate is None or r < rate else rate return up(I0), up(I1), rate def th0_base(sig_lo, sig_hi, nu_hi): if sig_hi < 3.0: return arb(2 * 1.2) / (1 + arb(nu_hi)) # v = (1+nu) th0/2 = 1.2, as in Section 9.4 return arb(0.45 * sig_lo) # l-form bounds, valid for all t def check_F_cell(sa, sb, na, nb): """Criterion (a) or (b) of Lemma 9.9 for all sigma in [sa, sb], nu in [na, nb] and every x >= xc.""" sa = A(sa); sb = A(sb); na = A(na); nb = A(nb) sig = sa.union(sb); nu = na.union(nb) xc = max(X0, int(float((sa / S_TOP).lower()))) # for smaller x this sigma-range has k <= KL eps = arb(0).union(arb(1) / xc).nonnegative_part() a2m, a2p = a2_bounds(sig, nu, eps) th0 = None for fac in ('1', '0.9', '0.8', '0.7', '0.6', '0.5'): t0 = th0_base(float(sa.lower()), float(sb.upper()), float(nb.upper())) * arb(fac) try: k4, k5, c3, c5 = deriv_consts(sig, nu, t0, eps) G0, G2, U = local_F(sig, a2m, a2p, k4, k5, c3, c5, t0, eps) th0 = t0 break except CertError: continue require(th0 is not None, "no admissible Taylor zone") a2p0 = arb(10)**30 if sa.lower() > 0 else a2_bounds(arb(0).union(sb), nu, eps)[1] # Riemann range [th0, th1], then Bernoulli blocks x = arb(xc) th1 = arb(max(5.0, 5.0 * float(sb.upper()))) th = A(th0) I0 = arb(0); I1 = arb(0); worstA = None while th.upper() < th1.upper(): w = arb('0.05') if th.upper() < 2 else (arb('0.1') if th.upper() < 6 else arb('0.25')) thb = th + w if thb.upper() > th1.upper(): thb = th1 Aa, Bb = pair_cell(sa, sb, na, nb, th, thb) worstA = Aa if worstA is None or Aa > worstA else worstA e = ((x * Aa + Bb) / 2).exp() I0 += (thb - th) / x * e I1 += (thb - th) * thb / x**2 * e th = thb J0, J1, rate = bernoulli_far(sa, sb, na, nb, th1, xc, a2p0) I0 += J0; I1 += J1 U0 = a2m.sqrt() * x**(arb(3) / 2) * I0 # sqrt(K) I0 at x = xc KU1 = a2m**(arb(3) / 2) * x**(arb(9) / 2) * I1 # K^{3/2} I1 # Monotonicity in x >= xc. Riemann part: sqrt(K) I0 and sqrt(K) I1/(1-rho) scale like x^{1/2} e^{xA/2} # (need A < -1/x), K^{3/2} I1 like x^{5/2} e^{xA/2} (need A < -5/x). Bernoulli part: at most # (x/xc)^7 e^{-rate (x - xc)} (x^{9/2} from K^{3/2}, x^{1/2} prefactor, x from the pole count, x from J or # 1/(1-rho)); decreasing when rate * xc > 7. monoB = bool((rate * xc).lower() > 7) critA = bool((G0 - U0) > 0) and bool((G2 - KU1) > 0) and bool(worstA < -arb(5) / xc) critB = False; marginB = None if sa.lower() > 0: onem = sa / x * (1 - sb / (2 * x)) # <= 1 - rho U1r = a2m.sqrt() * x**(arb(3) / 2) * I1 / onem marginB = float((G0 - U0 - U1r).lower()) critB = bool((G0 - U0 - U1r) > 0) and bool(worstA < -arb(1) / xc) ok = monoB and (critA or critB) return ok, dict(sig=(float(sa.lower()), float(sb.upper())), nu=(float(na.lower()), float(nb.upper())), xc=xc, G0=float(G0.lower()), G2=float(G2.lower()), U0=float(U0.upper()), KU1=float(KU1.upper()), critA=critA, critB=critB, marginB=marginB, A=float(worstA.upper()), rate_xc=float((rate * xc).lower())) def sigma_cells(): out = [(0.0, 1e-3)] v = 1e-3 while v < 0.3: w = min(v * 1.25, 0.3); out.append((v, w)); v = w while v < 6.0 - 1e-12: w = min(v * 1.02, 6.0); out.append((v, w)); v = w return out def _cell_worker(args): sa, sb, na, nb = args t = time.time() try: ok, r = check_F_cell(arb(sa), arb(sb), arb(na), arb(nb)) except CertError as e: ok, r = False, {'failed hypothesis': str(e)[:200]} r['time'] = time.time() - t return args, ok, r def part_F(jobs): """sigma-cells x nu-cells of width 0.1; a failing cell is split into nu-widths 0.02 and then 0.005.""" from multiprocessing import Pool nu_edges = [1 + 0.1 * i for i in range(10)] + [2.0] todo = [(a, b, nu_edges[i], nu_edges[i + 1]) for (a, b) in sigma_cells() for i in range(10)] done = []; failed = [] pool = Pool(jobs) if jobs > 1 else None for split in (5, 4, None): res = pool.map(_cell_worker, todo, chunksize=1) if pool else [_cell_worker(c) for c in todo] todo = [] for (args, ok, r) in res: if ok: done.append(r) elif split is not None: sa, sb, na, nb = args h = (nb - na) / split edges = [na + h * i for i in range(split)] + [nb] # exact endpoints: no gaps todo += [(sa, sb, edges[i], edges[i + 1]) for i in range(split)] else: failed.append((args, r)) if not todo: break if pool: pool.close() return done, failed # ============================================================================================= # [L] exact low degrees # ============================================================================================= def part_L_single(KL=KL): """y > YS = floor((KL-5)/2): at most one numerator factor below degree KL+1, so with z = coefficients of prod_{j=2}^{x} 1/(1-q^j): gamma(k) - gamma(k-1) = dz(k) - [z(k-y-2) - z(k-y-x-1)] >= dz(k) - zmax(k-YS-3). Check dz(k) >= zmax(k-YS-3) for 3 <= k <= KL and every x >= X0 (x = KL stands for all x >= KL).""" YS = (KL - 5) // 2 z = [0] * (KL + 1); z[0] = 1 for j in range(2, X0): for i in range(j, KL + 1): z[i] += z[i - j] bad = [] for x in range(X0, KL + 1): for i in range(x, KL + 1): z[i] += z[i - x] zmax = []; m = 0 for v in z: m = max(m, v); zmax.append(m) for k in range(3, KL + 1): u = k - YS - 3 if z[k] - z[k - 1] < (zmax[u] if u >= 0 else 0): bad.append((x, k)); break return bad def part_L_exact(KL=KL, report=250000, part=0, parts=1, ya=None, yb=None): """y <= YS: Gamma_{x,y} modulo q^(KL+1), incremental in x for each y; gamma(k-1) <= gamma(k) for 3 <= k <= KL. part/parts: this call handles the y with y % parts == part; ya <= y < yb restricts y further (each y is independent, so the work can be split and resumed by y-ranges).""" YS = (KL - 5) // 2 flint.ctx.cap = KL + 1 inv = {} bad = []; count = 0; t0 = time.time() lo_y = X0 if ya is None else max(X0, ya) hi_y = YS + 1 if yb is None else min(YS + 1, yb) for y in range(lo_y, hi_y): if y % parts != part: continue xa = max(X0, (y + 1) // 2) Num = fmpz_poly([1]); Den = fmpz_poly([1]) for j in range(2, xa + 1): if y + j <= KL: Num = (Num - Num.left_shift(y + j)).truncate(KL + 1) Den = (Den - Den.left_shift(j)).truncate(KL + 1) G = Num.mul_low(fmpz_poly((1 / fmpz_series(Den.coeffs(), prec=KL + 1)).coeffs()), KL + 1) for x in range(xa, y + 1): if x > xa: # Gamma_{x,y} = Gamma_{x-1,y} (1-q^{y+x})/(1-q^x) if y + x <= KL: G = (G - G.left_shift(y + x)).truncate(KL + 1) if x not in inv: inv[x] = fmpz_poly([1 if i % x == 0 else 0 for i in range(KL + 1)]) G = G.mul_low(inv[x], KL + 1) top = min(KL, ((x - 1) * y) // 2 - 1) d = (G - G.left_shift(1)).truncate(top + 1).coeffs() if any(d[k] < 0 for k in range(3, len(d))): bad.append((x, y)) count += 1 if count % report == 0: print(" [L] ... %d classes, y = %d, %.0fs" % (count, y, time.time() - t0), flush=True) return bad, count # ============================================================================================= def main(): ap = argparse.ArgumentParser() ap.add_argument('--jobs', type=int, default=1) ap.add_argument('--parts', default='LPF') ap.add_argument('--lsplit', default='0/1', help="i/n: run only the i-th of n interleaved pieces of [L]") ap.add_argument('--yrange', default=None, help="a:b: in [L], only the y with a <= y < b; 'single': only y > (KL-5)/2") args = ap.parse_args() li, ln = (int(v) for v in args.lsplit.split('/')) only_single = args.yrange == 'single' ya, yb = (int(v) for v in args.yrange.split(':')) if (args.yrange and not only_single) else (None, None) t0 = time.time() allok = True print("Certificate for Theorem 9.16 (ii) (rational q-Catalan monotonicity, x >= %d, x <= y <= 2x); arb prec = %d" % (X0, ctx.prec)) sK = PI / (6 * arb(KL + 1)).sqrt() require(sK < S_TOP, "s_K < S_TOP") print("k > KL = %d forces s < pi/sqrt(6(KL+1)) = %.6f < S_TOP = %s" % (KL, float(sK.upper()), S_TOP.str(4))) if 'P' in args.parts or 'F' in args.parts: t = time.time() W.part_A() # the constants A_r used below (Section 9.4) print("[A] A_r bounds verified %.0fs" % (time.time() - t)) if 'L' in args.parts: t = time.time() bad1 = part_L_single() if (li == 0 and ya is None) else [] if li == 0 and ya is None: print("[L] y > (KL-5)/2, all x >= %d: %s %.0fs" % (X0, bad1[:5] or "verified", time.time() - t), flush=True) t = time.time() bad2, cnt = ([], 0) if only_single else part_L_exact(part=li, parts=ln, ya=ya, yb=yb) if not only_single: print("[L] y <= (KL-5)/2 (piece %d/%d%s): %d series modulo q^(KL+1): %s %.0fs" % (li, ln, "" if ya is None else ", %d <= y < %d" % (ya, yb), cnt, bad2[:5] or "verified", time.time() - t), flush=True) allok &= not bad1 and not bad2 if 'P' in args.parts: t = time.time() try: ok, r = W.check_P(S_TOP) except CertError as e: ok, r = False, {'failed hypothesis': str(e)} allok &= ok print("[P] sigma >= 6, 0 < s <= %s: %s %s %.0fs" % (S_TOP.str(4), "verified" if ok else "FAILED", r, time.time() - t), flush=True) if 'F' in args.parts: t = time.time() done, failed = part_F(args.jobs) for (a, r) in failed: print(" FAILED cell sigma [%.6g, %.6g] nu [%.4g, %.4g]: %s" % (a + (r,))) allok &= not failed if done: mA = min((r['G0'] - r['U0'] for r in done if r['critA']), default=float('nan')) mB = min((r['marginB'] for r in done if r['critB'] and not r['critA']), default=float('nan')) print("[F] 0 <= sigma <= 6, 1 <= nu <= 2: %d cells verified, %d failed; min(G0 - sqrt(K) I0) = %.4f (criterion a)," " min margin %.4f (criterion b only); max A = %.4f; min rate*xc = %.1f %.0fs" % (len(done), len(failed), mA, mB, max(r['A'] for r in done), min(r['rate_xc'] for r in done), 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())