#!/usr/bin/env python3
"""Beyond the band (paper "Analytic q-inequalities and a proof of Bergeron's conjecture", Sections 8.8 and 9.6): local limit bounds for Gaussian
coefficients and the criterion of Lemma 9.19 for x >= xc and 2x(x-1) < A <= 16 x^2 (Theorem 9.24 (i); the paper uses
xc = 201, so the table holds for every m >= 199; for 22 <= x <= 200 see row_sweep.py).

Gaussian coefficient G(m, n) = [m+n choose m]_q = prod_{j=1}^{m} (1 - q^{n+j}) / (1 - q^j).  For s >= 0, rho = e^{-s},
phi(theta) = G(rho e^{i theta}) / G(rho) is the characteristic function of the tilted distribution, with mean kappa_1
and variance kappa_2.  Scaled variables: sigma = m s, nu = n/m, eps = 1/m, vartheta = m theta.

TABLE (part T).  For a cell sigma in [sa, sb], nu in [na, nb], and every m >= MC:
    a2m <= kappa_2 / m^3 <= a2p,
    pl <= sqrt(2 pi kappa_2) * pp(j)  for every j with |j - kappa_1| <= doff * m   (pp = tilted point probability),
    Lu >= sqrt(2 pi kappa_2) * Lambda,  Lambda = (1/2pi) int_{-pi}^{pi} |phi|  (bounds every point probability).
Taylor zone |vartheta| <= th0 (fourth-order Taylor formula with the bounds of Lemma 9.3), Riemann zone (pair factors
of Lemma 9.11(a), summed with an explicit error valid for every m), then the Bernoulli bound of Lemma 9.21 on blocks.
CRITERION (part C).  Cells of (sigma, lam = A/x^2, psi = r/x): the evaluation bound np_bound (the sum in the proof of
Lemma 9.20) against (min pl / max Lu) sqrt(1 - loss), with pl, Lu looked up in the table.
REGIME P (part P).  sigma >= 6: the value version of regime P of Section 9.4.  Degrees D < k < D + x are covered by
Lemma 9.19 (b), so the tilt satisfies s < pi/sqrt(6 D) <= pi/sqrt(6 (2 xc - 2)).
The routines here are also used by pn_families.py (Theorem 9.24, A >= 16 x^2).
"""
import sys, time, argparse, math
from flint import arb, ctx
import cww_certificate as W
import rq_certificate as R
from cww_certificate import (A, up, lo, require, CertError, PI, ZETA2, AR, NS, b_k, _q, Psibar, sinhc, sinc, H,
                             floor3)
from rq_certificate import lam_bounds, gauss_tail

ctx.prec = 96
MC = 199                 # table valid for every m >= MC
DOFF = arb('0.5')        # offsets |j - kappa_1| <= DOFF * m
NX = 200                 # t-cells in Riemann sums
THETA_M = 50
SQ2PI = (2 * PI).sqrt()


# ---------------------------------------------------------------------------------------------
# kappa_2 bracket for G(m, n):  kappa_2 = s^-2 [sum_{j=1}^m Psi((n+j)s/2) - sum_{j=1}^m Psi(js/2)],
# Psi increasing, so each sum lies between m int_0^1 and m int_eps^{1+eps} of the scaled integrand.
# ---------------------------------------------------------------------------------------------
def a2_bounds_G(sig, nu, eps):
    """(a2m, a2p) with a2m <= kappa_2/m^3 <= a2p.  With Pb = Psibar (increasing),
        lower L = (2/sigma^3)[Pb(sigma(nu+1)/2) - Pb(sigma nu/2) - Pb(sigma(1+eps)/2) + Pb(sigma eps/2)],
        upper U = (2/sigma^3)[Pb(sigma(nu+1+eps)/2) - Pb(sigma(nu+eps)/2) - Pb(sigma/2)];
    L increases in nu and decreases in eps, U increases in nu and eps, so they are evaluated at the endpoints
    na, e1 (lower) and nb, e1 (upper); sigma enters as a ball (series in sigma^2) or through endpoint bounds of
    the brackets (closed form, where each bracket is an integral of Psi over an interval moving with sigma)."""
    sig = A(sig); na = lo(nu); nb = up(nu); e1 = up(eps)
    if (up(sig) * (nb + 1 + e1) / 2).upper() < 1.2:
        def ser(al_list):
            tot = arb(0)
            for k in range(1, NS + 1):
                b = _q(b_k[k]); p = 2 * k + 1; c = b * sig**(2 * k - 2) / arb(4)**k
                tot += c * sum((sg * a**p for (sg, a) in al_list), arb(0))
            return tot
        lo_ = ser([(1, na + 1), (-1, na), (-1, 1 + e1), (1, e1)])
        hi_ = ser([(1, nb + 1 + e1), (-1, nb + e1), (-1, arb(1))])
        # tail: |b_k| <= 2 zeta(2) pi^{-2k}; |c_k a^p| <= 2 zeta(2) a^3 pi^{-2} (sigma a/(2 pi))^{2k-2} / 4
        al = nb + 1 + e1
        q = (up(sig) * al / 2 / PI)**2
        require(q.upper() < 1, "series ratio")
        tl = 4 * 2 * ZETA2 * al**3 / PI**2 / 4 * q**NS / (1 - q)
        rl, rh = lo(lo_ - tl), up(hi_ + tl)
        require(rl.is_finite() and rh.is_finite() and rl > 0, "a2 series bracket")
        return rl, rh
    sa = lo(sig); sb = up(sig)
    # bracket integrals: int_{sigma a/2}^{sigma b/2} Psi is increasing in sigma (Psi >= 0, the interval moves right
    # and grows), so its minimum is at sa and its maximum at sb
    lo_ = 2 * (Psibar(sa * (na + 1) / 2) - Psibar(sa * na / 2) - (Psibar(sb * (1 + e1) / 2) - Psibar(sb * e1 / 2))) / sb**3
    hi_ = 2 * (Psibar(sb * (nb + 1 + e1) / 2) - Psibar(sb * (nb + e1) / 2) - Psibar(sa / 2)) / sa**3
    require(lo_.is_finite() and hi_.is_finite() and lo(lo_) > 0, "a2 lower bracket positive")
    return lo(lo_), up(hi_)


def sum_tail(q, k0):
    """sum_{k >= k0} k q^k  for 0 <= q < 1."""
    q = A(q)
    require(q.upper() < 1, "tail ratio")
    return q**k0 * (k0 / (1 - q) + q / (1 - q)**2)


def deriv_consts_G(sig, nu, th0, eps):
    """(k4z, c3): sup over |vartheta| <= th0 of |L4(s - i theta)| <= m^5 k4z, and |kappa_3| <= m^4 sigma c3.
    L4 = sum_j [((n+j)/2)^4 h4((n+j)w/2) - (j/2)^4 h4(jw/2)], h(y) = log(sinh y / y), h4 = fourth derivative of h, w = s - i theta.
    Numerator terms: h-form H(4, v), v = (1+nu) th0/2 < pi (Lemma 9.3(b)), or the decayed form, valid for every
    theta: |h4(y)| <= 3!/|y|^4 + 2^4 Li_{-3}(e^{-2 Re y}), so ((n+j)/2)^4 |h4| <= 3!(1 + A_3)/s^4.
    Denominator terms: H(4, th0/2) (needs th0/2 < pi).  Or the l-form (Lemma 9.3(a)) for sigma > 0.5."""
    sig = A(sig); nu = A(nu); eps = A(eps)
    th0 = A(th0)
    v = (1 + up(nu)) * th0 / 2
    vd = th0 / 2
    k4c = []
    if vd.upper() < float(PI.lower()) - 0.05:
        den = H(4, vd) * (1 + up(eps))**5 / 80
        numc = []
        if v.upper() < float(PI.lower()) - 0.05:
            nbu = up(nu)                                            # (nu+1+eps)^5 - nu^5 increases in nu
            numc.append(up(H(4, v) * ((nbu + 1 + up(eps))**5 - nbu**5) / 80))
        if sig.lower() > 0:
            numc.append(up(6 * (1 + AR[3]) / lo(sig)**4 * (1 + up(eps))))
        if numc:
            k4c.append(up(min(numc) + den))
    if sig.lower() > 0.5:
        k4c.append(up(AR[4] * 24 * ZETA2 / lo(sig)**5))
    nbu = up(nu); e1 = up(eps)
    c3c = [up(H(4, 0) * ((nbu + 1 + e1)**5 - nbu**5 + (1 + e1)**5) / 80)]      # S_4 at nu = nb (increasing)
    if sig.lower() > 0.5:
        c3c.append(up(AR[3] * 6 * ZETA2 / lo(sig)**5))
    if sig.lower() > 0:
        # kappa_3 = s^-3 sum_j [T((n+j)s/2) - T(js/2)], T(y) = y^3 h'''(y) = 2((y/sinh y)^3 cosh y - 1) in [-2, 0]
        c3c.append(up(arb(2) / lo(sig)**4))
    require(k4c, "no fourth-derivative bound")
    return min(k4c), min(c3c)


def local_value(sig, a2m, k4z, c3, th0, eps, doff, ucell=0.02, Ucap=12):
    """(Ilow, Iup, U) with u = theta sqrt(kappa_2) (true variance): Ilow <= (1/sqrt(2pi)) int_{|u|<=uZ}
    Re(phi~ e^{-i delta theta}) du and Iup >= (1/sqrt(2pi)) int_{|u|<=uZ} |phi~| du, for every m >= MC and
    |delta| <= doff m.  log phi~ = -u^2/2 - i kappa_3 theta^3/6 + R4, |R4| <= E4 (Lemma 9.3 bounds)."""
    sig = A(sig); eps = A(eps)
    cond = k4z * th0**2 / (24 * a2m)
    require(cond.upper() < 0.25, "Taylor-zone condition")      # |R4| <= kappa_2 theta^2/4 on the zone
    uZ = th0 * (a2m * MC).sqrt()                               # uZ = th0 sqrt(kappa_2)/m >= th0 sqrt(a2m m)
    U = min(Ucap, floor3(float(uZ.lower())))
    require(U >= 1.5 and arb(U) <= uZ, "U")
    ncell = int(U / ucell) + 1
    du = arb(U) / ncell
    e05 = eps.sqrt()
    Il = arb(0); Iu = arb(0)
    for i in range(ncell):
        u = (du * i).union(du * (i + 1))
        E4 = k4z * eps * u**4 / (24 * a2m**2)
        Y = sig * c3 * u**3 * e05 / (6 * a2m * a2m.sqrt()) + doff * u * e05 / a2m.sqrt() + 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)                            # >= int_U^uZ e^{-u^2/4} (Re log phi~ <= -u^2/4)
    Ilow = 2 * (Il - tail) / SQ2PI
    Iup = 2 * (Iu + tail) / SQ2PI
    return lo(Ilow), up(Iup), U


def _absrange(f, a, b):
    """(min, max) of |f| on [a, b] from one ball evaluation (a narrow interval of a single variable)."""
    v = f(A(a).union(A(b)))
    av = abs(v)
    return lo(av).nonnegative_part(), up(av)


def F_upper_G(ta, tb, sa, sb, tha, thb, na, nb):
    """Upper bound of F(t) = log(1 + X(nu+t)) - log(1 + X(t)), X(w) = sin^2(w vt/2)/sinh^2(w sig/2), on the box.
    Every factor is a function of one variable evaluated on a narrow interval, so no dependency blow-up:
        S0 = |sinc(t vt/2)|, H0 = sinhc(t sig/2), S1 = |sinc((nu+t) vt/2)|, H1 = sinhc((nu+t) sig/2).
    Three exact forms, the smallest bound is used (R = 1/X = (sig/vt)^2 H^2/S^2, y1 = (nu+t) sig/2):
      (1) log(1 + X1) - log(1 + X0),  X = (vt/sig)^2 S^2/H^2;
      (2) 2 log(S1 H0 / (S0 H1)) + log(1 + R1) - log(1 + R0);
      (3) 2 log(2 (y1 coth y1) H0 / ((nu+t) vt S0)) - log(1 + R0)   (uses 1 + X1 <= coth^2 y1)."""
    ta, tb, sa, sb, tha, thb, na, nb = map(A, (ta, tb, sa, sb, tha, thb, na, nb))
    S0m, S0M = _absrange(lambda z: z.sinc(), ta * tha / 2, tb * thb / 2)
    H0m, H0M = _absrange(sinhc, ta * sa / 2, tb * sb / 2)
    S1m, S1M = _absrange(lambda z: z.sinc(), (na + ta) * tha / 2, (nb + tb) * thb / 2)
    H1m, H1M = _absrange(sinhc, (na + ta) * sa / 2, (nb + tb) * sb / 2)
    require(S0m > 0, "sinc(t vt/2) vanishes on the cell")
    cands = []
    R0m = (sa / thb)**2 * H0m**2 / S0M**2
    if sa > 0:
        X1M = (thb / sa)**2 * S1M**2 / H1m**2
        X0m = (tha / sb)**2 * S0m**2 / H0M**2
        cands.append(up((1 + X1M).log() - (1 + X0m).log()))
    if S1m > 0:
        R1M = (sb / tha)**2 * H1M**2 / S1m**2
        cands.append(up(2 * (S1M * H0M / (S0m * H1m)).log() + (1 + R1M).log() - (1 + R0m).log()))
    y1b = (nb + tb) * sb / 2
    # y coth y = cosh(y) / sinhc(y), increasing in y
    ycy = up((y1b).cosh() / sinhc(y1b))
    cands.append(up(2 * (2 * ycy * H0M / ((na + ta) * tha * S0m)).log() - (1 + R0m).log()))
    best = cands[0]
    for c in cands[1:]:
        if c < best:
            best = c
    return best


def pair_cell_G(sa, sb, na, nb, tha, thb):
    """(A, B) with log|phi(theta)|^2 = sum_{j=1}^m F(j/m) <= m A + B for vartheta in [tha, thb], every m >= 1:
    M_c = max of F on t-cell c, sum_j F(j/m) <= (m/NX) sum_c M_c + sum_{c>=1} (M_c - M_{c-1})_+ (Abel summation)."""
    ts = [arb(i) / NX for i in range(NX + 1)]
    M = [F_upper_G(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
    PV = sum(((M[c] - M[c - 1]).nonnegative_part() for c in range(1, NX)), arb(0))
    return up(Aval), up(PV)


def profile_bounds(sig, nu, mc):
    """Bernoulli profile of G(m, n) = [m+n choose m]_q at tilt s = sigma/m: p_i = 1/(1 + e^{s(i-c)}), 0 <= i < M'
    = m(1+nu), sum p_i = m.  Returns (v, pm): V >= v m for every m >= mc, and max_i p_i(1-p_i) <= pm.
    sc lies in lam_bounds - [0, sigma/m] (Lemma 9.17); V >= int_0^{M'-1} - 1/4 (unimodal terms):
    V/m >= (1/(2 sigma)) [tanh(a1/2) + tanh(a2/2)] - 1/(4 mc), a1 = sc, a2 = s(M'-1) - sc; also
    V >= M' min(pi_0, pi_{M'-1}), pi = (1 - tanh^2(a/2))/4 (no cosh^2 of wide balls)."""
    sig = A(sig); nu = A(nu)
    eps_all = arb(0).union(arb(1) / mc)
    lam = lam_bounds(sig, nu)
    require(lam.is_finite(), "Fermi level not finite")
    a1 = lam - sig * eps_all
    a2 = sig * (1 + nu) - sig * eps_all - a1
    w1 = (1 - (a1 / 2).tanh()**2) / 4; w2 = (1 - (a2 / 2).tanh()**2) / 4
    require(w1.is_finite() and w2.is_finite(), "profile weights not finite")
    v = (1 + lo(nu)) * min(lo(w1), lo(w2))
    if sig.lower() > 1e-4:
        v2 = lo(((a1 / 2).tanh() + (a2 / 2).tanh()) / (2 * sig) - arb(1) / (4 * mc))
        require(v2.is_finite(), "V bound not finite")
        if v2 > v:
            v = v2
    top = up(lam)
    pm = arb(1) / 4 if top >= 0 else up((1 - (top / 2).tanh()**2) / 4)
    return lo(v), pm


def abel_expo(v, pm, th, x):
    """upper bound of log[(|phi(theta)| / pref)] = -V + |Z| <= -v m + pm / sin(theta/2), at m = x."""
    return -v * x + pm / (A(th) / 2).sin()


def far_zone(sa, sb, na, nb, th0, mc, a2p):
    """(F, rateA, rateB): F >= sqrt(kappa_2/(2 pi)) * 2 int_{th0/m}^{pi} |phi| for every m >= mc.
    Riemann cells (pair products, Lemma 9.11(a)) while the Abel exponent is weak, then Bernoulli-Abel blocks:
    |phi(theta)| <= (4/3)(4 sqrt V + 1) exp(-V + pm / sin(theta/2)) (Bernoulli representation, Abel summation)."""
    sa = A(sa); sb = A(sb); na = A(na); nb = A(nb)
    sig = sa.union(sb); nu = na.union(nb)
    x = arb(mc)
    v, pm = profile_bounds(sig, nu, mc)
    VM = (1 + up(nu) + 1 / x) * x / 4
    pref = arb(4) / 3 * (4 * VM.sqrt() + 1)
    tot = arb(0); rateA = None; rateB = None
    th = A(th0)
    # Riemann zone: until the Abel exponent per m is below -max(0.05, 4.5/mc) (and at most up to vartheta = 8); the
    # Bernoulli blocks then have rate * mc > 4, which the m-monotonicity check below requires
    e_sw = max(0.05, 4.5 / mc)
    while True:
        e_here = 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_G(sa, sb, na, nb, 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
    # Bernoulli-Abel blocks
    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 = 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 = abel_expo(v, pm, a_, x)                              # the integrand decreases in theta
        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 taylor_candidates(sa, sb, nb):
    """candidate Taylor-zone ends th0 (m-scaled)"""
    c = []
    for v in ('1.2', '1.0', '0.8', '0.6'):
        c.append(2 * arb(v) / (1 + arb(nb)))
    for f in ('0.1', '0.15', '0.2', '0.3', '0.5', '0.8', '1.2', '1.8', '2.6', '3.6', '5.0'):
        c.append(arb(f))
    if sa > 0.5:
        for f in ('0.2', '0.3', '0.45', '0.6'):
            c.append(arb(f) * sa)
    return c


def doff_cell(sa, sb, na, nb):
    """offset allowance (in units of m) for G(x, D) at the saddle point of G(x, D-1): Delta/x =
    (1/sigma)[B(nu sigma) - B((1+nu) sigma)], B(u) = u/(e^u - 1) decreasing with B' in [-1/2, 0], so <= 1/2."""
    sa = A(sa); sb = A(sb); na = A(na); nb = A(nb)
    if sa.lower() <= 0:
        return arb(1) / 2
    Bf = lambda u: A(u) / A(u).expm1()
    d = (Bf(na * sa) - Bf((1 + nb) * sb)) / sa
    require(d.is_finite(), "offset bound")
    return min(arb(1) / 2, up(d))


def table_cell(sa, sb, na, nb, mc=MC, doff=None):
    """LLT bounds on the cell, valid for every m >= mc: plc <= sqrt(2 pi k2) pp (offsets up to doff m),
    Luc >= sqrt(2 pi k2) Lambda."""
    sa = A(sa); sb = A(sb); na = A(na); nb = A(nb)
    if doff is None:
        doff = doff_cell(sa, sb, na, nb)
    sig = sa.union(sb); nu = na.union(nb)
    eps = arb(0).union(arb(1) / mc).nonnegative_part()
    a2m, a2p = a2_bounds_G(sig, nu, eps)
    best = None; why = {}
    for t0 in taylor_candidates(float(sa.lower()), float(sb.upper()), float(nb.upper())):
        try:
            k4z, c3 = deriv_consts_G(sig, nu, t0, eps)
            Il, Iu, U = local_value(sig, a2m, k4z, c3, t0, eps, doff)
            F, rateA, rateB = far_zone(sa, sb, na, nb, t0, mc, a2p)
            # m-monotonicity: Riemann part ~ m^{1/2} e^{m A/2} (A < -2/mc); Abel part ~ m^2 e^{-rate m}
            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, Il, Iu, F, U, rateB)
    require(best is not None, "no admissible Taylor zone: %s" % why)
    score, th0, plc, Luc, Il, Iu, F, U, rate = best
    return dict(sig=(float(sa.lower()), float(sb.upper())), nu=(float(na.lower()), float(nb.upper())),
                a2m=a2m, a2p=a2p, pl=plc, Lu=Luc, Il=float(Il.lower()), Iu=float(Iu.upper()), far=float(F.upper()),
                U=U, rate_mc=float((rate * mc).lower()), th0=float(th0.lower()), score=score)


# =============================================================================================
# Evaluation bound:  N_A(rho) / P_A(rho) <= Q (1 - psi + F(t)) / B   (see the plan)
# lam = A/x^2, psi = r/x, eps = 1/x, s = sigma/x, D = lam x - psi, W = A/(x(x-1)) = lam/(1-eps), t = W - 1 + psi.
# =============================================================================================
def _Epoint(z):
    z = A(z)
    return (-z / 2).exp() * sinhc(z / 2)


def Efun(z):
    """E(z) = (1 - e^{-z})/z = int_0^1 e^{-zt} dt (E(0) = 1), decreasing in z: evaluated at the two endpoints of the
    ball, so that wide arguments do not produce sign-indefinite balls."""
    z = A(z)
    a = _Epoint(lo(z)); b = _Epoint(up(z))
    r = b.union(a)
    require(r.is_finite() and r.lower() > 0, "E: ball")
    return r


def np_bound(sig, lam, psi, eps):
    """Upper bound for N(rho)/P(rho) on the ball (sigma, lam, psi, eps): all factors in closed form."""
    sig = A(sig); lam = A(lam); psi = A(psi); eps = A(eps)
    Dx = lam - psi * eps                                   # D/x
    Q = ((1 - eps) / (Dx * (Dx + eps)) * Efun(sig) * Efun(sig * (1 - eps))
         / (Efun(sig * Dx) * Efun(sig * (Dx + eps))) * (-sig * (Dx + eps - 1)).exp())
    z = (-sig * eps).exp() * (Dx + 1) / (Dx + 2 * eps) * Efun(sig * (Dx + 1)) / Efun(sig * (Dx + 2 * eps))
    Rp = (Dx + 1) / Dx * Efun(sig * (Dx + 1)) / Efun(sig * Dx)      # [x+D]/[D]
    B = 1 + lo(psi) * (lo(Rp) - 1)                                  # (1-psi) + psi R', increasing in psi (R' >= 1)
    t = lam / (1 - eps) - 1 + psi
    tu = up(t)
    fl = int(float(tu.upper()) // 1)
    while arb(fl) > tu:                                   # exact floor of tu (float rounding could overshoot)
        fl -= 1
    while arb(fl + 1) <= tu:
        fl += 1
    zu = up(z)
    F = sum((zu**(u - 1) for u in range(2, fl + 2)), arb(0)) + (tu - fl) * zu**(fl + 1)
    # S(psi) = (1 - psi) + F(W - 1 + psi) increases in psi when z >= 1 (F' = z^{floor t + 1} >= 1)
    first = (1 - up(psi)) if lo(z) >= 1 else (1 - lo(psi))
    res = up(up(Q) * (first + F) / lo(B))
    require(res.is_finite() and lo(B) > 0, "np_bound not finite")
    return res


# =============================================================================================
# Drivers: the table (with adaptive splitting) and the criterion over (sigma, lam, psi) cells.
# =============================================================================================
def source_digest(names):
    """short sha256 digests of the given source files (next to this one), for the log header."""
    import hashlib, os
    here = os.path.dirname(os.path.abspath(__file__))
    out = []
    for nm in names:
        with open(os.path.join(here, nm), 'rb') as f:
            out.append("%s %s" % (nm, hashlib.sha256(f.read()).hexdigest()[:16]))
    return ", ".join(out)


def _table_worker(args):
    global MC
    sa, sb, na, nb, mc = args
    MC = mc                                        # local_value reads MC; set it here too (not only in the parent)
    t = time.time()
    try:
        r = table_cell(arb(sa), arb(sb), arb(na), arb(nb), mc=mc)
        out = dict(sa=sa, sb=sb, na=na, nb=nb, ok=True, pl=math.nextafter(float(r['pl'].lower()), -math.inf),
                   Lu=math.nextafter(float(r['Lu'].upper()), math.inf),
                   score=r['score'], far=r['far'], U=r['U'], th0=r['th0'])
    except CertError as e:
        out = dict(sa=sa, sb=sb, na=na, nb=nb, ok=False, why=str(e)[:200])
    out['time'] = time.time() - t
    return out


def build_table(sig_edges, nu_edges, mc, jobs=1, max_depth=3, log=None):
    """leaf cells covering [sig_edges[0], sig_edges[-1]] x [nu_edges[0], nu_edges[-1]]; a failing cell is split
    into 2 x 2 (exact binary endpoints) up to max_depth times."""
    from multiprocessing import Pool
    todo = [(sig_edges[i], sig_edges[i + 1], nu_edges[j], nu_edges[j + 1], mc)
            for i in range(len(sig_edges) - 1) for j in range(len(nu_edges) - 1)]
    leaves = []; failed = []
    pool = Pool(jobs) if jobs > 1 else None
    for depth in range(max_depth + 1):
        res = pool.map(_table_worker, todo, chunksize=4) if pool else [_table_worker(c) for c in todo]
        todo = []
        for r in res:
            if r['ok']:
                leaves.append(r)
            elif depth < max_depth:
                sa, sb, na, nb = r['sa'], r['sb'], r['na'], r['nb']
                sm = (sa + sb) / 2; nm = (na + nb) / 2
                todo += [(a, b, c, d, mc) for (a, b) in ((sa, sm), (sm, sb)) for (c, d) in ((na, nm), (nm, nb))]
            else:
                failed.append(r)
        if log:
            print("    table depth %d: %d leaves, %d to split, %d failed" % (depth, len(leaves), len(todo), len(failed)),
                  file=log, flush=True)
        if not todo:
            break
    if pool:
        pool.close()
    return leaves, failed


def lookup(leaves, s0, s1, n0, n1, key, worst):
    """worst (min or max) of leaf[key] over the leaves meeting [s0, s1] x [n0, n1]; None if not covered."""
    vals = [lf[key] for lf in leaves if lf['sa'] <= s1 and lf['sb'] >= s0 and lf['na'] <= n1 and lf['nb'] >= n0]
    return worst(vals) if vals else None


def nuN_hi(lb, xc):
    """upper end of nu_e = e/(x-2) for the N-products G(x-2, e), D < e <= E+1, when lam <= lb and x >= xc:
    e <= E + 1 <= lam x^2/(x-1) + 1, so nu_e <= lam x^2/((x-1)(x-2)) + 1/(x-2); both terms decrease in x.
    Rounded up (the query ranges are padded outward as well)."""
    v = lb * xc * xc / ((xc - 1) * (xc - 2)) + 1.0 / (xc - 2)
    return math.nextafter(math.nextafter(v, math.inf), math.inf)


def criterion_cell(leaves, sa, sb, la, lb, pa, pb, xc):
    """Remark 8.25 criterion on sigma in [sa, sb], lam in [la, lb], psi in [pa, pb], every x >= xc:
    N/P < (min pl / max Lu) sqrt(1 - loss).  Returns (ok, margin, info)."""
    e = 1.0 / xc
    np_ = np_bound(arb(sa).union(arb(sb)), arb(la).union(arb(lb)), arb(pa).union(arb(pb)),
                   arb(0).union(arb(1) / xc))
    # P-products G(x, D-1+i): m = x, sigma, nu in [la - 2/xc, lb]; N-products G(x-2, e): m = x - 2,
    # sigma_N = (x-2) s in [sa (1 - 2/xc), sb], nu_e in [la, nuN_hi(lb, xc)]
    pl = lookup(leaves, sa, sb, la - 2 * e, lb, 'pl', min)
    Lu = lookup(leaves, sa * (1 - 2 * e), sb, la, nuN_hi(lb, xc), 'Lu', max)
    if pl is None or Lu is None:
        return False, None, "table does not cover the cell"
    # kappa_2(G(x-2, e)) >= kappa_2(G(x, D)) (1 - loss), loss <= min(2/sigma^2, (1+lam)^2/6) / (x a2(G(x, D)))
    sig = arb(sa).union(arb(sb))
    a2m, _ = a2_bounds_G(sig, arb(la - 2 * e).union(arb(lb)), arb(0).union(arb(1) / xc))
    c = (1 + arb(lb))**2 / 6
    if sa > 0:
        c = min(c, up(arb(2) / arb(sa)**2))
    loss = up(c / (xc * a2m))
    require(loss < 1, "variance loss")
    rhs = arb(pl) / arb(Lu) * (1 - loss).sqrt()
    margin = float((rhs - np_).lower())
    return margin > 0, margin, dict(np=float(np_.upper()), pl=pl, Lu=Lu, loss=float(loss.upper()))


class LeafIndex:
    """leaves grouped by the sigma-strip of the base grid (a leaf lies in exactly one strip)."""
    def __init__(self, leaves, sig_edges):
        import bisect
        self.edges = sig_edges
        self.strips = [[] for _ in range(len(sig_edges) - 1)]
        for lf in leaves:
            i = bisect.bisect_right(sig_edges, lf['sa']) - 1
            i = min(max(i, 0), len(self.strips) - 1)
            self.strips[i].append(lf)

    def query(self, s0, s1, n0, n1, key, worst):
        import bisect
        i0 = max(bisect.bisect_right(self.edges, s0) - 2, 0)
        i1 = min(bisect.bisect_right(self.edges, s1), len(self.strips) - 1)
        vals = []
        covered_s = []
        for i in range(i0, i1 + 1):
            for lf in self.strips[i]:
                if lf['sa'] <= s1 and lf['sb'] >= s0 and lf['na'] <= n1 and lf['nb'] >= n0:
                    vals.append(lf[key])
        return worst(vals) if vals else None


def _pad(v, sgn):
    return v + sgn * 1e-12 * max(1.0, abs(v))        # query ranges padded outward against float rounding


def criterion_cell_idx(idx, sa, sb, la, lb, pa, pb, xc):
    e = 1.0 / xc
    pl = idx.query(_pad(sa, -1), _pad(sb, 1), _pad(la - 2 * e, -1), _pad(lb, 1), 'pl', min)
    Lu = idx.query(_pad(sa * (1 - 2 * e), -1), _pad(sb, 1), _pad(la, -1), _pad(nuN_hi(lb, xc), 1), 'Lu', max)
    if pl is None or Lu is None:
        return False, None, "table does not cover the cell"
    np_ = np_bound(arb(sa).union(arb(sb)), arb(la).union(arb(lb)), arb(pa).union(arb(pb)), arb(0).union(arb(1) / xc))
    sig = arb(sa).union(arb(sb))
    a2m, _ = a2_bounds_G(sig, arb(la - 2 * e).union(arb(lb)), arb(0).union(arb(1) / xc))
    c = (1 + arb(lb))**2 / 6
    if sa > 0:
        c = min(c, up(arb(2) / arb(sa)**2))
    loss = up(c / (xc * a2m))
    if not (loss < 1):
        return False, None, "variance loss"
    rhs = arb(pl) / arb(Lu) * (1 - loss).sqrt()
    margin = float((rhs - np_).lower())
    return margin > 0, margin, dict(np=float(np_.upper()), pl=pl, Lu=Lu, loss=float(loss.upper()))


def run_criterion(idx, sig_edges, lam_edges, psi_edges, xc, max_depth=8, log=None):
    """every (sigma, lam, psi) box of the grid, refined on failure: lam and psi halved, and sigma halved as well
    from the second refinement on (every sub-box still lies in the union of the table leaves it queries)."""
    todo = [(sig_edges[i], sig_edges[i + 1], lam_edges[j], lam_edges[j + 1], psi_edges[k], psi_edges[k + 1])
            for i in range(len(sig_edges) - 1) for j in range(len(lam_edges) - 1) for k in range(len(psi_edges) - 1)]
    nok = 0; failed = []; minmargin = None
    for depth in range(max_depth + 1):
        nxt = []
        for (sa, sb, la, lb, pa, pb) in todo:
            try:
                ok, mg, info = criterion_cell_idx(idx, sa, sb, la, lb, pa, pb, xc)
            except CertError as ex:
                ok, mg, info = False, None, str(ex)
            if ok:
                nok += 1
                minmargin = mg if minmargin is None or mg < minmargin else minmargin
            elif depth < max_depth and info != "table does not cover the cell":
                lm = (la + lb) / 2; pm_ = (pa + pb) / 2
                sgs = ((sa, sb),) if depth == 0 else ((sa, (sa + sb) / 2), ((sa + sb) / 2, sb))
                nxt += [(s0, s1, a, b, c, d) for (s0, s1) in sgs for (a, b) in ((la, lm), (lm, lb))
                        for (c, d) in ((pa, pm_), (pm_, pb))]
            else:
                failed.append(((sa, sb, la, lb, pa, pb), info))
        if log:
            print("    criterion depth %d: %d ok, %d to refine, %d failed" % (depth, nok, len(nxt), len(failed)),
                  file=log, flush=True)
        todo = nxt
        if not todo:
            break
    return nok, failed, minmargin


# =============================================================================================
# Regime P: sigma = m s >= SIGP0 and n s >= NUMIN * SIGP0, 0 < s <= sb (value version of Section 9.4, [P]).
# s-scaled: s^3 kappa_2 in [k2m, pi^2/3]; l-form derivative bounds |L^(r)(s - i t)| <= A_r r! zeta(2)/s^{r+1}.
# =============================================================================================
SIGP0 = arb('5.85')      # sigma_N = (x-2) s >= 6 (1 - 2/x) >= SIGP0 for x >= 81 when sigma = x s >= 6
NUMIN = arb('1.96')      # n s >= (lam - 2/x) sigma >= (2 - 4/x) sigma >= NUMIN sigma for x >= 100
LAM3 = AR[3] * 6 * ZETA2; LAM4 = AR[4] * 24 * ZETA2
K2P = PI**2 / 3


def Gtail(a):
    """>= int_a^oo y^2 / (4 sinh^2(y/2)) dy  (a >= 1)."""
    a = A(a)
    return (a**2 + 2 * a + 2) * (-a).exp() / (1 - (-a).exp())**2


def local_P_value(sb, gam, ucell=0.02, Ucap=12):
    """(Ilow, Iup, U, k2m) for the Taylor zone theta <= gam s, every s in (0, sb], u = theta sqrt(kappa_2)."""
    s = up(sb)
    k2m = K2P - s - Gtail(SIGP0) - Gtail(NUMIN * SIGP0)
    cond = LAM4 * gam**2 / (24 * k2m)
    require(cond.upper() < 0.25, "P: Taylor-zone condition")
    uZ = gam * (k2m / s).sqrt()
    U = min(Ucap, floor3(float(uZ.lower())))
    require(U >= 3 and arb(U) <= uZ, "P: U")
    ncell = int(U / ucell) + 1
    du = arb(U) / ncell
    s05 = s.sqrt(); s15 = s * s05
    # offset for G(x, D): Delta <= beta(D) <= (1/s) y e^{-y}/(1 - e^{-y}), y = D s >= NUMIN SIGP0 (y e^{-y} decreases)
    y0 = NUMIN * SIGP0
    dl = y0 * (-y0).exp() / (1 - (-y0).exp())                    # s * Delta <= dl
    Il = arb(0); Iu = arb(0)
    for i in range(ncell):
        u = (du * i).union(du * (i + 1))
        E4 = LAM4 * u**4 * s / (24 * k2m**2)
        Y = LAM3 * u**3 * s05 / (6 * k2m * k2m.sqrt()) + dl * u * s05 / k2m.sqrt() + 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, k2m


def far_P_value(sb, gam, cx=arb(20), ratio=arb('1.1'), ny=240):
    """zones for log|phi(theta)| <= -kappa/s on theta >= gam s (Lemma 9.11(b) with parts j >= 1):
    denominators sum_j log(1 + sin^2(j th/2)/sinh^2(j s/2)) >= (1/s) int_s^Y ..., numerators <= NUM/s."""
    sb = up(sb)
    y0 = NUMIN * SIGP0 - sb                                       # (n+1) s - s >= n s - s
    NUMs = 4 * (-y0).exp() / (1 - (-y0).exp())                    # s * sum_j 1/sinh^2((n+j)s/2) <= NUMs
    zone1 = []
    c = A(gam)
    while c.upper() < cx.upper():
        cb = c * ratio
        if cb.upper() > cx.upper():
            cb = cx
        Y = min(lo(PI / cb), lo(SIGP0))
        J = W.lower_riemann(lambda yy: (1 + (c * yy / 2).sin()**2 / (yy / 2).sinh()**2).log(), sb, Y, ny)
        zone1.append((c, cb, lo(J / 2 - NUMs / 2)))
        c = cb
    SIN2 = (PI / 8).sin()**2
    Bg = W.lower_riemann(lambda yy: (1 + SIN2 / (yy / 2).sinh()**2).log(), up(PI / cx + sb),
                         lo(SIGP0 - 2 * PI / cx - sb), 400, geometric=True)
    kap2 = lo(Bg / 4 - NUMs / 2)
    return zone1, kap2


def check_P_value(sb, gam=arb('0.55')):
    """(pl, Lu) for regime P, every s in (0, sb]: pl <= sqrt(2 pi kappa_2) pp (with the G(x,D) offset), Lu >=
    sqrt(2 pi kappa_2) Lambda."""
    Il, Iu, U, k2m = local_P_value(sb, gam)
    zone1, kap2 = far_P_value(sb, gam)
    ks = [k for (_, _, k) in zone1] + [kap2]
    kmin = ks[0]
    for k in ks[1:]:
        if k < kmin:
            kmin = k
    s = up(sb)
    # far terms ~ s^{-1/2} e^{-k/s} (zone 1) and s^{-3/2} e^{-k/s} (zone 2) increase with s when k/s >= 3/2
    require((kmin / s).lower() > 1.5, "P: kappa_min / s > 3/2")
    Ek = lambda k: (-k / s).exp()
    S0 = sum(((cb - c) * s * Ek(k) for (c, cb, k) in zone1), arb(0)) + PI * Ek(kap2)
    far = (K2P / (2 * PI * s**3)).sqrt() * 2 * S0                  # sqrt(kappa_2/(2 pi)) * 2 int |phi|
    return lo(Il - far), up(Iu + far), dict(Il=float(Il.lower()), Iu=float(Iu.upper()), far=float(far.upper()),
                                            U=U, kmin=float(kmin.lower()))


def V_P(sb):
    """regime P (m s >= SIGP0, n s >= NUMIN SIGP0): V >= 0.99/s - 1/4 from the Fermi level bounds
    m - 2 - 1/s <= c <= m + 1 + 1/s, so s c >= SIGP0 - 2s - 1 and s(M'-1-c) >= NUMIN SIGP0 - 2s - 1; V <= 2/s + 2."""
    s = up(sb)
    a1 = SIGP0 - 2 * s - 1
    a2 = NUMIN * SIGP0 - 2 * s - 1
    require(a1.lower() > 0 and a2.lower() > 0, "V_P")
    vcoef = lo(((a1 / 2).tanh() + (a2 / 2).tanh()) / 2)          # V >= vcoef/s - 1/4
    return vcoef


def check_P_cell(sa, sb, gam, cx=arb(3), ratio=arb('1.1'), ny=240):
    """(pl, Lu, info) for every s in [sa, sb] (sa > 0), or every s in (0, sb] if sa == 0 (then the far terms are
    shown to increase with s, so s = sb is the worst case)."""
    sa = A(sa); sb = up(sb)
    Il, Iu, U, k2m = local_P_value(sb, gam)
    y0 = NUMIN * SIGP0 - sb
    NUMs = 4 * (-y0).exp() / (1 - (-y0).exp())
    # zone 1: theta = c s, c in [gam, cx]: log|phi| <= -kappa_c / s
    zone1 = []
    c = A(gam)
    while c.upper() < cx.upper():
        cb = c * ratio
        if cb.upper() > cx.upper():
            cb = cx
        Y = min(lo(PI / cb), lo(SIGP0))
        J = W.lower_riemann(lambda yy: (1 + (c * yy / 2).sin()**2 / (yy / 2).sinh()**2).log(), sb, Y, ny)
        zone1.append((c, cb, lo(J / 2 - NUMs / 2)))
        c = cb
    kmin = min((k for (_, _, k) in zone1), key=lambda k: float(k.lower()))
    vcoef = V_P(sb)
    # zone 2: theta in [cx s, pi]: Bernoulli-Abel, |phi| <= pref exp(-V + 1/(4 sin(theta/2)))
    #   1/(4 sin(theta/2)) <= 1/(2 theta (1 - theta^2/24)) on [cx s, 2], <= 1/(4 sin 1) on [2, pi]
    def far2(s):
        s = A(s)
        pref = arb(4) / 3 * (4 * (2 / s + 2).sqrt() + 1)
        tot = arb(0)
        th = cx * s
        while th.upper() < 2:
            tb = th * arb('1.1')
            if tb.upper() > 2: tb = arb(2)
            tot += (tb - th) * (1 / (2 * th * (1 - th**2 / 24))).exp()
            th = tb
        tot += (PI - 2) * (1 / (4 * arb(1).sin())).exp()
        return (K2P / (2 * PI * s**3)).sqrt() * 2 * pref * (-vcoef / s + arb(1) / 4).exp() * tot
    if sa.lower() > 0:
        # s in [sa, sb]: zone-1 terms bounded at the worst ends, zone 2 evaluated on the ball
        sball = sa.union(sb)
        S0 = sum(((cb - c) * sb * (-k / sb).exp() for (c, cb, k) in zone1), arb(0))
        far = up((K2P / (2 * PI * sa**3)).sqrt() * 2 * S0) + up(far2(sball))
    else:
        # (0, sb]: each far term ~ s^{-1/2} e^{-k/s} or s^{-2} e^{-b/s} increases with s when k/s >= 1/2, b/s >= 2
        require((kmin / sb).lower() > 0.5, "P: kappa_min/s")
        b2 = vcoef - 1 / (2 * cx * (1 - (cx * sb)**2 / 24))        # zone-2 exponent at theta = cx s is <= -b2/s
        require(b2 / sb > 2, "P: zone-2 exponent")
        S0 = sum(((cb - c) * sb * (-k / sb).exp() for (c, cb, k) in zone1), arb(0))
        # zone 2 in closed form: exp(1/(4 sin(theta/2))) <= exp(1/(2 cx s (1 - (cx sb)^2/24))) on [cx s, pi], so the
        # zone-2 term is <= sqrt(K2P/(2 pi s^3)) 2 pref(s) pi e^{1/4} e^{-b2/s}, which increases in s when b2/s > 2
        pref_b = arb(4) / 3 * (4 * (2 / sb + 2).sqrt() + 1)
        far2c = (K2P / (2 * PI * sb**3)).sqrt() * 2 * pref_b * PI * arb('0.25').exp() * (-b2 / sb).exp()
        far = up((K2P / (2 * PI * sb**3)).sqrt() * 2 * S0) + up(far2c)
    return lo(Il - far), up(Iu + far), dict(Il=float(Il.lower()), Iu=float(Iu.upper()), far=float(far.upper()), U=U,
                                           kmin=float(kmin.lower()), vcoef=float(vcoef))


# =============================================================================================
def np_sup_P(la, lb, pa, pb, xc, sig_top=300):
    """sup of the N/P bound over sigma >= 6 (cells up to sig_top, then the tail bound)."""
    best = arb(0)
    v = 6.0
    while v < sig_top:
        w = min(v * 1.1, sig_top)
        b = np_bound(arb(v).union(arb(w)), arb(la).union(arb(lb)), arb(pa).union(arb(pb)), arb(0).union(arb(1) / xc))
        require(b.is_finite(), "np_sup_P")
        if b > best: best = up(b)
        v = w
    # sigma >= sig_top: Q <= e^{-sigma(D/x + eps - 1)} / ((1 - e^{-sigma D/x})(1 - e^{-sigma(D/x+eps)})) and
    # S <= (1 - psi) + sum_{u=2}^{floor(t)+2} z^{u-1} with z <= 1 + 4/sigma... bounded crudely by (t + 2) 2^{t+1}
    Dx = arb(la) - arb(1) / xc
    t = arb(lb) / (1 - arb(1) / xc) + 1
    tail = (-arb(sig_top) * (Dx - 1)).exp() / (1 - (-arb(sig_top) * Dx).exp())**2 * (t + 2) * arb(2)**(t + 1)
    require(tail.is_finite(), "np_sup_P tail")
    if tail > best: best = up(tail)
    return best


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--jobs', type=int, default=1)
    ap.add_argument('--xc', type=int, default=201)
    ap.add_argument('--lammax', type=float, default=16.0)
    ap.add_argument('--parts', default='TCP')
    ap.add_argument('--test', default=None)
    ap.add_argument('--table', default=None, help="part C from a table saved by part T (json) instead of rebuilding it")
    args = ap.parse_args()
    global MC
    MC = args.xc - 2                                          # local_value uses MC: the table holds for every m >= MC
    if args.test:
        sa, sb, na, nb = (float(v) for v in args.test.split(','))
        t = time.time()
        r = table_cell(arb(sa), arb(sb), arb(na), arb(nb), mc=args.xc - 2)
        print({k: (float(v.lower()) if isinstance(v, arb) else v) for k, v in r.items()}, "%.1fs" % (time.time() - t))
        return 0
    xc = args.xc; mc = xc - 2; e = 1.0 / xc
    require(6 * (1 - arb(2) / xc) >= SIGP0 and 2 - arb(4) / xc >= NUMIN, "regime P constants need a larger xc")
    log = sys.stdout
    t0 = time.time()
    allok = True
    print("sources: " + source_digest(['pn_certificate.py', 'rq_certificate.py', 'cww_certificate.py']))
    print("Theorem 9.24 (i): R(x, A) for x >= %d, 2x(x-1) < A <= %.1f x^2; arb prec = %d"
          % (xc, args.lammax, ctx.prec), flush=True)
    lam0 = 2 - 2 * e                                          # A > 2x(x-1) and x >= xc give lam > 2 - 2/xc
    sig_edges = [0.0, 0.001]
    v = 0.001
    while v < 6:
        w = min(v * (1.25 if v < 0.05 else 1.1), 6.0); sig_edges.append(w); v = w
    nu_edges = [lam0 - 2 * e]
    v = nu_edges[0]
    numax = nuN_hi(args.lammax, xc)
    while v < numax:
        w = min(v * 1.05, numax); nu_edges.append(w); v = w
    if args.table and 'T' not in args.parts and 'C' in args.parts:
        import json
        leaves = json.load(open(args.table))
        # the stored values are already rounded outward; round once more (harmless, and covers older files)
        for lf in leaves:
            lf['pl'] = math.nextafter(lf['pl'], -math.inf); lf['Lu'] = math.nextafter(lf['Lu'], math.inf)
        require(all(lf['ok'] for lf in leaves), "saved table has failed cells")
        # the leaves must tile the base grid: total area check
        area = sum((lf['sb'] - lf['sa']) * (lf['nb'] - lf['na']) for lf in leaves)
        full = (sig_edges[-1] - sig_edges[0]) * (nu_edges[-1] - nu_edges[0])
        require(abs(area - full) < 1e-9 * full and min(lf['sa'] for lf in leaves) <= sig_edges[0]
                and max(lf['sb'] for lf in leaves) >= sig_edges[-1] and min(lf['na'] for lf in leaves) <= nu_edges[0]
                and max(lf['nb'] for lf in leaves) >= nu_edges[-1], "saved table does not tile the domain")
        print("[T] table loaded from %s: %d leaves, min score %.3f" % (args.table, len(leaves),
                                                                     min(lf['score'] for lf in leaves)), flush=True)
    elif 'T' in args.parts or 'C' in args.parts:
        t = time.time()
        leaves, failed = build_table(sig_edges, nu_edges, mc, jobs=args.jobs, log=log)
        print("[T] table: %d sigma x %d nu base cells -> %d leaves, %d failed; min score %.3f  %.0fs"
              % (len(sig_edges) - 1, len(nu_edges) - 1, len(leaves), len(failed),
                 min(lf['score'] for lf in leaves), time.time() - t), flush=True)
        for f in failed[:10]:
            print("   failed:", f, flush=True)
        allok &= not failed
        import json
        json.dump(leaves, open('pn_table_%d.json' % xc, 'w'))
    if 'C' in args.parts:
        t = time.time()
        idx = LeafIndex(leaves, sig_edges)
        lam_edges = [lam0]
        v = lam0
        while v < args.lammax:
            w = min(v * 1.05, args.lammax); lam_edges.append(w); v = w
        psi_edges = [i / 10 for i in range(11)]
        nok, cfailed, mm = run_criterion(idx, sig_edges, lam_edges, psi_edges, xc, log=log)
        print("[C] criterion, sigma <= 6: %d boxes verified, %d failed; min margin %.4f  %.0fs"
              % (nok, len(cfailed), mm if mm is not None else float('nan'), time.time() - t), flush=True)
        for f in cfailed[:10]:
            print("   failed:", f, flush=True)
        allok &= not cfailed
    if 'P' in args.parts:
        t = time.time()
        # degrees D < k < D + x are covered by the low-degree lemma; k >= D + x means that the tilted mean k - x of
        # G(x, D-1) is >= D >= 2x - 2, so s < pi/sqrt(6 D) <= pi/sqrt(6 (2 xc - 2))
        smax = PI / (6 * arb(2 * xc - 2)).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:
            ok = False
            for g in ('0.5', '0.55', '0.6', '0.45'):
                try:
                    pl, Lu, info = check_P_cell(a, b, arb(g))
                    ok = True
                    break
                except CertError:
                    continue
            require(ok, "regime P cell failed")
            r = pl / Lu
            require(r.is_finite(), "regime P ratio")
            worstP = r if worstP is None or r < worstP else worstP
        # criterion: sup_{sigma >= 6} N/P < worstP * sqrt(1 - loss), loss <= 2 s_max / k2m
        lossP = up(2 * smax / (K2P - smax - Gtail(SIGP0) - Gtail(NUMIN * SIGP0)))
        rhs = lo(worstP) * (1 - lossP).sqrt()
        lam_edges = [lam0]
        v = lam0
        while v < args.lammax:
            w = min(v * 1.1, args.lammax); lam_edges.append(w); v = w
        worst_np = arb(0)
        for j in range(len(lam_edges) - 1):
            for k in range(10):
                b = np_sup_P(lam_edges[j], lam_edges[j + 1], k / 10, (k + 1) / 10, xc)
                if b > worst_np: worst_np = b
        okP = bool(worst_np < rhs)
        allok &= okP
        print("[P] sigma >= 6: %d s-cells, min pl/Lu = %.3f; sup N/P = %.3g < %.3f: %s  %.0fs"
              % (len(cells), float(worstP.lower()), float(worst_np.upper()), float(rhs.lower()),
                 "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())