#!/usr/bin/env python3
"""Certificate for Theorem C' of the paper "Analytic q-inequalities and a proof of Bergeron's conjecture" (Section 9.4).

    C_n(q) = [2n choose n]_q / [n+1]_q = prod_{j=2}^{n} (1 - q^{n+j}) / (1 - q^j) = sum_k c_k q^k,
    mu = n(n-1)/2.   Theorem C':  c_{k-1} <= c_k for 2 <= k <= mu, for every n >= 16.

This script proves the part n >= 1001 (corner_check.py covers n <= 1000):
  [L] Exact integers.  For 16 <= n <= K0 = 10^4 and 2 <= k <= min(K0, mu): c_{k-1} <= c_k, and
      c_{k-1} < c_k if also k >= n+2.  Modulo q^{K0+1}, C_n is the same polynomial for all n >= K0,
      so this covers every n >= 16.
  [F] Degrees whose saddle point s satisfies sigma = n s <= 6, for every n >= 1001: 179 sigma-cells.
  [P] Degrees with sigma >= 6 and 0 < s <= S_TOP = 0.01283, for every n: one s-cell.
A degree k > K0 has s < pi/sqrt(6 (K0+1)) < S_TOP (the tilted mean is below pi^2/(6 s^2)), so for
n >= 1001 every degree K0 < k <= mu falls under [F] or [P], and there c_k > c_{k-1} (Lemma 9.2).

[F] and [P] use ball arithmetic from python-flint (pip install python-flint): an inequality is accepted
only if it holds at every point of the balls involved.  The inequalities checked, and why they suffice,
are Lemmas 9.9-9.12 and the proof of Theorem C' in the paper.

Every hypothesis is checked with require(), which raises CertError; nothing depends on Python's assert.

Usage:  python3 cww_certificate.py [--jobs J] [--parts ALFP]      (F and P always run A first)
"""
import sys, time, argparse
from fractions import Fraction
from math import comb, factorial
from flint import arb, ctx

ctx.prec = 96
PI = arb.pi()
ZETA2 = PI**2 / 6
N0 = 1001                    # [F] holds for all n >= N0
K0 = 10000                   # [L] degree bound
S_TOP = arb('0.01283')       # [P] covers 0 < s <= S_TOP


class CertError(Exception):
    pass


def require(cond, msg):
    """a hypothesis of the proof; failure makes the cell (or the run) fail."""
    if not cond:
        raise CertError(msg)


def A(x):
    return x if isinstance(x, arb) else arb(x)


def up(x):                   # exact upper endpoint, as an arb
    return arb(A(x).upper())


def lo(x):                   # exact lower endpoint, as an arb
    return arb(A(x).lower())


# ---------------------------------------------------------------------------------------------
# Special functions.  h(y) = log(sinh y / y) = sum_k a_k y^{2k},  a_k = (-1)^{k+1} zeta(2k)/(k pi^{2k}).
# Psi(y) = 1 - y^2/sinh^2 y = y^2 h''(y) (increasing, Psi(0) = 0),  Psibar(y) = int_0^y Psi.
# ---------------------------------------------------------------------------------------------
def _bernoulli(m):
    B = [Fraction(1)]
    for k in range(1, m + 1):
        B.append(-sum(comb(k + 1, j) * B[j] for j in range(k)) / Fraction(k + 1))
    return B


NS = 120
_B = _bernoulli(2 * NS + 2)
a_k = [Fraction(0)] + [_B[2 * k] * 2**(2 * k) / (2 * k * factorial(2 * k)) for k in range(1, NS + 1)]
b_k = [None] + [a_k[k] * (2 * k) * (2 * k - 1) / (2 * k + 1) for k in range(1, NS + 1)]   # Psibar = sum b_k y^{2k+1}


def _q(x):
    return arb(x.numerator) / arb(x.denominator)


def Psibar(y):
    """int_0^y (1 - t^2/sinh^2 t) dt.  Series for y < 1 (|a_k| 2k(2k-1) <= 8.8 k pi^{-2k} bounds the tail)."""
    y = A(y)
    if y.upper() < 1:
        s = arb(0)
        for k in range(1, NS + 1):
            s += _q(b_k[k]) * y**(2 * k + 1)
        q = (up(y) / PI)**2
        tail = arb('8.8') * (NS + 1) * q**(NS + 1) / (1 - q)**2 * up(y)
        return s + arb(0, 1) * tail
    e = (-2 * y).exp()
    return y + y * y / y.tanh() - y * y - 2 * y * (1 - e).log() + e.polylog(2) - PI**2 / 6


def sinhc(y):                # sinh(y)/y
    y = A(y)
    if y.upper() < 0.5 and y.lower() > -0.5:
        s = arb(0); t = arb(1)
        for k in range(40):
            s += t
            t = t * y * y / ((2 * k + 2) * (2 * k + 3))
        return s + arb(0, 1) * abs(t) * 2
    return y.sinh() / y


def sinc(y):                 # sin(y)/y
    return A(y).sinc()


def tanhc(y):                # tanh(y)/y
    y = A(y)
    if y.upper() < 0.5:
        return sinhc(y) / y.cosh()
    return y.tanh() / y


def H(r, v):
    """(r-1)! * 2 sum_{l>=1} (pi l - v)^{-r}: bound for |h^{(r)}(y)| when |Im y| <= v < pi (Lemma 9.3(b))."""
    return arb(factorial(r - 1)) * 2 * PI**(-r) * arb(r).zeta(1 - A(v) / PI)


def sig_r(r, eps):
    """sum_{m in M} (m/2)^r <= n^{r+1} (2 + eps)^{r+1} / (2^r (r+1)),  eps = 1/n."""
    return (2 + A(eps))**(r + 1) / (arb(2)**r * (r + 1))


# ---------------------------------------------------------------------------------------------
# Constants A_r = sup_{x>0} x^{r+1} Li_{-r}(e^{-x}) / r!  (Lemma 9.3(a) then gives
# |L^{(r)}(s - it)| <= A_r r! zeta(2) / s^{r+1} for all real t).
# x^{r+1} Li_{-r}(e^{-x}) = (x/(1-e^{-x}))^{r+1} * e^{-x} E_r(e^{-x}), E_r = Eulerian polynomial:
# the first factor increases and the second decreases, so cell endpoints give rigorous bounds.
# ---------------------------------------------------------------------------------------------
EULER = {3: [1, 4, 1], 4: [1, 11, 11, 1], 5: [1, 26, 66, 26, 1], 6: [1, 57, 302, 302, 57, 1]}
AR = {3: arb('1.011'), 4: arb('1.039'), 5: arb('1.08'), 6: arb('1.13')}   # values used below


def A_r_bound(r, x1=120, cells=60000):
    def f1(x):
        x = A(x)
        if x.upper() < 1e-6:
            return arb(1) + arb(0, 1e-6)
        return x / (-(-x).expm1())

    def f2(x):
        z = (-A(x)).exp()
        return z * sum(c * z**i for i, c in enumerate(EULER[r]))

    best = arb(0)
    for i in range(cells):
        v = f1(arb(x1) * (i + 1) / cells)**(r + 1) * f2(arb(x1) * i / cells)
        require(v.is_finite(), ("A_r cell not finite", r, i))
        if v.upper() > best.upper():
            best = up(v)
    tail = arb(x1)**(r + 1) * (-arb(x1)).exp() * 2**(r + 1) * sum(EULER[r])   # x >= x1
    require(tail.upper() < best.upper(), ("A_r tail", r))
    return best / factorial(r)


def part_A():
    out = {}
    for r in (3, 4, 5, 6):
        b = A_r_bound(r)
        require(b.upper() < AR[r].lower(), ("A_r", r, b))
        out[r] = float(b.upper())
    return out


# =============================================================================================
# [F]  0 <= sigma <= 6, n >= N0.  Scaled variables: sigma = n s, eps = 1/n in [0, 1/N0],
# vartheta = n theta, u = theta sqrt(K) with K = n^3 a2m <= kappa_2.
# =============================================================================================
EPS = arb(0).union(arb(1) / N0).nonnegative_part()


def a2_bounds(sig, eps=EPS):
    """(a2m, a2p) with a2m <= kappa_2/n^3 <= a2p for sigma in sig and 1/n in eps (Lemma 9.10)."""
    sig = A(sig); 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 - 2 * ((1 + eps) / 2)**p + eps**p)
            hi_ += b * sig**(2 * k - 2) * ((1 + eps / 2)**p - ((1 + 2 * eps) / 2)**p - arb(2)**(-p) + (eps / 2)**p)
        q = (up(sig) * arb('1.01') / PI)**2
        tail = 2 * arb('8.8') * (NS + 1) * q**NS / (1 - q)**2
        return lo(2 * (lo_ - arb(0, 1) * tail)), up(2 * (hi_ + arb(0, 1) * tail))
    # closed form at the cell endpoints (Psibar increases; both brackets are positive)
    sa = lo(sig); sb = up(sig); e1 = up(eps)
    lo_ = 2 * (Psibar(sa) - 2 * Psibar(sb * (1 + e1) / 2)) / sb**3
    hi_ = 2 * (Psibar(sb * (1 + e1 / 2)) - 2 * Psibar(sa / 2) + Psibar(sb * e1 / 2)) / sa**3
    require(lo_.lower() > 0 and hi_.lower() > 0, "a2 brackets must be positive")
    return lo(lo_), up(hi_)


def deriv_consts(sig, th0, eps=EPS):
    """k4 >= |kappa_4|/n^5;  k5 >= sup_{|t| <= th0/n} |L^(5)(s - it)|/n^6;
    kappa_3 <= n^4 sigma c3;  sup_{|t| <= th0/n} |Re L^(5)(s - it)| <= n^6 sigma c5   (Lemmas 9.3, 9.4)."""
    sig = A(sig); eps = A(eps)
    h4 = H(4, 0) * sig_r(4, eps)
    k4 = h4
    k5 = H(5, th0) * sig_r(5, eps)
    c3 = h4
    c5 = H(6, th0) * sig_r(6, eps)
    if sig.lower() > 0.5:                     # l-form bounds, valid for s > 0 and all t
        k4 = min(up(k4), up(AR[4] * 24 * ZETA2 / sig**5))
        k5 = min(up(k5), up(AR[5] * 120 * ZETA2 / sig**6))
        c3 = min(up(c3), up(AR[3] * 6 * ZETA2 / sig**5))
        c5 = min(up(c5), up(k5 / sig))
    return up(k4), up(k5), up(c3), up(c5)


def gauss_tail(U, j):
    """I_j >= int_U^oo u^j e^{-u^2/4} du   (integration by parts; needs U^2 > 2(j-1))."""
    U = A(U)
    if j == 0:
        return 2 * (-U**2 / 4).exp() / U
    require((U**2).lower() > 2 * (j - 1), ("gauss_tail", j))
    return 2 * U**(j - 1) * (-U**2 / 4).exp() / (1 - 2 * arb(j - 1) / U**2)


def floor3(x):
    return int(x * 1000) / 1000.0


def local_F(sig, a2m, a2p, k4, k5, c3, c5, th0, ucell=0.02, eps=EPS, Ucap=18):
    """Lower bounds G0 for int_0^{uZ} e^{Re f}(1 - b) du and G2 for int_0^{uZ} (1 - cos theta) K e^{Re f} du,
    uZ = th0 sqrt(a2m n) >= U (Lemma 9.9, Taylor zone vartheta <= th0)."""
    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", cond))    # gives Re f <= -u^2/4 on the zone
    uZ = th0 * (a2m * N0).sqrt()                                        # u_Z at n = N0; it grows with n
    U = min(Ucap, floor3(float(uZ.lower())))
    require(U >= 5 and arb(U) <= uZ, ("U", 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())
    # [U, uZ]: drop the positive part, bound the negative part by int_U^oo e^{-u^2/4} b(u) du
    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 F_pair(x, sig, th):
    """Ball bounds (upper, lower) of F(x) = log(1 + X(1+x)) - log(1 + X(x)),
    X(y) = sin^2(y vartheta/2)/sinh^2(y sigma/2), written as in Lemma 9.11 so that sigma = 0 is allowed."""
    S1 = sinhc((1 + x) * sig / 2)**2; S0 = sinhc(x * sig / 2)**2
    s1 = sinc((1 + x) * th / 2)**2;  s0 = sinc(x * th / 2)**2
    t1 = sig**2 * S1 + th**2 * s1
    t0 = sig**2 * S0 + th**2 * s0
    upper = (up(t1).log() - lo(t0).log()
             + 2 * (up(sinhc(x * sig / 2)).log() - lo(sinhc((1 + x) * sig / 2)).log()))
    if t1.lower() > 0:
        lower = (lo(t1).log() - up(t0).log()
                 + 2 * (lo(sinhc(x * sig / 2)).log() - up(sinhc((1 + x) * sig / 2)).log()))
        lower = lo(lower)
    else:
        lower = None
    return up(upper), lower


def _Gparts(y, sig, th):
    """partial derivatives of G(y) = log(1 + sin^2(y vartheta/2)/sinh^2(y sigma/2)) in y, sigma, vartheta."""
    D = (y * sig / 2).sinh()**2 + (y * th / 2).sin()**2
    cth = 1 / (y * sig / 2).tanh()
    dy = ((sig / 2) * (y * sig).sinh() + (th / 2) * (y * th).sin()) / D - sig * cth
    ds = y * ((y * sig).sinh() / (2 * D) - cth)
    dt = (y / 2) * (y * th).sin() / D
    return dy, ds, dt


def F_box_upper(xa, xb, sa, sb, ta, tb):
    """upper bound of F on the box [xa,xb] x [sa,sb] x [ta,tb]: the smaller of the direct ball bound and
    the mean-value (centred) bound F(centre) + sum_i sup|dF/dz_i| * halfwidth_i."""
    xa, xb, sa, sb, ta, tb = map(A, (xa, xb, sa, sb, ta, tb))
    x = xa.union(xb); sig = sa.union(sb); th = ta.union(tb)
    direct, _ = F_pair(x, sig, th)
    if xa.lower() <= 0 or sa.lower() < 0.3:
        return direct
    fc, _ = F_pair((xa + xb) / 2, (sa + sb) / 2, (ta + tb) / 2)
    dy1, ds1, dt1 = _Gparts(1 + x, sig, th)
    dy0, ds0, dt0 = _Gparts(x, sig, th)
    cen = fc + abs(dy1 - dy0) * (xb - xa) / 2 + abs(ds1 - ds0) * (sb - sa) / 2 + abs(dt1 - dt0) * (tb - ta) / 2
    return min(direct, up(cen))


def pair_cell(sa, sb, tha, thb, nx=200):
    """(A, B) with log|psi(theta)|^2 <= n A + B for vartheta = n theta in [tha, thb], sigma in [sa, sb],
    and every n >= nx (Lemma 9.11(a): Riemann sum with Abel summation)."""
    xs = [arb(i) / nx for i in range(nx + 1)]
    M = [F_box_upper(xs[c], xs[c + 1], sa, sb, tha, thb) for c in range(nx)]
    _, L0 = F_pair(xs[0].union(xs[1]), A(sa).union(A(sb)), A(tha).union(A(thb)))
    require(L0 is not None and L0.is_finite(), "lower bound on the first x-cell")
    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))
    B = TV if L0 >= 0 else TV - L0
    return up(Aval), up(B)


def bernoulli_tail(sig, th1, n):
    """(beta, rate): |psi(theta)| <= e^beta for th1/n <= theta <= pi (Lemma 9.12), with
    beta = log(3n+1) + log(4/3 (4 sqrt(n/2) + 1)) - rate*n + (3/8) omega, and rate = omega*Cp."""
    sig = A(sig); n = A(n); eps = 1 / n
    s = sig * eps
    omega = 1 / (4 * ((PI * sig / th1 + s / 2) / 2).cosh()**2)          # <= 1 - R
    z = sig * (1 - eps / 2) / 2
    Cp = (1 - eps / 2) * tanhc(z) / 2 - PI / (4 * th1)                   # sum pi_i >= n Cp - 3/8
    V_lo = n * (1 - eps / 2) * tanhc(z) / 2 - arb(1) / 4
    require(V_lo.lower() > 12, "V_- > 12")
    require((PI * n / th1 + arb(1) / 2).upper() <= n.lower(), "d <= n")
    require((n * Cp).lower() > 3 / 8, "n Cp > 3/8")
    beta = (3 * n + 1).log() + (arb(4) / 3 * (4 * (n / 2).sqrt() + 1)).log() - omega * (n * Cp - arb(3) / 8)
    return beta, omega * Cp


def th0_for(sig_hi):
    return arb('1.2') if sig_hi < 3.0 else arb(0.4 * sig_hi)


def far_F(sa, sb, th0, nx=200):
    """At n = N0: I0 >= int_{th0/n}^pi |psi|, I1 >= int_{th0/n}^pi theta |psi|, for sigma in [sa, sb];
    also the largest pair-cell rate A and the Bernoulli rate (the n-monotonicity conditions)."""
    n = arb(N0)
    th1 = arb(max(5.0, 5.0 * float(A(sb).upper())))
    th = A(th0)
    I0 = arb(0); I1 = arb(0); worstA = None
    while th.upper() < th1.upper():
        w = arb('0.02') if th.upper() < 2 else (arb('0.05') if th.upper() < 6 else arb('0.25'))
        thb = th + w
        if thb.upper() > th1.upper():
            thb = th1
        Aa, Bb = pair_cell(sa, sb, th, thb, nx)
        worstA = Aa if worstA is None or Aa > worstA else worstA
        e = ((n * Aa + Bb) / 2).exp()
        I0 += (thb - th) / n * e
        I1 += (thb - th) * thb / n**2 * e
        th = thb
    beta, rate = bernoulli_tail(A(sa).union(A(sb)), th1, n)
    eb = beta.exp()
    I0 += PI * eb
    I1 += PI**2 / 2 * eb
    return I0, I1, worstA, rate


def check_F_cell(sa, sb):
    """Criterion of Lemma 9.9 for all sigma in [sa, sb] and all n >= N0."""
    sa = A(sa); sb = A(sb)
    sig = sa.union(sb)
    th0 = th0_for(float(sb.upper()))
    a2m, a2p = a2_bounds(sig)
    k4, k5, c3, c5 = deriv_consts(sig, th0)
    G0, G2, U = local_F(sig, a2m, a2p, k4, k5, c3, c5, th0)
    I0, I1, worstA, rate = far_F(sa, sb, th0)
    n = arb(N0)
    # every n >= N0: the far bounds times the scale factors decrease in n when A < -5/n and rate > 6/n
    mono = (worstA < -arb(5) / N0) and (rate > arb(6) / N0)
    U0 = a2m.sqrt() * n**(arb(3) / 2) * I0                     # sqrt(K) I0
    KU1 = a2m**(arb(3) / 2) * n**(arb(9) / 2) * I1             # K^{3/2} I1
    okA = (G0 - U0) > 0
    okB = (G2 - KU1) > 0
    okC = False
    margin_C = None
    if not okB and sa.lower() > 0:
        onem = sa / n * (1 - sb / (2 * n))                     # <= 1 - rho
        U1r = a2m.sqrt() * n**(arb(3) / 2) * I1 / onem
        okC = (G0 - U0 - U1r) > 0
        margin_C = float((G0 - U0 - U1r).lower())
    ok = bool(mono and okA and (okB or okC))
    return ok, dict(sig=(float(sa.lower()), float(sb.upper())), G0=float(G0.lower()), G2=float(G2.lower()),
                    U0=float(U0.upper()), KU1=float(KU1.upper()), okB=bool(okB), okC=bool(okC), margin_C=margin_C,
                    A=float(worstA.upper()), rate=float(rate.lower()), U=U, a2=(float(a2m), float(a2p)))


def F_cells():
    out = [(0.0, 1e-3)]
    x = 1e-3
    while x < 0.3:
        y = min(x * 1.25, 0.3); out.append((x, y)); x = y
    while x < 6.0 - 1e-12:
        y = min(x * 1.02, 6.0); out.append((x, y)); x = y
    return out


def _F_worker(ab):
    a, b = ab
    t = time.time()
    try:
        ok, r = check_F_cell(arb(a), arb(b))
    except (CertError, AssertionError) as e:
        ok, r = False, {'failed hypothesis': str(e)[:200]}
    r['time'] = time.time() - t
    return a, b, ok, r


# =============================================================================================
# [P]  sigma >= 6, 0 < s <= S_TOP, all n.  Scale K = k2m/s^3 <= kappa_2, u = theta sqrt(K).
# =============================================================================================
LAM3 = AR[3] * 6 * ZETA2; LAM4 = AR[4] * 24 * ZETA2; LAM5 = AR[5] * 120 * ZETA2
E6 = arb(-6).exp()
G6 = 50 * E6 / (1 - E6)**2          # >= int_6^oo x^2 / (4 sinh^2(x/2)) dx
GAM = arb('0.5')                    # Taylor zone theta <= GAM s
K2P = PI**2 / 3                     # s^3 kappa_2 <= pi^2/3
NU = 2 / (arb(6).exp() - 1)         # (1/2) sum_{m >= n+2} X_m <= NU / s


def local_P(sb, ucell=0.02, Ucap=18):
    """G0 as in local_F, valid for every s in (0, sb] (each ingredient is monotone in s)."""
    s = up(sb)
    k2m = K2P - 2 * s - 2 * G6
    cond = LAM4 * GAM**2 / (24 * k2m) + LAM5 * GAM**3 / (120 * k2m)
    require(cond.upper() < 0.25, ("Taylor-zone condition", cond))
    uZ = GAM * (k2m / s).sqrt()                                         # u_Z at s = sb; larger for smaller s
    U = min(Ucap, floor3(float(uZ.lower())))
    require(U >= 5 and arb(U) <= uZ, ("U", U))
    ncell = int(U / ucell) + 1
    du = arb(U) / ncell
    s15 = s * s.sqrt(); k15 = k2m * k2m.sqrt(); k25 = k2m**2 * k2m.sqrt()
    G0 = arb(0)
    for i in range(ncell):
        u = (du * i).union(du * (i + 1))
        E = LAM4 * s * u**4 / (24 * k2m**2) + LAM5 * s15 * u**5 / (120 * k25)
        lo_ = -(K2P / k2m) * u**2 / 2 - E
        hi_ = -u**2 / 2 + E
        Phi2 = LAM3 * u**3 / (6 * k15) + LAM5 * s * u**5 / (120 * k25)
        m1 = min(up(s * Phi2**2 / 2), arb(2))
        b = m1 + u * s * Phi2 / k2m.sqrt() + u**2 * s**3 * Phi2**2 / (4 * k2m)
        G0 += du * lo(lo_.exp() - hi_.exp() * b)
    C3 = LAM3 / (6 * k15); C5 = LAM5 * s / (120 * k25)
    I = lambda j: gauss_tail(arb(U), j)
    neg = (2 * I(0) + s / k2m.sqrt() * (C3 * I(4) + C5 * I(6))
           + s**3 / (4 * k2m) * (C3**2 * I(8) + 2 * C3 * C5 * I(10) + C5**2 * I(12)))
    return lo(G0 - up(neg)), U


def lower_riemann(f, y0, y1, ny, geometric=False):
    """lower bound of int_{y0}^{y1} f for f >= 0 decreasing: sum of (y_{i+1} - y_i) f(y_{i+1})."""
    y0 = A(y0); y1 = A(y1)
    if y1.lower() <= y0.upper():
        return arb(0)
    if geometric:
        r = (y1 / y0)**(arb(1) / ny)
        ys = [y0 * r**i for i in range(ny)] + [y1]
    else:
        ys = [y0 + (y1 - y0) * i / ny for i in range(ny)] + [y1]
    tot = arb(0)
    for i in range(ny):
        tot += (ys[i + 1] - ys[i]) * lo(f(ys[i + 1]))
    return lo(tot)


def far_P(sb, cx=arb(20), ratio=arb('1.1'), ny=240):
    """log|psi(theta)| <= -kappa/s for GAM s <= theta <= pi, every s in (0, sb], every n with n s >= 6
    (Lemma 9.11(b)).  Zone 1: theta = c s with c in [ca, cb] (cells up to cx); zone 2: theta >= cx s."""
    sb = up(sb)
    zone1 = []
    c = GAM
    while c.upper() < cx.upper():
        cb = c * ratio
        if cb.upper() > cx.upper():
            cb = cx
        Y = min(lo(PI / cb), arb(6))
        J = lower_riemann(lambda y: (1 + (c * y / 2).sin()**2 / (y / 2).sinh()**2).log(), up(2 * sb), Y, ny)
        zone1.append((c, cb, lo(J / 2 - NU)))
        c = cb
    SIN2 = (PI / 8).sin()**2
    Bg = lower_riemann(lambda y: (1 + SIN2 / (y / 2).sinh()**2).log(), up(PI / cx + 2 * sb),
                       lo(6 - 2 * PI / cx - sb), 400, geometric=True)
    kap2 = lo(Bg / 4 - NU)
    return zone1, kap2


def check_P(sb=S_TOP):
    G0, U = local_P(sb)
    zone1, kap2 = far_P(sb)
    ks = [lo(k) for (_, _, k) in zone1] + [lo(kap2)]
    kmin = ks[0]
    for k in ks[1:]:
        if k < kmin:                       # exact endpoints: exact comparison
            kmin = k
    s = up(sb)
    # the far terms below increase with s on (0, sb] when kmin/sb > 3, so s = sb is the worst case
    require((kmin / s).lower() > 3, "kappa_min / s > 3")
    Ek = lambda k: (-k / s).exp()
    S0 = sum(((cb - c) * s * Ek(k) for (c, cb, k) in zone1), arb(0)) + PI * Ek(kap2)
    S1 = sum(((cb - c) * s * cb * s * Ek(k) for (c, cb, k) in zone1), arb(0)) + PI**2 / 2 * Ek(kap2)
    sqK = K2P.sqrt() / (s * s.sqrt())                 # >= sqrt(kappa_2)
    onem = s * (1 - s / 2)                            # <= 1 - rho
    val = G0 - sqK * S0 - sqK * S1 / onem
    return bool(val > 0), dict(G0=float(G0.lower()), U=U, kmin=float(kmin.lower()),
                               far0=float((sqK * S0).upper()), far1=float((sqK * S1 / onem).upper()),
                               margin=float(val.lower()))


# =============================================================================================
# [L]  exact low degrees
# =============================================================================================
def part_L(K=K0, nmin=16):
    """C_n = C_{n-1} (1 + q^n)(1 - q^{2n-1})/(1 - q^{n+1}) modulo q^{K+1}, for n = 2..K."""
    C = [0] * (K + 1); C[0] = 1                        # C_1 = 1
    bad = []; weak = []
    for n in range(2, K + 1):
        if n <= K:
            for i in range(K, n - 1, -1):
                C[i] += C[i - n]
        m = 2 * n - 1
        if m <= K:
            for i in range(K, m - 1, -1):
                C[i] -= C[i - m]
        m = n + 1
        if m <= K:
            for i in range(m, K + 1):
                C[i] += C[i - m]
        if n >= nmin:
            top = min(K, n * (n - 1) // 2)
            for k in range(2, top + 1):
                d = C[k] - C[k - 1]
                if d < 0:
                    bad.append((n, k)); break
                if d == 0 and k >= n + 2:
                    weak.append((n, k))
    # sanity: for n = K the truncation is prod_{j=2}^K 1/(1-q^j), i.e. partitions into parts >= 2
    P = [0] * (K + 1); P[0] = 1
    for j in range(2, K + 1):
        for i in range(j, K + 1):
            P[i] += P[i - j]
    require(P == C, "recurrence check")
    return bad, weak


# =============================================================================================
def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--jobs', type=int, default=1)
    ap.add_argument('--parts', default='ALFP')
    args = ap.parse_args()
    t0 = time.time()
    allok = True
    print("Certificate for Theorem C' (q-Catalan unimodality for n >= 1001); python-flint arb, prec = %d bits"
          % ctx.prec)
    # the reduction to s <= S_TOP for k > K0: s < pi/sqrt(6 (K0+1))
    sK = PI / (6 * arb(K0 + 1)).sqrt()
    require(sK < S_TOP, "s_K < S_TOP")
    print("k > K0 = %d forces s < pi/sqrt(6(K0+1)) = %.6f <= S_TOP = %s" % (K0, float(sK.upper()), S_TOP.str(5)))
    if 'F' in args.parts or 'P' in args.parts:
        args.parts += 'A'                  # F and P use the constants A_r
    if 'A' in args.parts:
        t = time.time()
        try:
            res = part_A()
        except CertError as e:
            print("[A] FAILED:", e)
            print("NOT VERIFIED")
            return 1
        print("[A] A_r upper bounds %s (values used: %s)  %.0fs" % (
            {r: round(v, 5) for r, v in res.items()}, {r: AR[r].str(4) for r in AR}, time.time() - t))
    if 'L' in args.parts:
        t = time.time()
        try:
            bad, weak = part_L()
        except CertError as e:
            bad, weak = ["recurrence check failed"], []
        ok = not bad and not weak
        allok &= ok
        print("[L] exact, K0 = %d, 16 <= n <= K0: failures of c_{k-1} <= c_k: %s; equalities with k >= n+2: %s  %.0fs"
              % (K0, bad[:5] or "none", weak[:5] or "none", time.time() - t))
    if 'P' in args.parts:
        t = time.time()
        try:
            ok, r = check_P()
        except CertError as e:
            ok, r = False, None
            print("[P] FAILED:", e)
        allok &= ok
        if r is not None:
            print("[P] 0 < s <= %s, sigma >= 6: %s  G0 >= %.4f, far terms <= %.2e + %.2e, margin %.4f, kappa_min = %.4f, U = %g  %.0fs"
                  % (S_TOP.str(5), "verified" if ok else "FAILED", r['G0'], r['far0'], r['far1'], r['margin'], r['kmin'],
                     r['U'], time.time() - t))
    if 'F' in args.parts:
        t = time.time()
        cells = F_cells()
        if args.jobs > 1:
            from multiprocessing import Pool
            with Pool(args.jobs) as pool:
                results = pool.map(_F_worker, cells, chunksize=1)
        else:
            results = [_F_worker(c) for c in cells]
        nok = sum(1 for (_, _, ok, _) in results if ok)
        for (a, b, ok, r) in results:
            if not ok:
                print("   FAILED cell [%.6g, %.6g]: %s" % (a, b, r))
        allok &= (nok == len(cells))
        good = [r for (_, _, ok, r) in results if ok]
        if good:
            mA = min(r['G0'] - r['U0'] for r in good)
            mB = min(r['G2'] / r['KU1'] for r in good if r['okB'])
            nC = sum(1 for r in good if r['okC'])
            mC = min((r['margin_C'] for r in good if r['okC']), default=float('nan'))
            print("[F] 0 <= sigma <= 6, n >= %d: %d/%d cells verified; min(G0 - sqrt(K) I0) = %.4f, min G2/(K^1.5 I1) = %.3f,"
                  " %d cells via (1-rho)-form (min margin %.4f); max A = %.4f (< -5/N0), min rate = %.4f (> 6/N0)  %.0fs"
                  % (N0, nok, len(cells), mA, mB, nC, mC, max(r['A'] for r in good), min(r['rate'] for r in good),
                     time.time() - t))
    print("ALL VERIFIED" if allok else "NOT VERIFIED", "(%.0fs)" % (time.time() - t0))
    return 0 if allok else 1


if __name__ == "__main__":
    sys.exit(main())