#!/usr/bin/env python3
r"""row33_eps0.py (vendored) -- row 33 eps^0 layer as a SUBTRACTED FINITE INTEGRAL.

I(eps) = c_-2/eps^2 + c_-1/eps + c0 + O(eps) at (s,t,m2)=(-1,-1/3,1):
  c_-2 = -K2^(0)
  c_-1 = -K2^(1) + 2 K2^(0)(gamma_E + lnL) + J0 + Tm1
  c0   = -[ghat(e) K2(e) L^{-2e}]_2 + ghat1*J0 + J1 + B0 + T0
No eps grid anywhere: K2 jets = exact parametric 2-folds (GL n-doubling until
bound), kernel jet K^(1)(w) = 3-layer block jet transport of the exact box DE
(seeded by analytic Laurent jets of the Feynman-parameter masters), rho jets =
closed 2-body form + Frobenius-jet series + 2-layer self-energy block jet
transport, finite integrals via detransport.quad.quad_refine, all fail-closed.
The recorded reference strings (pole layers, eps^0) are comparisons only,
never inputs.

This file is the research engine, vendored statement for statement; comments
and message strings were edited for release, and two things were added after
the first release (every other unit is the research code unchanged):

  * step controller (JetMarch._march): a rejected step re-sums the SAME local
    Taylor series at the halved step instead of recomputing it -- RatDE.step's
    coefficients depend on the start point and the order N only, never on h,
    so the accepted steps, their tail bounds and every printed value are the
    same as before; only the wasted recomputations (one per halving, i.e. one
    per accepted step at dps 50 and two per accepted step at dps 110) are gone.
  * stage checkpoint + resume: after stage [1] (kernel transport) and stage
    [2] (density transport) the stage's product -- every stored Taylor segment
    as exact decimal strings at the working precision, the threshold moments,
    the statistics and the check results -- is written to a directory with a
    MANIFEST (sha256 per file, run parameters, producer); run(..., resume=DIR)
    verifies every pin and parameter before computing anything, rebuilds the
    stage objects from the stored series and re-enters at the first missing
    stage.  The acceptance criterion of a step (trailing-window tail bound
    below tol_step, fail-closed) is untouched.

Its host module is eval_row33.py beside it (the exact-rational systems +
RatDE); the served entry point with pinned files and exit codes is
../lbl3vp-eps0-evaluate.py.
"""
import os, sys, math, time, bisect
import mpmath as mp
import sympy as sp

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
sys.path.insert(0, os.path.join(HERE, '..'))
import k33lib
import lbl3disp as LD
from detransport.quad import quad_refine

ORD = 3  # jet truncation order: coefficients of e^0, e^1, e^2


# ---------------- order-3 jet arithmetic (Taylor in eps at 0) ----------------
class J:
    __slots__ = ('c',)

    def __init__(self, c):
        self.c = list(c) + [mp.mpf(0)] * (ORD - len(c))

    @staticmethod
    def const(x):
        return J([mp.mpf(x) if not isinstance(x, (mp.mpf, mp.mpc)) else x])

    def __add__(s, o):
        o = o if isinstance(o, J) else J.const(o)
        return J([a + b for a, b in zip(s.c, o.c)])

    __radd__ = __add__

    def __neg__(s):
        return J([-a for a in s.c])

    def __sub__(s, o):
        return s + (-(o if isinstance(o, J) else J.const(o)))

    def __rsub__(s, o):
        return (-s) + o

    def __mul__(s, o):
        if not isinstance(o, J):
            return J([a * o for a in s.c])
        c = [mp.mpf(0)] * ORD
        for i, a in enumerate(s.c):
            for k in range(ORD - i):
                c[i + k] += a * o.c[k]
        return J(c)

    __rmul__ = __mul__

    def inv(s):
        a0 = s.c[0]
        r = [1 / a0]
        for m in range(1, ORD):
            acc = mp.mpf(0)
            for k in range(1, m + 1):
                acc += s.c[k] * r[m - k]
            r.append(-acc / a0)
        return J(r)

    def __truediv__(s, o):
        return s * (o.inv() if isinstance(o, J) else J.const(o).inv())

    def __rtruediv__(s, o):
        return (o if isinstance(o, J) else J.const(o)) * s.inv()

    def log(s):
        # log(a0) + log(1 + t), t = (s - a0)/a0
        a0 = s.c[0]
        t = [x / a0 for x in s.c[1:]]  # coeffs of e^1..
        out = [mp.log(a0), mp.mpf(0), mp.mpf(0)]
        if ORD >= 2:
            out[1] = t[0]
        if ORD >= 3:
            out[2] = t[1] - t[0] * t[0] / 2
        return J(out)

    def exp(s):
        e0 = mp.exp(s.c[0])
        out = [e0, mp.mpf(0), mp.mpf(0)]
        if ORD >= 2:
            out[1] = e0 * s.c[1]
        if ORD >= 3:
            out[2] = e0 * (s.c[2] + s.c[1] * s.c[1] / 2)
        return J(out)

    def pow_c(s, a):          # s^a, a constant
        return (s.log() * a).exp()


def jpow(x, al):
    """x^al for positive scalar x and jet exponent al: exp(al*log x)."""
    return (al * mp.log(x)).exp() if isinstance(al, J) else mp.mpf(x) ** al


EPSJ = None  # set at working precision: J([0,1,0])


def gamma_jet(a):
    """Jet of Gamma(a + e) to order 2 (a > 0 scalar)."""
    g = mp.gamma(a)
    p0 = mp.psi(0, a)
    p1 = mp.psi(1, a)
    return J([g, g * p0, g * (p0 * p0 + p1) / 2])


def ghat_jet():
    """ghat(e) = Gamma(1+e)Gamma(1-e)/Gamma(1-2e) = 1 - 2g e + (2g^2-z2)e^2."""
    g = mp.euler
    return J([mp.mpf(1), -2 * g, 2 * g * g - mp.pi ** 2 / 6])


def gamma_at(argj):
    """Gamma(argj) for jet argument with argj.c0 > 0."""
    a0, a1, a2 = argj.c
    g = mp.gamma(a0)
    p0 = mp.psi(0, a0)
    p1 = mp.psi(1, a0)
    # Gamma(a0 + d), d = a1 e + a2 e^2: G*(1 + p0 d + (p0^2+p1)d^2/2)
    c0 = g
    c1 = g * p0 * a1
    c2 = g * (p0 * a2 + (p0 * p0 + p1) * a1 * a1 / 2)
    return J([c0, c1, c2])


# ---------------- parametric powdiff jets (kernel 2-folds) ----------------
def pd_jets(A, B, X, p0, korder=3):
    """[d^k/dp^k of (A^{-p} - C^{-p})/B]_{p=p0} / k!  for k=0..korder-1,
    C = A + B*X; cancellation-free (log1p forms).  A, C > 0."""
    C = A + B * X
    z = B * X / A
    pd0 = k33lib._powdiff(A, B, X, p0)
    out = [pd0]
    if korder == 1:
        return out
    lA = mp.log(A)
    l1 = mp.log1p(z)                      # ln(C/A)
    lC = lA + l1
    Amp = A ** (-p0)
    # l1/B stable: X*(log1p(z)/z)/A
    if z == 0:
        l1oB = X / A
    else:
        l1oB = X * (l1 / z) / A
    pd1 = Amp * l1oB - pd0 * lC
    out.append(pd1)
    if korder >= 3:
        pd2 = (-Amp * l1oB * (lA + lC) + pd0 * lC * lC) / 2
        out.append(pd2)
    return out


def kn_jets(nmax, n, work, k2_order=3):
    """Jets (order 2 in e, order-3 for K2) of threshold moments K_j(e),
    j=0..nmax-1, k33lib.kn_moments structure with pd_jets. Returns list of J."""
    c = mp.mpf(4) / 3
    nodes, weights = k33lib.gl_nodes(n, work)
    half = mp.mpf('0.5')
    tots = [[mp.mpf(0)] * ORD for _ in range(nmax)]
    for t1, w1 in zip(nodes, weights):
        x1 = half * (t1 + 1)
        L = 1 - x1
        accs = [[mp.mpf(0)] * ORD for _ in range(nmax)]
        for t2, w2 in zip(nodes, weights):
            x2 = L * half * (t2 + 1)
            X = 1 - x1 - x2
            A1 = 1 - c * x2 * X
            B = x1 + c * x2
            C = A1 + B * X
            lA = mp.log(A1)
            z = B * X / A1
            l1 = mp.log1p(z)
            lC = lA + l1
            l1oB = X / A1 if z == 0 else X * (l1 / z) / A1
            Ainv = 1 / A1
            Cinv = 1 / C
            Ampow = Ainv                  # A1^{-p} at p = 1+j
            Cpow = Cinv                   # C^{-p}
            pd0 = k33lib._powdiff(A1, B, X, mp.mpf(1))
            lAlC = lA + lC
            lC2 = lC * lC
            x1p = mp.mpf(1)
            for j in range(nmax):
                # stable recursion pd0(p+1) = (pd0(p) + X*C^{-(p+1)})/A1
                pd1 = Ampow * l1oB - pd0 * lC
                pd2 = (-Ampow * l1oB * lAlC + pd0 * lC2) / 2
                accs[j][0] += w2 * x1p * pd0
                accs[j][1] += w2 * x1p * pd1
                accs[j][2] += w2 * x1p * pd2
                x1p *= x1
                Cpow *= Cinv
                pd0 = (pd0 + X * Cpow) * Ainv
                Ampow *= Ainv
        for j in range(nmax):
            for k in range(ORD):
                tots[j][k] += w1 * accs[j][k] * (L * half)
    # K_j(e) = Gamma(2+e) * (2+e)_j/j! * (-1)^j * tots_j(e)*half/(1+e+j)
    # pd was taken at p = 1+j; true p = 1+e+j -> tots_j(e) = jet in e via pd_jets
    # NOTE: pd_jets above are d/dp jets == d/de jets since p = (1+j) + e.
    G2 = gamma_jet(mp.mpf(2))
    out = []
    poch = J.const(1)
    e = J([mp.mpf(0), mp.mpf(1)])
    for j in range(nmax):
        if j > 0:
            poch = poch * (J.const(1 + j) + e) / j
        pref = G2 * poch * ((-1) ** j) * half
        tj = J(tots[j])
        out.append(pref * tj / (J.const(1 + j) + e))
    return out


def kn_jets_refine(nmax, dps, guard, n0, label, verbose=True):
    """n-doubling refine loop for kn_jets: accept when successive-n agreement
    over all returned coefficients beats 10^-(dps+guard). Fail-closed."""
    work = dps + guard + 20
    tol = mp.mpf(10) ** (-(dps + guard))
    n = n0
    prev = None
    for it in range(6):
        cur = kn_jets(nmax, n, work)
        if prev is not None:
            agr = max(max(abs(a - b) for a, b in zip(x.c, y.c))
                      for x, y in zip(cur, prev))
            if agr <= tol:
                if verbose:
                    print(f"    [certified] {label}: GL n={n} vs n={npr} agreement "
                          f"{mp.nstr(agr, 3)} <= tol {mp.nstr(tol, 3)}", flush=True)
                return cur
        prev, npr = cur, n
        n = int(n * 1.5) + 8
    raise RuntimeError(f"{label}: kn_jets not converged at GL n={npr} "
                       f"(agreement {mp.nstr(agr, 3)} > tol {mp.nstr(tol, 3)}) -- fail-closed")


# ---------------- box-family master Laurent jets at w'=5 (seeds) -------------
def _lnF_ints(f, n, work, kmax=2):
    """int_0^1 ln^k F(x) dx for k=0..kmax via GL (F analytic > 0)."""
    nodes, weights = k33lib.gl_nodes(n, work)
    half = mp.mpf('0.5')
    tots = [mp.mpf(0)] * (kmax + 1)
    for t, w in zip(nodes, weights):
        x = half * (t + 1)
        lf = mp.log(f(x))
        p = mp.mpf(1)
        for k in range(kmax + 1):
            tots[k] += w * p
            p *= lf
    return [half * t for t in tots]


def master_laurent_jets(wp, n, work):
    """8 box1eq masters as Laurent jets: list of (k0, J) meaning
    sum_i J.c[i] * e^(k0+i). Master order = k33lib.masters_at."""
    wp = mp.mpf(wp)
    c = mp.mpf(4) / 3
    e = J([mp.mpf(0), mp.mpf(1)])
    G1p = gamma_jet(mp.mpf(1))            # Gamma(1+e)
    # tadpole(m2) = Gamma(e)/(1-e) m2^{1-e} = (1/e)*Gamma(1+e)/(1-e)*m2^{1-e}
    def tad(m2):
        return (-1, G1p / (J.const(1) - e) * jpow(m2, J.const(1) - e))

    # bubbles: Gamma(e) * int F^{-e} = (1/e)*Gamma(1+e)*[1 - e I1 + e^2 I2/2]
    def bub(F):
        I0, I1, I2 = _lnF_ints(F, n, work, 2)
        return (-1, G1p * J([I0, -I1, I2 / 2]))

    Fs = lambda x: 1 + x * (wp - 1) + x * (1 - x)
    Fu = lambda x: 1 - c * x * (1 - x)

    # triangles: -Gamma(1+e) * int pd(A,Cc,X,e)/e ; pd(p0=0)=0 ->
    # pd(...,e)/e = pd1 + e*pd2 (+O(e^2))
    def tri(Afun, Cfun, Xfun):
        nodes, weights = k33lib.gl_nodes(n, work)
        half = mp.mpf('0.5')
        t1 = mp.mpf(0)
        t2 = mp.mpf(0)
        for t, w in zip(nodes, weights):
            x = half * (t + 1)
            A, B, X = Afun(x), Cfun(x), Xfun(x)
            pds = pd_jets(A, B, X, mp.mpf(0), 3)
            t1 += w * pds[1]
            t2 += w * pds[2]
        t1 *= half
        t2 *= half
        return (0, -G1p * J([t1, t2, mp.mpf(0)]))

    tri_s = tri(lambda x: 1 + x * (wp - 1), lambda x: x, lambda x: 1 - x)
    tri_u = tri(lambda x: 1 - c * x * (1 - x), lambda x: wp - 1 + c * x,
                lambda x: 1 - x)
    tri_c = tri(lambda x: mp.mpf(1), lambda x: -c * x, lambda x: 1 - x)

    # box: Gamma(2+e)/(1+e) * II pd(A,B,X,1+e)
    nodes, weights = k33lib.gl_nodes(n, work)
    half = mp.mpf('0.5')
    tb = [mp.mpf(0)] * ORD
    for t1_, w1 in zip(nodes, weights):
        x1 = half * (t1_ + 1)
        L = 1 - x1
        acc = [mp.mpf(0)] * ORD
        for t2_, w2 in zip(nodes, weights):
            x2 = L * half * (t2_ + 1)
            X = 1 - x1 - x2
            A = 1 + x1 * (wp - 1) - c * x2 * X
            B = x1 + c * x2
            pds = pd_jets(A, B, X, mp.mpf(1), 3)
            for k in range(ORD):
                acc[k] += w2 * pds[k]
        for k in range(ORD):
            tb[k] += w1 * acc[k] * (L * half)
    boxj = gamma_jet(mp.mpf(2)) / (J.const(1) + e) * J(tb) * half
    return [(0, boxj), tad(wp), tad(1), bub(Fs), bub(Fu), tri_s, tri_u, tri_c]


def laurent_to_layers(lj, kmin=-1, nlay=3):
    """(k0, J) -> [layer_kmin, ..., layer_{kmin+nlay-1}] scalar list."""
    k0, jet = lj
    out = []
    for k in range(kmin, kmin + nlay):
        i = k - k0
        out.append(jet.c[i] if 0 <= i < ORD else mp.mpf(0))
    return out


# ---------------- block jet DE (layers of the eps-Taylor of A) ---------------
def de_eps_taylor(de_json, korder):
    """A(w, 4-2e) = A0 + A1 e + ... : list of korder sympy matrices (nested
    lists of exact rational-in-w expressions)."""
    W, D, E = sp.symbols('w d epsx')
    nm = len(de_json['masters'])
    out = [[[sp.Integer(0)] * nm for _ in range(nm)] for _ in range(korder)]
    for i in range(nm):
        for j in range(nm):
            ex = sp.sympify(de_json['A'][i][j], locals={'w': W, 'd': D})
            if ex == 0:
                continue
            ex = sp.cancel(ex.subs(D, 4 - 2 * E))
            ser = sp.series(ex, E, 0, korder).removeO()
            for k in range(korder):
                ck = sp.cancel(ser.coeff(E, k))
                out[k][i][j] = ck
    return out


def block_de_json(de_json, nlay):
    """Lower-triangular block DE for Laurent layers M_l' = sum_k A_k M_{l-k}.
    Returns a de_json-like dict consumable by eval_row33.RatDE at eps=0."""
    Ak = de_eps_taylor(de_json, nlay)
    nm = len(de_json['masters'])
    N = nm * nlay
    A = [['0'] * N for _ in range(N)]
    masters = []
    for l in range(nlay):
        for m in de_json['masters']:
            masters.append(tuple(m) + (l,))
    for li in range(nlay):
        for lj in range(li + 1):
            k = li - lj
            for i in range(nm):
                for j in range(nm):
                    ex = Ak[k][i][j]
                    if ex != 0:
                        A[li * nm + i][lj * nm + j] = str(ex)
    return {'masters': masters, 'A': A}


# ---------------- stage checkpoint: exact number strings + the stage store ----
def _num_str(x):
    """Exact decimal form of an mpf (string) or mpc ([re, im]) at the working
    precision: mpmath's repr digits, which parse back to the same value."""
    from mpmath.libmp import to_str, repr_dps
    nd = repr_dps(mp.mp.prec)
    if isinstance(x, mp.mpc):
        return [to_str(x.real._mpf_, nd), to_str(x.imag._mpf_, nd)]
    return to_str(mp.mpf(x)._mpf_, nd)


def _num_val(s):
    """Inverse of _num_str at the same working precision."""
    if isinstance(s, list):
        return mp.mpc(mp.mpf(s[0]), mp.mpf(s[1]))
    return mp.mpf(s)


def _sha256_file(path):
    import hashlib
    h = hashlib.sha256()
    with open(path, 'rb') as f:
        for chunk in iter(lambda: f.read(1 << 20), b''):
            h.update(chunk)
    return h.hexdigest()


def _plain(d):
    """Check-result dict with mp numbers as short strings (for the MANIFEST)."""
    return {k: (mp.nstr(v, 20) if isinstance(v, (mp.mpf, mp.mpc)) else v) for k, v in d.items()}


class CheckpointRefused(Exception):
    """A resume directory that cannot be used, named before anything computes:
    exit_code 4 when it or a file it names is absent, 3 when a pin, the engine
    identity or a run parameter does not match."""

    def __init__(self, msg, exit_code):
        super().__init__(msg)
        self.exit_code = exit_code


class StageStore:
    """The on-disk stage checkpoint of one run: DIR/MANIFEST.json (format,
    producer, run parameters, one entry per saved stage with the file's sha256)
    plus one JSON file per stage.  Writes are atomic (temp file + rename), and
    a stage file is named by the MANIFEST only after its series were re-read
    and compared with the in-memory values.  open() verifies every pin, the
    engine identity and every parameter before anything is computed; each
    mismatch is a CheckpointRefused naming the file or parameter."""

    FORMAT = 'row33-eps0-stage-checkpoint-1'
    MANIFEST = 'MANIFEST.json'
    FILES = {1: 'stage1_kernel.json', 2: 'stage2_rho.json'}

    def __init__(self, path, params, producer):
        self.dir = path
        self.params = params
        self.manifest = {'format': self.FORMAT, 'producer': producer, 'params': params, 'stages': {}}

    @staticmethod
    def producer(VP):
        """Engine identity + launch line: the sha256 of the engine, its host module and the data file."""
        return {'engine': os.path.basename(__file__), 'engine_sha256': _sha256_file(os.path.abspath(__file__)),
                'host_sha256': _sha256_file(os.path.abspath(VP.__file__)),
                'data_sha256': _sha256_file(os.path.join(os.path.dirname(os.path.abspath(VP.__file__)), 'row33_data.json')),
                'argv': list(sys.argv), 'started_utc': time.strftime('%Y-%m-%dT%H:%M:%SZ', time.gmtime())}

    def has(self, stage):
        # stages count as an unbroken prefix 1, 2, ...: a stage 2 without stage 1 is not used
        return all(str(s) in self.manifest['stages'] for s in range(1, stage + 1))

    def load(self, stage):
        import json
        ent = self.manifest['stages'][str(stage)]
        with open(os.path.join(self.dir, ent['file']), encoding='utf-8') as f:
            return json.load(f)

    def save(self, stage, payload, wall_s, verbose=True):
        import json
        os.makedirs(self.dir, exist_ok=True)
        if not self.manifest['stages']:          # first save of a fresh run: an earlier run's files go
            for name in [self.MANIFEST] + list(self.FILES.values()):
                p = os.path.join(self.dir, name)
                if os.path.exists(p):
                    os.remove(p)
        name = self.FILES[stage]
        path = os.path.join(self.dir, name)
        t0 = time.time()
        with open(path + '.tmp', 'w', encoding='utf-8') as f:
            json.dump(payload, f, separators=(',', ':'))
        os.replace(path + '.tmp', path)
        sha = _sha256_file(path)
        self.manifest['stages'][str(stage)] = {
            'file': name, 'sha256': sha, 'bytes': os.path.getsize(path),
            'saved_utc': time.strftime('%Y-%m-%dT%H:%M:%SZ', time.gmtime()), 'stage_wall_s': round(wall_s, 1),
            'stats': payload.get('stats'), 'checks': payload.get('checks')}
        mpath = os.path.join(self.dir, self.MANIFEST)
        with open(mpath + '.tmp', 'w', encoding='utf-8') as f:
            json.dump(self.manifest, f, indent=1)
        os.replace(mpath + '.tmp', mpath)
        if verbose:
            print(f"    [checkpoint] stage [{stage}] saved in {time.time()-t0:.1f}s: {path} "
                  f"({os.path.getsize(path)} bytes, sha256 {sha[:16]}...); MANIFEST updated; "
                  f"resume with --resume {self.dir}", flush=True)

    @classmethod
    def open(cls, path, params, producer):
        """Open DIR for resume; every refusal names its cause and nothing has been computed yet."""
        import json
        if not os.path.isdir(path):
            raise CheckpointRefused(f"REFUSED: --resume {path}: no such directory; nothing computed", 4)
        mpath = os.path.join(path, cls.MANIFEST)
        if not os.path.exists(mpath):
            raise CheckpointRefused(f"REFUSED: --resume {path}: {cls.MANIFEST} is missing (not a stage "
                                    f"checkpoint directory); nothing computed", 4)
        try:
            with open(mpath, encoding='utf-8') as f:
                man = json.load(f)
        except ValueError as e:
            raise CheckpointRefused(f"REFUSED: --resume {path}: {cls.MANIFEST} is not readable JSON ({e}); "
                                    f"nothing computed", 3)
        if man.get('format') != cls.FORMAT:
            raise CheckpointRefused(f"REFUSED: --resume {path}: format {man.get('format')!r} is not "
                                    f"{cls.FORMAT!r}; nothing computed", 3)
        prod = man.get('producer', {})
        for k in ('engine_sha256', 'host_sha256', 'data_sha256'):
            if prod.get(k) != producer[k]:
                raise CheckpointRefused(f"REFUSED: --resume {path}: {k} {prod.get(k)} in the MANIFEST is not "
                                        f"this engine's {producer[k]} (written by other code or data); "
                                        f"nothing computed", 3)
        mp_ = man.get('params', {})
        bad = [k for k in params if mp_.get(k) != params[k]]
        if bad:
            detail = '; '.join(f"{k}: checkpoint {mp_.get(k)!r}, requested {params[k]!r}" for k in bad)
            raise CheckpointRefused(f"REFUSED: --resume {path}: run parameters differ ({detail}); "
                                    f"nothing computed", 3)
        for st, ent in sorted(man.get('stages', {}).items()):
            p = os.path.join(path, ent['file'])
            if not os.path.exists(p):
                raise CheckpointRefused(f"REFUSED: --resume {path}: stage [{st}] file {ent['file']} named by "
                                        f"the MANIFEST is missing; nothing computed", 4)
            sha = _sha256_file(p)
            if sha != ent['sha256']:
                raise CheckpointRefused(f"REFUSED: --resume {path}: {ent['file']} sha256 {sha} does not match "
                                        f"the MANIFEST pin {ent['sha256']} -- the stage [{st}] file was "
                                        f"altered; nothing computed", 3)
        self = cls(path, params, prod)
        self.manifest = man
        return self


class JetMarch:
    """Stored-series transport: march a RatDE from w0 outward, storing the
    local Taylor series of selected components; eval anywhere in the covered
    range by Horner on the containing step's series. Per-step acceptance =
    trailing-window tail bound < tol_step, fail-closed.

    Step controller: RatDE.step returns the Taylor coefficients C of the local
    solution and their sum at the proposed h.  C depends on the start point,
    the state and the order N only -- not on h -- so a rejected step is
    re-summed at h/2 from the same C (bit for bit what RatDE.step would return
    for h/2) instead of being recomputed; only an N-raise recomputes.  The
    accepted steps and their bounds are therefore identical to the earlier
    recompute-on-halving form; stats['evals'] counts the RatDE.step calls.

    Checkpoint form: to_record() / from_record() turn a built march into (and
    back from) exact decimal strings at the working precision -- the whole
    product of a transport stage; a restored march evaluates without its
    differential system."""

    def __init__(self, de, M0, w0, N, sf, comps, dps, guard, nsteps_guess=400):
        self.de = de
        self.N = N
        self.sf = mp.mpf(sf)
        self.comps = comps            # component indices whose series we keep
        self.tol_step = mp.mpf(10) ** (-(dps + guard)) / (4 * nsteps_guess)
        self.segsL = []               # descending marches: list (wstart, h, {c: ser})
        self.segsR = []
        self.M0 = list(M0)
        self.w0 = mp.mpf(w0)
        self.bound_acc = mp.mpf(0)    # accumulated per-step tail bounds
        self.stats = {'steps': 0, 'halved': 0, 'nraised': 0, 'evals': 0}

    def _tail_bound(self, C, h):
        r = self.sf / (1 - self.sf)
        win = mp.mpf(0)
        ah = abs(h)
        for m in range(max(0, len(C) - 8), len(C)):
            t = max(abs(x) for x in C[m]) * ah ** m
            if t > win:
                win = t
        return win * r / (1 - r) if r < 1 else mp.inf

    def _resum(self, C, h):
        """Sum the local Taylor series C at step h in RatDE.step's own loop
        order, so the result is bit for bit the value RatDE.step returns for h."""
        n = len(C[0])
        v = [mp.mpc(0)] * n
        hp = mp.mpc(1)
        for m in range(len(C)):
            cm = C[m]
            for i in range(n):
                v[i] += cm[i] * hp
            hp *= h
        return v

    def _march(self, w_end, store):
        M = list(self.M0)
        w = self.w0
        sgn = 1 if w_end > w else -1
        N = self.N
        while (w_end - w) * sgn > abs(w_end) * mp.mpf('1e-45') + mp.mpf('1e-45'):
            d = self.de.nearest(w)
            h = sgn * min(self.sf * d, abs(w_end - w))
            tries = 0
            Mn, C = self.de.step(M, w, h, N)
            self.stats['evals'] += 1
            scale = max(mp.mpf(1), max(abs(x) for x in M))
            while True:
                bnd = self._tail_bound(C, h)
                if bnd <= self.tol_step * scale:
                    break
                tries += 1
                if tries <= 4:
                    h = h / 2
                    self.stats['halved'] += 1
                    Mn = self._resum(C, h)          # same series, halved step: no recomputation
                elif tries <= 6:
                    N = int(N * 3 / 2)
                    self.stats['nraised'] += 1
                    Mn, C = self.de.step(M, w, h, N)  # a higher order needs the series recomputed
                    self.stats['evals'] += 1
                else:
                    raise RuntimeError(
                        f"JetMarch: step at w={mp.nstr(w,8)} h={mp.nstr(h,4)} "
                        f"tail bound {mp.nstr(bnd,3)} > tol {mp.nstr(self.tol_step,3)} "
                        f"after 4 halvings + 2 N-raises (N={N}) -- fail-closed")
            self.bound_acc += bnd
            store.append((w, h, {c: [C[m][c] for m in range(len(C))] for c in self.comps}))
            M = Mn
            w = w + h
            self.stats['steps'] += 1
        return M

    def build(self, w_lo, w_hi):
        if mp.mpf(w_lo) < self.w0:
            self._march(mp.mpf(w_lo), self.segsL)
        if mp.mpf(w_hi) > self.w0:
            self._march(mp.mpf(w_hi), self.segsR)
        self._index(w_lo, w_hi)

    def _index(self, w_lo, w_hi):
        # unified interval table: (lo_float, w0, sers), sorted by lo
        iv = []
        for (w0, h, sers) in self.segsL + self.segsR:
            lo = w0 + h if h < 0 else w0
            iv.append((float(lo), w0, sers))
        iv.sort(key=lambda t: t[0])
        self.iv = iv
        self.ivlo = [t[0] for t in iv]
        self.lo = mp.mpf(w_lo)
        self.hi = mp.mpf(w_hi)

    def to_record(self):
        """Checkpoint form of a built march: parameters, statistics and every
        stored segment (w, h, per-component series) as exact decimal strings."""
        def segs(lst):
            return [[_num_str(w), _num_str(h),
                     {str(c): [_num_str(x) for x in ser] for c, ser in sers.items()}]
                    for (w, h, sers) in lst]
        return {'w0': _num_str(self.w0), 'N': self.N, 'sf': _num_str(self.sf),
                'comps': [int(c) for c in self.comps], 'tol_step': _num_str(self.tol_step),
                'bound_acc': _num_str(self.bound_acc), 'M0': [_num_str(x) for x in self.M0],
                'lo': _num_str(self.lo), 'hi': _num_str(self.hi), 'stats': dict(self.stats),
                'segsL': segs(self.segsL), 'segsR': segs(self.segsR)}

    @classmethod
    def from_record(cls, rec):
        """Rebuild a march from its checkpoint form.  No differential system is
        needed (the stored series are the whole product); the interval table is
        rebuilt exactly as build() does."""
        self = cls.__new__(cls)
        self.de = None
        self.N = int(rec['N'])
        self.sf = _num_val(rec['sf'])
        self.comps = [int(c) for c in rec['comps']]
        self.tol_step = _num_val(rec['tol_step'])
        self.M0 = [_num_val(x) for x in rec['M0']]
        self.w0 = _num_val(rec['w0'])
        self.bound_acc = _num_val(rec['bound_acc'])
        self.stats = dict(rec['stats'])

        def segs(lst):
            return [(_num_val(w), _num_val(h),
                     {int(c): [_num_val(x) for x in ser] for c, ser in sers.items()})
                    for (w, h, sers) in lst]
        self.segsL = segs(rec['segsL'])
        self.segsR = segs(rec['segsR'])
        self._index(_num_val(rec['lo']), _num_val(rec['hi']))
        return self

    def record_matches(self, rec):
        """True iff every number in the checkpoint form rec parses back to the
        in-memory value with the same type -- the round-trip check made before
        a stage file is named by the MANIFEST."""
        def same(x, s):
            y = _num_val(s)
            return type(y) is type(x) and y == x
        if not (same(self.w0, rec['w0']) and same(self.sf, rec['sf']) and same(self.tol_step, rec['tol_step'])
                and same(self.bound_acc, rec['bound_acc']) and same(self.lo, rec['lo']) and same(self.hi, rec['hi'])
                and len(rec['M0']) == len(self.M0) and all(same(x, s) for x, s in zip(self.M0, rec['M0']))):
            return False
        for mem, stored in ((self.segsL, rec['segsL']), (self.segsR, rec['segsR'])):
            if len(mem) != len(stored):
                return False
            for (w, h, sers), (ws, hs, serss) in zip(mem, stored):
                if not (same(w, ws) and same(h, hs)) or set(serss) != {str(c) for c in sers}:
                    return False
                for c, ser in sers.items():
                    ss = serss[str(c)]
                    if len(ss) != len(ser) or not all(same(x, s) for x, s in zip(ser, ss)):
                        return False
        return True

    def eval(self, w, comp):
        w = mp.mpf(w)
        if not (self.lo <= w <= self.hi):
            raise RuntimeError(f"JetMarch.eval: w={mp.nstr(w,8)} outside covered "
                               f"[{mp.nstr(self.lo,8)}, {mp.nstr(self.hi,8)}]")
        i = bisect.bisect_right(self.ivlo, float(w)) - 1
        i = min(max(i, 0), len(self.iv) - 1)
        _, w0, sers = self.iv[i]
        x = w - w0
        acc = mp.mpf(0)
        for cm in reversed(sers[comp]):
            acc = acc * x + cm
        return mp.re(acc)


# ---------------- rho jets (2-layer: e^-1, e^0) -------------------------------
class RhoJets:
    """rho(w)/pi Laurent layers r_-1(w), r_0(w).
    (1,9): closed 2-body jets; [9,wsw]: - c*psi0/pi Frobenius-jet correction;
    (wsw, WMAX]: 2-layer sigvp block jet transport."""

    def __init__(self, VP, dps, guard, K=None, wsw=10):
        self.VP = VP
        self.dps = dps
        self.guard = guard
        self.wsw = mp.mpf(wsw)
        self.K = K if K is not None else int(mp.mpf('1.45') * (dps + guard)) + 8
        e = J([mp.mpf(0), mp.mpf(1)])
        self.e = e
        self.Ghat = ghat_jet() / (J.const(1) - 2 * e)     # G(e)*e
        self.G1m = gamma_at(J([mp.mpf(1), mp.mpf(-1)]))   # Gamma(1-e)
        self.G2m2 = gamma_at(J([mp.mpf(2), mp.mpf(-2)]))  # Gamma(2-2e)
        self.N2 = self.G1m / self.G2m2
        self.Tj = (-1, gamma_jet(mp.mpf(1)) / (J.const(1) - e))  # tadpole(1) Laurent
        self._Ak = de_eps_taylor(VP.DATA['sigvp_de'], 2)  # A0, A1 exact sympy
        self._frob()
        self._c_jet()
        self.march = None

    # ---- closed 2-body jets: rho/pi = -(Ghat/e) P u^{-2e} w^{e-1} ----
    def rho2b_layers(self, w):
        w = mp.mpf(w)
        u = w - 1
        e = self.e
        P = -(J.const(1) - 2 * e) * (w + 1) / u + e * (u / 6)
        anal = -(self.Ghat * P * jpow(u, -2 * e) * jpow(w, e - 1))
        return anal.c[0], anal.c[1]        # r_-1, r_0

    def drho2b_layers(self, w):
        """d/dw of the 2-body rho/pi layers (jets of drho_exact/pi... via
        analytic derivative of the closed form)."""
        w = mp.mpf(w)
        u = w - 1
        e = self.e
        P = -(J.const(1) - 2 * e) * (w + 1) / u + e * (u / 6)
        Pp = 2 * (J.const(1) - 2 * e) / (u * u) + e * (mp.mpf(1) / 6)
        anal = -(self.Ghat * (Pp + P * ((-2) * e / u + (e - J.const(1)) / w))
                 * jpow(u, -2 * e) * jpow(w, e - 1))
        return anal.c[0], anal.c[1]

    def ImTB_layers(self, w):
        w = mp.mpf(w)
        u = w - 1
        e = self.e
        ImB = mp.pi * self.N2 * jpow(u, J.const(1) - 2 * e) * jpow(w, e - J.const(1))
        k0, tj = self.Tj
        lau = tj * ImB                     # Laurent from e^-1
        return lau.c[0], lau.c[1]          # layers -1, 0

    def _Aat_jet(self, i, j, w):
        """Jet (2 coeffs used) of A[i][j](w, 4-2e) at scalar w."""
        W = sp.symbols('w')
        out = []
        for k in range(2):
            ex = self._Ak[k][i][j]
            if ex == 0:
                out.append(mp.mpf(0))
                continue
            num, den = sp.fraction(sp.together(ex))
            pn = sp.Poly(num, W).all_coeffs() if num.has(W) else [num]
            pdn = sp.Poly(den, W).all_coeffs() if den.has(W) else [den]
            nv = mp.mpf(0)
            for cc in pn:
                nv = nv * w + mp.mpf(sp.Rational(cc).p) / mp.mpf(sp.Rational(cc).q)
            dv = mp.mpf(0)
            for cc in pdn:
                dv = dv * w + mp.mpf(sp.Rational(cc).p) / mp.mpf(sp.Rational(cc).q)
            out.append(nv / dv)
        return J(out)

    # ---- Frobenius psi jets at w=9 (the Frobenius construction of the density, in jets) ----
    def _taylor_at9_jet(self, i, j, K):
        Rk = []
        Sk = []
        for k in range(2):
            ex = self._Ak[k][i][j]
            R, S = _taylor_at9_expr(ex, K)
            Rk.append(R)
            Sk.append(S)
        R = J([Rk[0], Rk[1]])
        S = [J([Sk[0][m], Sk[1][m]]) for m in range(len(Sk[0]))]
        return R, S

    def _frob(self):
        K = self.K
        R = [[None] * 2 for _ in range(2)]
        S = [[None] * 2 for _ in range(2)]
        for i in range(2):
            for j in range(2):
                R[i][j], S[i][j] = self._taylor_at9_jet(i, j, K)
        alpha = J([mp.mpf(2), mp.mpf(-2)])
        v0 = [R[0][1], alpha - R[0][0]]
        ps = [v0]
        for k in range(1, K + 1):
            rhs = [J.const(0), J.const(0)]
            for l in range(k):
                for i in range(2):
                    rhs[i] = rhs[i] + S[i][0][k - 1 - l] * ps[l][0] \
                        + S[i][1][k - 1 - l] * ps[l][1]
            a = alpha + k - R[0][0]
            b = J.const(0) - R[0][1]
            cc = J.const(0) - R[1][0]
            dd = alpha + k - R[1][1]
            det = a * dd - b * cc
            ps.append([(dd * rhs[0] - b * rhs[1]) / det,
                       (a * rhs[1] - cc * rhs[0]) / det])
        self.alpha = alpha
        self.v0 = v0
        self.ps = ps

    def psi_jets(self, w):
        """(psi0, psi1) as jets at scalar w (|w-9| <= ~1.5), with tail check."""
        uu = mp.mpf(w) - 9
        a0 = J.const(0)
        a1 = J.const(0)
        up = mp.mpf(1)
        lastmax = mp.mpf(0)
        for k in range(self.K + 1):
            a0 = a0 + self.ps[k][0] * up
            a1 = a1 + self.ps[k][1] * up
            t = max(abs(x) for x in self.ps[k][0].c[:2]) * abs(up)
            if k > self.K - 8:
                lastmax = max(lastmax, t)
            up *= uu
        tol = mp.mpf(10) ** (-(self.dps + self.guard))
        r = abs(uu) / 8
        tail = lastmax * r / (1 - r)
        sc = max(mp.mpf(1), max(abs(x) for x in a0.c[:2]))
        if tail > tol * sc * 100:
            raise RuntimeError(
                f"psi_jets: Frobenius tail bound {mp.nstr(tail,3)} > tol at w={w}; "
                f"raise K (K={self.K}) -- fail-closed")
        f = jpow(uu, self.alpha)
        return f * a0, f * a1

    def _c_jet(self):
        e = self.e
        d = J([mp.mpf(4), mp.mpf(-2)])
        halfd = (d - J.const(1)) / 2
        Bt = gamma_at(halfd) * gamma_at(halfd) / gamma_at(d - J.const(1))
        A0c = -(mp.pi / 16) * self.N2 * self.N2 * jpow(9, -(d - J.const(2)) / 2) \
            * jpow(48, (d - J.const(3)) / 2) * jpow(4, -(d - J.const(2)) / 2) \
            * Bt * jpow(mp.mpf(2) / 3, d - J.const(2))
        self.cj = A0c / self.v0[0]

    def rho_layers(self, w):
        """(r_-1, r_0) of rho/pi at scalar w in (1, wsw]."""
        w = mp.mpf(w)
        rm1, r0 = self.rho2b_layers(w)
        if w >= 9:
            p0, _ = self.psi_jets(w)
            corr = self.cj * p0 * (1 / mp.pi)
            # c*psi0 is ANALYTIC: its e^-1 layer must vanish (2-body-only pole)
            rm1 = rm1 - 0
            r0 = r0 - corr.c[0]
        return rm1, r0

    def seed_jet_at(self, w0):
        """8-vector seed for the 2-layer sigvp block at w0 (layers -1,0 of
        (ImM0, ImM1, 0, ImTB))."""
        w0 = mp.mpf(w0)
        i0m1, i00 = self.rho2b_layers(w0)
        d0m1, d00 = self.drho2b_layers(w0)
        # Im0 = -pi*(rho/pi) ... rho33lib works with rho itself: Im0 = -rho
        Im0 = [-mp.pi * i0m1, -mp.pi * i00]
        dIm0 = [-mp.pi * d0m1, -mp.pi * d00]
        tbm1, tb0 = self.ImTB_layers(w0)
        A00 = self._Aat_jet(0, 0, w0)
        A01 = self._Aat_jet(0, 1, w0)
        A03 = self._Aat_jet(0, 3, w0)
        num = J([dIm0[0], dIm0[1]]) - A00 * J([Im0[0], Im0[1]]) \
            - A03 * J([tbm1, tb0])
        Im1 = num / A01
        p0, p1 = self.psi_jets(w0)
        cp0 = self.cj * p0
        cp1 = self.cj * p1
        # cp0/cp1 analytic -> contribute to layer 0 only (layer -1 = 0)
        Mm1 = [Im0[0], Im1.c[0], mp.mpf(0), tbm1]
        M00 = [Im0[1] + cp0.c[0], Im1.c[1] + cp1.c[0], mp.mpf(0), tb0]
        return Mm1 + M00

    def build_tail(self, WMAX, N, sf):
        bj = block_de_json(self.VP.DATA['sigvp_de'], 2)
        de = self.VP.RatDE(bj, sp.Integer(0), [0, 1, 9, -3])
        seed = self.seed_jet_at(self.wsw)
        self.march = JetMarch(de, seed, self.wsw, N, sf, comps=[0, 4],
                              dps=self.dps, guard=self.guard)
        self.march.build(self.wsw, WMAX)

    def attach_march(self, rec):
        """Constructor path for stage [2] from its checkpoint form: the tail
        transport rebuilt from the stored series in place of build_tail."""
        if int(rec.get('K', -1)) != self.K or _num_val(rec['wsw']) != self.wsw:
            raise CheckpointRefused(f"REFUSED: stage [2] checkpoint was written with K={rec.get('K')} "
                                    f"wsw={rec.get('wsw')}, this run has K={self.K} wsw={mp.nstr(self.wsw, 6)}", 3)
        self.march = JetMarch.from_record(rec['march'])

    def to_record(self):
        """Checkpoint form of stage [2]: the built tail transport."""
        return {'stage': 2, 'K': self.K, 'wsw': _num_str(self.wsw), 'march': self.march.to_record()}

    def rho_layers_tail(self, w):
        """(r_-1, r_0) for w > wsw from the transported Im-vector."""
        m_m1 = self.march.eval(w, 0)
        m_0 = self.march.eval(w, 4)
        return -m_m1 / mp.pi, -m_0 / mp.pi


def _taylor_at9_expr(ex, K):
    """Residue + analytic Taylor coeffs (length K+2, mpf) of a rational
    expression in w at w=9 (simple pole allowed)."""
    W = sp.symbols('w')
    if ex == 0:
        return mp.mpf(0), [mp.mpf(0)] * (K + 2)
    num, den = sp.fraction(sp.together(ex))
    nc = [mp.mpf(sp.Rational(c).p) / mp.mpf(sp.Rational(c).q)
          for c in (sp.Poly(num, W).all_coeffs() if num.has(W) else [num])]
    dc = [mp.mpf(sp.Rational(c).p) / mp.mpf(sp.Rational(c).q)
          for c in (sp.Poly(den, W).all_coeffs() if den.has(W) else [den])]
    M = K + 4

    def shift(coefs):
        p = [coefs[0]]
        for c in coefs[1:]:
            pn = [mp.mpf(0)] * (len(p) + 1)
            for m, a in enumerate(p):
                pn[m] += a * 9
                pn[m + 1] += a
            pn[0] += c
            p = pn
        return (p + [mp.mpf(0)] * M)[:M]

    ns = shift(nc)
    ds = shift(dc)
    if abs(ds[0]) > mp.mpf('1e-30'):
        r = [ns[0] / ds[0]]
        for m in range(1, M):
            s = ns[m]
            for l in range(1, m + 1):
                s -= ds[l] * r[m - l] if l < len(ds) else 0
            r.append(s / ds[0])
        return mp.mpf(0), r[:K + 2]
    d1 = ds[1:] + [mp.mpf(0)]
    r = [ns[0] / d1[0]]
    for m in range(1, M):
        s = ns[m]
        for l in range(1, m + 1):
            s -= d1[l] * r[m - l] if l < len(d1) else 0
        r.append(s / d1[0])
    return r[0], r[1:K + 3][:K + 2]


# ---------------- kernel jets K^(0), K^(1) + threshold series ----------------
class KernelJets:
    """K_P layers 0,1 everywhere in (1, WMAX]:
      |u| < vs : threshold moment-jet series (kn_jets, tail-bounded)
      u >= vs  : (K^(j)(w) - K0^(j) - K1^(j) u)/u^2 with K^(j) from the 3-layer
                 block jet transport of the exact box1eq DE (JetMarch cache).
    Gates: layer -1 of the box component ~ 0; K^(0) vs closed dilog
    (lbl3disp.box1); K^(1) vs pointwise parametric 2-fold jet; series-vs-
    transport crossover agreement at u = vs."""

    B_SINGS = [0, 1, mp.mpc(0, 1), mp.mpc(0, -1), mp.mpc(0, 2), mp.mpc(0, -2)]

    def __init__(self, VP, dps, guard, vs=sp.Rational(1, 64)):
        self.VP = VP
        self.dps = dps
        self.guard = guard
        self.vs_rat = sp.Rational(vs)
        self.vs = mp.mpf(self.vs_rat.p) / self.vs_rat.q
        tol = mp.mpf(10) ** (-(dps + guard))
        # threshold series length from geometric tail at |u|=vs (radius 1)
        NS = int((dps + guard) / (-mp.log10(self.vs))) + 10
        n0 = max(80, int(mp.mpf('1.35') * (dps + guard)))
        self.kn = kn_jets_refine(NS, dps, guard, n0, f"kn_jets(NS={NS})")
        # tail bound of the series at vs from trailing terms
        lastmax = max(max(abs(x) for x in kj.c) * self.vs ** (j - 2)
                      for j, kj in enumerate(self.kn[NS - 8:], start=NS - 8))
        r = self.vs * mp.mpf('1.3')
        self.ser_tail = lastmax * r / (1 - r)
        if self.ser_tail > tol * 100:
            raise RuntimeError(
                f"KernelJets: threshold-series tail {mp.nstr(self.ser_tail,3)} "
                f"> tol {mp.nstr(tol,3)} at vs={self.vs_rat}; raise NS -- fail-closed")
        self.NS = NS
        self.march = None

    def seed_and_march(self, WMAX, N, sf, verbose=True):
        b_sings = self.B_SINGS + [
            mp.mpc(mp.mpf(1) / 3, mp.sqrt(mp.mpf(8)) / 3),
            mp.mpc(mp.mpf(1) / 3, -mp.sqrt(mp.mpf(8)) / 3)]
        bj = block_de_json(self.VP.DATA['box1eq_de'], 3)
        de = self.VP.RatDE(bj, sp.Integer(0), b_sings)
        n_gl = max(100, int(mp.mpf('1.35') * (self.dps + self.guard)))
        mj = None
        prev = None
        tol = mp.mpf(10) ** (-(self.dps + self.guard))
        work = self.dps + self.guard + 20
        for it in range(5):
            cur = master_laurent_jets(5, n_gl, work)
            if prev is not None:
                agr = max(max(abs(a - b) for a, b in zip(x[1].c, y[1].c))
                          for x, y in zip(cur, prev))
                if agr <= tol:
                    mj = cur
                    if verbose:
                        print(f"    [certified] master seeds: GL n={n_gl} agreement "
                              f"{mp.nstr(agr,3)} <= tol {mp.nstr(tol,3)}", flush=True)
                    break
            prev = cur
            n_gl = int(n_gl * 3 // 2) + 8
        if mj is None:
            raise RuntimeError("KernelJets: master seed jets not converged -- fail-closed")
        seed = []
        for lay in range(3):                       # layers -1, 0, 1
            for m in mj:
                seed.append(laurent_to_layers(m)[lay])
        self.march = JetMarch(de, seed, mp.mpf(5), N, sf, comps=[0, 8, 16],
                              dps=self.dps, guard=self.guard)
        w_lo = 1 + self.vs / 2
        self.march.build(w_lo, WMAX)
        self.K01 = [J([self.kn[0].c[0], self.kn[0].c[1]]),
                    J([self.kn[1].c[0], self.kn[1].c[1]])]

    def to_record(self):
        """Checkpoint form of stage [1]: the certified threshold moments and
        the built kernel transport."""
        return {'stage': 1, 'NS': self.NS, 'vs': str(self.vs_rat), 'ser_tail': _num_str(self.ser_tail),
                'kn': [[_num_str(x) for x in kj.c] for kj in self.kn], 'march': self.march.to_record()}

    @classmethod
    def from_record(cls, VP, dps, guard, rec):
        """Constructor path for stage [1] from its checkpoint form: the
        threshold moments and the prebuilt march restored, nothing recomputed.
        The gates are not restored -- the caller re-runs them on the result."""
        self = cls.__new__(cls)
        self.VP = VP
        self.dps = dps
        self.guard = guard
        self.vs_rat = sp.Rational(rec['vs'])
        self.vs = mp.mpf(self.vs_rat.p) / self.vs_rat.q
        self.NS = int(rec['NS'])
        self.kn = [J([_num_val(x) for x in c]) for c in rec['kn']]
        self.ser_tail = _num_val(rec['ser_tail'])
        self.march = JetMarch.from_record(rec['march'])
        self.K01 = [J([self.kn[0].c[0], self.kn[0].c[1]]),
                    J([self.kn[1].c[0], self.kn[1].c[1]])]
        return self

    # ---- raw kernel layers ----
    def K_layer(self, w, j):
        """K^(j)(w), j = 0 or 1, from the transport cache (comp 8/16)."""
        return self.march.eval(w, 8 if j == 0 else 16)

    def K_layer_m1(self, w):
        return self.march.eval(w, 0)

    # ---- threshold-subtracted layers ----
    def KP(self, w, j):
        w = mp.mpf(w)
        u = w - 1
        if abs(u) < self.vs:
            acc = mp.mpf(0)
            for kj in reversed(self.kn[2:]):
                acc = acc * u + kj.c[j]
            return acc
        return (self.K_layer(w, j) - self.kn[0].c[j] - self.kn[1].c[j] * u) / (u * u)

    # ---- gates ----
    def gates(self, verbose=True):
        dps, g = self.dps, self.guard
        rep = {}
        # (a) box layer -1 ~ 0 along the range
        m1max = max(abs(self.K_layer_m1(w)) for w in
                    [mp.mpf('1.05'), mp.mpf(2), mp.mpf(8), mp.mpf(50), mp.mpf(1e6)])
        rep['layer_m1_max'] = m1max
        # (b) K^(0) vs closed dilog at spot points
        d0 = min(_agree(self.K_layer(w, 0), LD.box1(mp.mpf(-1), mp.mpf(-1) / 3,
                                                    mp.mpf(1), w, dps + g))
                 for w in [mp.mpf(3) / 2, mp.mpf(30)])
        rep['K0_vs_dilog_d'] = d0
        # (c) K^(1) vs pointwise parametric 2-fold jet at spot points
        dmin = mp.inf
        for w in [mp.mpf(3) / 2, mp.mpf(30)]:
            ptw = _box_jet_pointwise(w, dps, g)
            dmin = min(dmin, _agree(self.K_layer(w, 1), ptw.c[1]))
            dmin = min(dmin, _agree(self.K_layer(w, 0), ptw.c[0]))
        rep['K1_vs_2fold_d'] = dmin
        # (d) series-vs-transport crossover at u = vs
        wv = 1 + self.vs
        dser = min(_agree((self.K_layer(wv, j) - self.kn[0].c[j]
                           - self.kn[1].c[j] * self.vs) / self.vs ** 2,
                          sum(kj.c[j] * self.vs ** (m - 2)
                              for m, kj in enumerate(self.kn[2:], start=2)))
                   for j in (0, 1))
        rep['crossover_d'] = dser
        if verbose:
            print(f"    kernel gates: |box layer -1| <= {mp.nstr(m1max,3)}; "
                  f"K0 vs dilog {d0:.1f} d; K jets vs pointwise 2-fold {dmin:.1f} d; "
                  f"series/transport crossover {dser:.1f} d", flush=True)
        bar = dps - 8
        if not (d0 >= bar and dmin >= bar and dser >= bar
                and m1max < mp.mpf(10) ** (-(dps - 5))):
            raise RuntimeError(
                f"KernelJets gates FAILED (bar {bar} d): K0/dilog {d0:.1f}, "
                f"jets/2fold {dmin:.1f}, crossover {dser:.1f}, "
                f"layer-1 {mp.nstr(m1max,3)} -- fail-closed")
        return rep


def _agree(a, b):
    if a == b:
        return mp.inf
    if b == 0:
        return float(-mp.log10(abs(a - b)))
    return float(-mp.log10(abs(a - b) / abs(b)))


def _box_jet_pointwise(w, dps, guard):
    """Jet of K(w;e) by the direct parametric 2-fold (independent route,
    n-doubling until agreement)."""
    tol = mp.mpf(10) ** (-(dps + guard))
    work = dps + guard + 20
    c = mp.mpf(4) / 3
    n = max(80, int(mp.mpf('1.35') * (dps + guard)))
    prev = None
    for it in range(5):
        nodes, weights = k33lib.gl_nodes(n, work)
        half = mp.mpf('0.5')
        tb = [mp.mpf(0)] * ORD
        for t1_, w1 in zip(nodes, weights):
            x1 = half * (t1_ + 1)
            L = 1 - x1
            acc = [mp.mpf(0)] * ORD
            for t2_, w2 in zip(nodes, weights):
                x2 = L * half * (t2_ + 1)
                X = 1 - x1 - x2
                A = 1 + x1 * (w - 1) - c * x2 * X
                B = x1 + c * x2
                pds = pd_jets(A, B, X, mp.mpf(1), 3)
                for k in range(ORD):
                    acc[k] += w2 * pds[k]
            for k in range(ORD):
                tb[k] += w1 * acc[k] * (L * half)
        e = J([mp.mpf(0), mp.mpf(1)])
        cur = gamma_jet(mp.mpf(2)) / (J.const(1) + e) * J(tb) * half
        if prev is not None and max(abs(a - b) for a, b in zip(cur.c, prev.c)) <= tol:
            return cur
        prev = cur
        n = int(n * 3 // 2) + 8
    raise RuntimeError("_box_jet_pointwise: not converged -- fail-closed")


# ======================= assembly: c_-2, c_-1, c0 ============================
# Reference value for the final comparison: the eps^0 layer as recorded on
# 2026-07-05 by a different route (pole-subtracted Neville extrapolation of
# fixed-eps dispersive values over the grid eps = 2^-20..2^-25 at dps 150),
# whose two depth-varied runs agreed to 39.16 digits (trust ~39 d).
# Comparison only -- it enters no computation here:
I0_ORACLE = "-0.1466053426200695281803505147176612969729"

L_RAT = sp.Rational(2, 5)      # s0 - 1, s0 = 7/5 exact rational split
S0_RAT = sp.Rational(7, 5)


def run(dps, sabotage=None, verbose=True, checkpoint=None, resume=None):
    """Subtracted-integral eps^0 extraction (plus full exact pole layers).
    Returns dict with c_m2, c_m1, c0, gates, budget. Fail-closed everywhere.
    checkpoint=DIR saves stage [1] and stage [2] there as they complete
    (None: no saving); resume=DIR verifies that directory's pins and
    parameters BEFORE any computation (CheckpointRefused on any mismatch),
    restores the saved stages and re-enters at the first missing one; a
    resumed run keeps saving into the same directory unless checkpoint=DIR
    names another."""
    t00 = time.time()
    guard = 12
    wdps = dps + guard + 30
    mp.mp.dps = wdps
    VP = sys.modules.get('eval_row33') or __import__('eval_row33')
    tol = mp.mpf(10) ** (-(dps + guard))
    WMAX = mp.mpf(10) ** max(35, dps + 10)
    N_TAY = max(80, int(mp.mpf('0.95') * (dps + guard + 20)))
    SF = mp.mpf('0.28')
    L = mp.mpf(L_RAT.p) / L_RAT.q
    s0 = mp.mpf(S0_RAT.p) / S0_RAT.q
    lnL = mp.log(L)

    if verbose:
        print(f"row 33 eps^0 layer, subtracted-integral path: dps={dps} "
              f"guard={guard} WMAX=1e{int(mp.log10(WMAX))} N_TAY={N_TAY} "
              f"vs=1/64 s0=7/5", flush=True)

    # stage checkpoint store: opened (and verified) for resume before anything computes
    params = {'dps': dps, 'guard': guard, 'wdps': wdps, 'prec_bits': mp.mp.prec,
              'WMAX_log10': int(mp.log10(WMAX)), 'N_TAY': N_TAY, 'SF': mp.nstr(SF, 10),
              'vs': str(sp.Rational(1, 64)), 's0': str(S0_RAT), 'nsteps_guess': 400}
    src = None
    if resume is not None:
        src = StageStore.open(resume, params, StageStore.producer(VP))
        if verbose:
            print(f"    [resume] {resume}: MANIFEST pins verified for "
                  f"{len(src.manifest['stages'])} saved stage file(s), engine identity and run "
                  f"parameters match dps={dps}", flush=True)
    store = src
    if checkpoint is not None and (src is None or os.path.abspath(checkpoint) != os.path.abspath(src.dir)):
        store = StageStore(checkpoint, params, StageStore.producer(VP))
    if verbose and store is not None and src is None:
        print(f"    [checkpoint] stages [1] and [2] will be saved under {checkpoint}", flush=True)

    # [1] kernel jets: threshold moments + 3-layer box block transport
    t0 = time.time()
    if src is not None and src.has(1):
        rec = src.load(1)
        KJ = KernelJets.from_record(VP, dps, guard, rec)
        if verbose:
            print(f"    [resume] stage [1] restored from {src.manifest['stages']['1']['file']}: "
                  f"kn_jets(NS={KJ.NS}), march {KJ.march.stats['steps']} steps; kernel gates "
                  f"re-run below on the restored series", flush=True)
    else:
        KJ = KernelJets(VP, dps, guard)
        if sabotage == 'k2':
            KJ.kn[2] = J([KJ.kn[2].c[0] + mp.mpf('1e-30'), KJ.kn[2].c[1], KJ.kn[2].c[2]])
            print("    !! SABOTAGE: K2^(0) perturbed by 1e-30 (mutation test)", flush=True)
        KJ.seed_and_march(WMAX * mp.mpf('1.01'), N_TAY, SF, verbose)
    kg = KJ.gates(verbose)
    if store is not None and not (src is not None and src.has(1)):
        rec = KJ.to_record()
        if not KJ.march.record_matches(rec['march']):
            raise RuntimeError("stage [1] checkpoint: a stored series does not parse back to its "
                               "in-memory value -- not saved, fail-closed")
        rec['stats'] = dict(KJ.march.stats)
        rec['checks'] = _plain(kg)
        store.save(1, rec, time.time() - t0, verbose)
    if verbose:
        st = KJ.march.stats
        print(f"    [1] kernel: {time.time()-t0:.0f}s, march {st['steps']} steps "
              f"(halved {st['halved']}, N-raised {st['nraised']}), "
              f"step-bound acc {mp.nstr(KJ.march.bound_acc,3)}", flush=True)

    # [2] rho jets: closed 2-body + Frobenius + 2-layer sigvp block transport
    t0 = time.time()
    RJ = RhoJets(VP, dps, guard)
    if src is not None and src.has(2):
        RJ.attach_march(src.load(2))
        if verbose:
            print(f"    [resume] stage [2] restored from {src.manifest['stages']['2']['file']}: "
                  f"march {RJ.march.stats['steps']} steps; density gates re-run below", flush=True)
    else:
        RJ.build_tail(WMAX * mp.mpf('1.01'), N_TAY, SF)
    # gates: transported r_-1 vs closed 2/u - 1/w; r_0 series-vs-transport at 10.5
    rg = {}
    dm1 = min(_agree(RJ.rho_layers_tail(w)[0], 2 / (w - 1) - 1 / w)
              for w in [mp.mpf(12), mp.mpf(100), mp.mpf(1e6)])
    w105 = mp.mpf('10.5')
    dser = _agree(RJ.rho_layers_tail(w105)[1], RJ.rho_layers(w105)[1])
    rg['rm1_transport_d'] = dm1
    rg['r0_series_vs_transport_d'] = dser
    bar = dps - 8
    if dm1 >= bar and dser >= bar and store is not None and not (src is not None and src.has(2)):
        rec = RJ.to_record()                 # saved only when the density gates pass
        if not RJ.march.record_matches(rec['march']):
            raise RuntimeError("stage [2] checkpoint: a stored series does not parse back to its "
                               "in-memory value -- not saved, fail-closed")
        rec['stats'] = dict(RJ.march.stats)
        rec['checks'] = _plain(rg)
        store.save(2, rec, time.time() - t0, verbose)
    if verbose:
        st = RJ.march.stats
        print(f"    [2] rho: {time.time()-t0:.0f}s, march {st['steps']} steps, "
              f"r_-1 transport vs closed {dm1:.1f} d, r_0 series/transport {dser:.1f} d",
              flush=True)
    if not (dm1 >= bar and dser >= bar):
        raise RuntimeError(f"rho gates FAILED (bar {bar} d): r_-1 {dm1:.1f}, "
                           f"r_0 {dser:.1f} -- fail-closed")

    # [3] finite integrals (quad_refine, fail-closed)
    t0 = time.time()
    vs = KJ.vs

    def q0(u):
        return (2 + u) / (1 + u) * KJ.KP(1 + u, 0) - 2 * KJ.kn[2].c[0]

    def q1(u):
        return (2 + u) / (1 + u) * (KJ.KP(1 + u, 1)
                                    + mp.log1p(u) * KJ.KP(1 + u, 0)) \
            - 2 * KJ.kn[2].c[1]

    certs = {}

    def QR(f, iv, label, scale=1):
        c = quad_refine(f, iv, dps, guard=guard, scale=scale, full_output=True)
        certs[label] = c
        if verbose:
            print("    " + c.bound_line(label), flush=True)
        return c.value

    J0 = QR(lambda u: q0(u) / u, [mp.mpf(0), vs, L], "J0")
    J1 = QR(lambda u: (-2 * mp.log(u) * q0(u) + q1(u)) / u,
            [mp.mpf(0), vs, L], "J1")
    B0 = QR(lambda u: -(u / (1 + u)) * KJ.KP(1 + u, 0) / 6, [mp.mpf(0), L], "B0")

    def rm1(w):
        return 2 / (w - 1) - 1 / w

    def r0f(w):
        return RJ.rho_layers(w)[1] if w <= RJ.wsw else RJ.rho_layers_tail(w)[1]

    def fTm1(w):
        return rm1(w) * KJ.KP(w, 0)

    def fT0(w):
        return rm1(w) * KJ.KP(w, 1) + r0f(w) * KJ.KP(w, 0)

    Tm1 = QR(fTm1, [s0, mp.mpf(9), RJ.wsw, mp.mpf(200)], "Tm1.body")
    T0 = QR(fT0, [s0, mp.mpf(9), RJ.wsw, mp.mpf(200)], "T0.body")

    # far tail via w = 200 + 200(1+t)/(1-t), truncated at WMAX + measured bound
    R = (WMAX - 200) / 200
    tW = (R - 1) / (R + 1)

    def mapped(f):
        def g(t):
            w = 200 + 200 * (1 + t) / (1 - t)
            return f(w) * 400 / (1 - t) ** 2
        return g

    Tm1 += QR(mapped(fTm1), [mp.mpf(-1), tW], "Tm1.far")
    T0 += QR(mapped(fT0), [mp.mpf(-1), tW], "T0.far")

    # beyond-WMAX remainder bound from measured local decay power
    tail_rem = mp.mpf(0)
    for f, lab in ((fTm1, 'Tm1'), (fT0, 'T0')):
        f2 = abs(f(WMAX))
        f1 = abs(f(WMAX / 4))
        if f2 == 0:
            continue
        p = float(mp.log(f1 / f2) / mp.log(4))
        if p < 1.5:
            raise RuntimeError(f"far-tail decay power {p:.2f} < 1.5 for {lab} "
                               f"-- WMAX tail not certifiably negligible")
        tail_rem += f2 * WMAX / (p - 1)
    if verbose:
        print(f"    [3] integrals: {time.time()-t0:.0f}s; beyond-WMAX remainder "
              f"bound {mp.nstr(tail_rem,3)}", flush=True)

    # [4] exact Gamma-series pole product + assembly
    e = J([mp.mpf(0), mp.mpf(1)])
    gh = ghat_jet()
    PT = gh * J([KJ.kn[2].c[0], KJ.kn[2].c[1], KJ.kn[2].c[2]]) * jpow(L, -2 * e)
    c_m2 = -PT.c[0]
    c_m1 = -PT.c[1] + J0 + Tm1
    c0 = -PT.c[2] + gh.c[1] * J0 + J1 + B0 + T0

    budget = sum(c.err_bound for c in certs.values()) + tail_rem \
        + KJ.march.bound_acc + RJ.march.bound_acc + KJ.ser_tail
    out = {'c_m2': c_m2, 'c_m1': c_m1, 'c0': c0, 'budget': budget,
           'kernel_gates': kg, 'rho_gates': rg, 'wall_s': time.time() - t00}

    # [5] smoke comparisons of the exact pole layers with the recorded
    # 50-digit strings (comparison only; those strings are honest to about
    # 45.8 / 39.6 digits, so at dps >= 60 they, not the layers, set the figure)
    o2 = mp.mpf(VP.ORC['laurent']['eps-2'])
    o1 = mp.mpf(VP.ORC['laurent']['eps-1'])
    d2 = _agree(c_m2, o2)
    d1 = _agree(c_m1, o1)
    bar2 = min(44, dps - 8)
    bar1 = min(38, dps - 8)
    if verbose:
        print(f"    [smoke] c_-2 vs stored string {d2:.2f} d (bar {bar2}; string "
              f"honest ~45.8); c_-1 {d1:.2f} d (bar {bar1}; string honest ~39.6)",
              flush=True)
    if not (d2 >= bar2 and d1 >= bar1):
        raise SystemExit(f"BLOCKED: exact pole layers disagree with the recorded "
                         f"strings: c_-2 {d2:.2f} d (bar {bar2}), c_-1 {d1:.2f} d "
                         f"(bar {bar1}) -- stopping")
    out['smoke_d2'] = d2
    out['smoke_d1'] = d1

    # [6] the final comparison: c0 vs the recorded grid-route reference (trust ~39 d)
    d0 = _agree(c0, mp.mpf(I0_ORACLE))
    bar0 = min(38, dps - 10)
    out['oracle_d'] = d0
    out['oracle_bar'] = bar0
    if verbose:
        mp.mp.dps = min(dps, 60)
        print(f"  c_-2 = {mp.nstr(c_m2, min(dps, 50))}")
        print(f"  c_-1 = {mp.nstr(c_m1, min(dps, 50))}")
        print(f"  c0   = {mp.nstr(c0, min(dps, 50))}   (SUBTRACTED FINITE INTEGRAL)")
        mp.mp.dps = wdps
        print(f"  error budget (quad certs + march tail bounds + series tails + "
              f"beyond-WMAX) <= {mp.nstr(budget, 3)}")
        print(f"  ORACLE GATE: c0 vs the recorded reference I0 (grid route, 2026-07-05): {d0:.2f} d "
              f"(bar {bar0}; reference trust ~39 d)   [{out['wall_s']:.0f}s total]",
              flush=True)
    if d0 < bar0:
        print(f"GATE FAILED: {d0:.2f} d < bar {bar0} d", flush=True)
        raise SystemExit(1)
    return out


def main():
    import argparse
    ap = argparse.ArgumentParser(description="row 33 eps^0 subtracted-integral path (engine CLI; the served, pinned entry point is ../lbl3vp-eps0-evaluate.py)")
    ap.add_argument('--dps', type=int, default=60)
    ap.add_argument('--sabotage-k2', action='store_true',
                    help="mutation test: perturb K2^(0) by 1e-30 (gate must FAIL, rc=1); implies --no-checkpoint")
    ap.add_argument('--json', default=None, help="write results JSON here")
    ap.add_argument('--checkpoint', metavar='DIR', default=None,
                    help="save stage [1] and stage [2] under DIR as they complete (default: ./row33_ckpt_dps<dps>)")
    ap.add_argument('--no-checkpoint', action='store_true', help="save nothing")
    ap.add_argument('--resume', metavar='DIR', default=None,
                    help="verify DIR's MANIFEST pins and parameters, restore its saved stages and re-enter at the "
                         "first missing one; saving continues into DIR (--checkpoint may not name another directory)")
    args = ap.parse_args()
    if args.resume and args.sabotage_k2:
        ap.error("--sabotage-k2 needs a fresh run (the perturbation is applied when the moments are computed)")
    if args.resume and args.checkpoint and os.path.abspath(args.checkpoint) != os.path.abspath(args.resume):
        ap.error("--checkpoint with --resume: a resumed run saves into the resume directory")
    ckpt = None if (args.no_checkpoint or args.sabotage_k2) else (args.checkpoint or args.resume or f"./row33_ckpt_dps{args.dps}")
    try:
        out = run(args.dps, sabotage=('k2' if args.sabotage_k2 else None), checkpoint=ckpt, resume=args.resume)
    except CheckpointRefused as e:
        print(str(e), flush=True)
        sys.exit(e.exit_code)
    if args.json:
        import json as _json
        mp.mp.dps = args.dps + 20
        _json.dump({k: (mp.nstr(v, args.dps + 15) if isinstance(v, (mp.mpf, mp.mpc))
                        else (v if not isinstance(v, dict) else
                              {kk: (mp.nstr(vv, 20) if isinstance(vv, (mp.mpf, mp.mpc))
                                    else vv) for kk, vv in v.items()}))
                    for k, v in out.items()}, open(args.json, 'w'), indent=1)
        print(f"[json] {args.json}")


if __name__ == '__main__':
    main()