#!/usr/bin/env python3
r"""LBL3VP (vacuum-polarisation-dressed QED light-by-light box, 3 loops):
LIVE fully-analytic evaluation at s=-1, t=-1/3, m^2=1 — round 2, compliant
final form (no numeric seeds, no node caches, no DE-transported densities).

FINAL FORM.  The integral is the one-fold convolution in the internal
photon mass w' (the VP-channel virtuality),

  I_VP(eps) = (1/pi) Integral_1^inf  rho_VP(w';eps) K_P(w';eps) dw',

in which EVERY factor is now an analytic closed form or a one-fold iterated
integral over closed-form kernels, exact in d = 4-2eps:

  * two-body density (all w' > 1, exact in d, elementary):
      rho_2b(w') = -pi G(eps) [ -(1-2eps)(w'+1)/(w'-1) + eps(w'-1)/6 ]
                   (w'-1)^{-2eps} w'^{eps-1},
      G(eps) = Gamma(1+eps)Gamma(1-eps)/(eps Gamma(2-2eps));
  * three-body density (w' > 9) — NEW closed form of this round, the
    sunrise-elliptic Gamma_1(6) piece as an exact-in-d PERIOD ONE-FOLD:
      rho_3b(w') = pi c1(eps)^2 w'^{eps-1}
                   Integral_4^{(sqrt(w')-1)^2} dq  (q-4)^{1/2-eps} q^{-5/2}
                     [ ((sqrt(w')-1)^2-q) ((sqrt(w')+1)^2-q) ]^{1/2-eps},
      c1(eps) = Gamma(1-eps)/Gamma(2-2eps).
    Its branch points {0, 4, (sqrt(w')-+1)^2} are the equal-mass-sunrise
    quartic — the Gamma_1(6) curve (Bloch–Vanhove); at eps=0 the integrand
    reduces to the R_2 elliptic integrand of SIGVP_CLOSED_FORM.json, i.e.
    rho_3b^{(0)} is a Gamma_1(6) period of that curve. This replaces round
    1's DE-TRANSPORTED (rejected-form) Sigma_VP density. Derivation:
    dispersing the VP sub-bubble B(q^2) inside Sigma[1,2,1,1] and partial-
    fractioning 1/((q^2)^2(q^2-sigma)) reduces Sigma to bubble integrals,
    whose cuts are elementary powers exact in d; the 2-body formula above
    is re-obtained and the 3-body one-fold follows. Validated per eps-order
    (eps^0..eps^4) at 58-108 digits against banked 2-loop AMFlow probes at
    w=12, 50 (<archive>/lbl3vp/r2-eichler/v1b_orders.py).
  * kernel K(w';eps) (the 1-loop box with photon-line mass^2 = w'), exact
    in eps and SEED-FREE. Its defining analytic form is the exact 2-fold
    Feynman-parametric integral (one of the three parametric integrations
    done analytically, elementary integrand, exact in eps):
      K(w') = Gamma(2+eps) Int_0^1 dx1 Int_0^{1-x1} dx2
              [ A^{-1-eps} - (A+B X)^{-1-eps} ] / ((1+eps) B),
      A = w' x1 + (1-x1) - u x2(1-x1-x2),  B = -s x1 + u x2,  X = 1-x1-x2.
    Evaluating this 2-fold at every quadrature node would be slow, so the
    script evaluates it ONCE per eps (the anchor, at w0 = 200) and carries
    it everywhere with the exact rational IBP connection A(w',d) (shipped
    in lbl3vp-data.json, reconstructed from 46 exact-rational Kira samples,
    0 mismatches) by radius-controlled Taylor-series stepping — the
    standard series evaluation of the DE-defined iterated integral. All
    subsector masters are Gamma-function closed forms (tadpoles) or
    elementary one-folds exact in eps (bubbles/triangles). Threshold
    Taylor coefficients K_n(eps) come from an arc-VoP Cauchy circle
    |w'-1| = 3/8 (loop-closure defect printed, ~1e-60). In addition the
    kernel is RE-DERIVED each run from the w'->infinity VACUUM boundary
    (K -> -Tri_u(eps)/w') via the variation-of-parameters one-fold
      K(w') = H(w') [ K(W)/H(W) + Integral_W^{w'} (sum_j A_0j M_j)(v)/H(v) dv ],
      H(w') = (1+w'^2)^{-(1+2eps)/2} = exp Int A_00,
    and the two independent constructions are compared LIVE (printed
    "vacuum VoP vs anchor", ~30-32 d at the defaults; the vacuum fold's
    quadrature carries a small ~1/eps-scaling numeric residual, which is
    why the parametric anchor is the primary). Round 1's 1-loop AMFlow
    boundary Laurent seeds and threshold series are now LIVE CROSS-CHECKS
    only, never inputs.

Assembly (threshold peel of the u^{-2eps} mode, tanh-sinh panels, one
Richardson step) is unchanged from round 1. The gate compares, at fixed
rational eps, against the held-out 3-loop AMFlow oracle grid and recomputes
agreement live as -log10 |I-I_AMF|/|I_AMF|.

Laurent pole layers eps^-2, eps^-1 are also computed at runtime from the
closed-form dilog kernel (round-1 machinery, unchanged); since 2026-07-05 the
1-loop constant K_2^{(1)} (order-eps threshold Taylor coefficient) is ALSO
derived at runtime from the exact parametric 2-fold (K2_threshold_eps01 —
threshold-analytic differentiation under the integral, no seeds, no fit), so
the artifact ships NO numeric inputs: the stored AMFlow string Kn['2']['1']
demotes to a held-out cross-check like the rest of Kn.  The layer depth knobs
(fit degree nfit, quad dps, kernel dps) are wired to --dps (frozen at the
shipped demo values for dps <= 60), so the printed layer digits track the
requested precision instead of the old fixed-depth ~45/39.6 d demo caps.

Shipped data files (same directory, provenance inside each file):
  lbl3vp-data.json          box1eq_de: the exact rational connection A(w',d)
                            (analytic formula — the only data the fixed-eps
                            computation uses); box1eq_bnd_w5 / Kn: 1-loop
                            AMFlow values, now cross-checks only;
                            sigvp_de / sigvp_bnd_w5: round-1 transport data,
                            RETIRED from the computation (kept for
                            provenance; the density is closed-form now).
  lbl3vp-oracle-amflow.json held-out oracle: independent 3-loop AMFlow solve
                            of the full 8-propagator integral (fixed-eps grid
                            k=4..18 + Laurent layers). Comparison only.

Default demo: eps in {2^-6, 2^-9, 2^-13} at dps 60, one forked worker per
eps point (wall time printed; measured ~9 min wall / ~6 min CPU per point
on the build host under load; ~25 min serial with VP_JOBS=1). Round-1's
archived DISPERSION gate (31.70-32.79 d) is superseded by this seed-free
form; the live gate prints its own digits (measured at dps 60: 37.7 d at 2^-6,
40.1 d at 2^-13; raise --dps for more).

Interface (evaluate at a different point / precision, no code editing):
  python3 lbl3vp-evaluate.py                      default gate demo
  python3 lbl3vp-evaluate.py --point 2^-8         another eps (gated live if
                                                  on the k=4..18 oracle grid)
  python3 lbl3vp-evaluate.py --point 1/300 --dps 90
  python3 lbl3vp-evaluate.py --no-fastcache       force the full live run on
                                                  the banked demo grid
  python3 lbl3vp-evaluate.py --double [--point 2^-10] [--dps 60]
                                                  dps-doubling demo: one gate
                                                  point at dps and 2*dps
Domain limits (stated, not faked):
  * Kinematics are FIXED at s=-1, t=-1/3, m^2=1: the shipped rational
    connection A(w',d) was reconstructed on this kinematic slice, and the
    held-out 3-loop oracle is specific to it. Other (s,t) require
    regenerating the connection (the production pipeline,
    <archive>/phys_lbl3_qed_parents/vp/) — stated gap, so the runtime inputs
    of the interface are (eps, dps).
  * eps: rational, 0 < eps <= 1/16. Everything is exact in eps (round 1's
    eps^20 seed-truncation cap is GONE — there are no seeds); the working-
    precision overhead ~ 2*eps*log10(W) of the vacuum boundary is added
    internally. eps <= 1/16 keeps that overhead small and matches the
    below-threshold oracle grid.
  * precision: arbitrary — raising --dps raises every internal order
    (Taylor order, VoP panels, circle points, quadrature levels, inner
    one-fold nodes all scale with dps) until the 320-digit oracle strings
    cap the PRINTED comparison. No stored numeric table caps the value.

FAST-START CACHE (2026-07-06): lbl3vp-fastcache.json in
this directory holds the certified demo-grid values (I_VP at the banked eps
grid + the Laurent pole layers) BANKED from one full live run of THIS script
(provenance + sha pins inside the file).  Default behavior when every
requested eps is banked and the requested dps <= the banked depth: the banked
values print instantly, the held-out AMFlow-grid gates are re-asserted NOW on
the banked strings (any mismatch, incl. a 1e-30 mutation of a cached string,
exits nonzero), and a live cheap cross-check runs (exact parametric 2-fold
K_2^(1) vs its banked value AND the stored held-out AMFlow string).  The
cache is a fast start ONLY -- the live machinery in this file remains the
definition: any other eps, deeper dps, --double, --no-fastcache, or a
VP_TAYLOR_N/VP_LMAX/VP_CIRCLE_M/VP_STEP_FRAC/VP_VS/VP_NGL_INF env override
runs it unchanged.

Env knobs: VP_EPS_KS="6,9,13"  VP_DPS=70  VP_TAYLOR_N  VP_LMAX  VP_JOBS
           VP_STEP_FRAC=0.28   VP_VS      VP_CIRCLE_M  VP_NGL_INF
Requires python3 + mpmath + sympy (connection parsing only).
"""
import os, json, time, bisect, math
import mpmath as mp
import sympy as sp

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = json.load(open(os.path.join(HERE, 'lbl3vp-data.json')))
ORC = json.load(open(os.path.join(HERE, 'lbl3vp-oracle-amflow.json')))

# ---- kinematic point (fixed by the shipped connection; see docstring). ----
def kin():
    return mp.mpf(-1), mp.mpf(-1) / 3, mp.mpf(1)


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


# ============== fast-start cache (2026-07-06) =================================
# Doctrine: cache = sha-pinned + compare-gated fast start; the live machinery
# in this file remains the DEFINITION.  lbl3vp-fastcache.json holds the
# certified demo-grid values (I_VP at eps=2^-6,2^-9,2^-13 at the banked dps,
# plus the Laurent pole layers) BANKED from one full live run.  Default
# behavior when every requested eps is banked and --dps <= the banked depth:
# banked values print instantly, the held-out AMFlow-grid gates are
# re-asserted NOW on the banked strings (any mismatch, incl. a 1e-30 mutation
# of a cached string, exits nonzero), and a live cheap cross-check runs (the
# exact parametric 2-fold K_2^(1), vs BOTH its banked value and the stored
# held-out AMFlow string).  Any other eps, deeper --dps, --double,
# --no-fastcache, or VP_TAYLOR_N/VP_LMAX/VP_CIRCLE_M/VP_STEP_FRAC/VP_VS/
# VP_NGL_INF env override runs the live machinery unchanged.
FASTCACHE_PATH = os.path.join(HERE, 'lbl3vp-fastcache.json')
_FC_ENV_OVERRIDES = ('VP_TAYLOR_N', 'VP_LMAX', 'VP_CIRCLE_M', 'VP_STEP_FRAC',
                     'VP_VS', 'VP_NGL_INF')
_FC_PIN_FILES = ('lbl3vp-evaluate.py', 'lbl3vp-data.json',
                 'lbl3vp-oracle-amflow.json')


def _fc_sha(fname):
    import hashlib
    with open(os.path.join(HERE, fname), 'rb') as fh:
        return hashlib.sha256(fh.read()).hexdigest()


def _fc_str_sha(s):
    import hashlib
    return hashlib.sha256(s.encode()).hexdigest()


def _fc_eligible(args, DPS, pts):
    live_note = ("live machinery run (the definition) -- expect the full "
                 "wall (~7-8 min for the default demo at dps 60; "
                 "~2-6 min per eps at lower dps)")
    if args.no_fastcache or args.fastcache_bank:
        return None, None
    if args.double:
        return None, f"--double is a live demo by design; {live_note}"
    envs = [k for k in _FC_ENV_OVERRIDES if k in os.environ]
    if envs:
        return None, f"env override {envs} set; {live_note}"
    if not os.path.exists(FASTCACHE_PATH):
        return None, f"no lbl3vp-fastcache.json; {live_note}"
    try:
        with open(FASTCACHE_PATH) as fh:
            fc = json.load(fh)
    except Exception as ex:
        return None, f"fast cache unreadable ({ex!r}); {live_note}"
    if DPS > fc['cached_dps']:
        return None, (f"requested dps {DPS} exceeds the banked depth "
                      f"{fc['cached_dps']}; {live_note}")
    missing = [eps_label(ee) for ee in pts if eps_label(ee) not in fc['points']]
    if missing:
        return None, (f"requested eps {missing} not in the banked grid "
                      f"{sorted(fc['points'])}; {live_note}")
    if args.point is None and 'pole_layers' not in fc:
        return None, f"pole layers not banked; {live_note}"
    return fc, None


def _fc_fail(msg):
    raise RuntimeError(f"[fast-cache] COMPARE GATE FAIL: {msg} -- refusing "
                       "to serve the cache; rerun with --no-fastcache for "
                       "the live machinery (the definition)")


def _fc_gate(label, d_now, d_banked, floor):
    ok = (float(d_now) >= floor) and (abs(float(d_now) - float(d_banked)) <= 0.05)
    print(f"    [fast-cache gate] {label}: {float(d_now):.4f} d "
          f"(banked {d_banked:.4f} d, floor {floor} d) -- "
          f"{'OK' if ok else 'FAIL'}")
    if not ok:
        _fc_fail(f"{label}: recomputed {float(d_now):.4f} d vs banked "
                 f"{d_banked:.4f} d (floor {floor} d)")


def _fc_fast_path(fc, args, DPS, pts, t_start):
    prov = fc['_provenance']
    mp.mp.dps = max(80, DPS + 20)
    print("LBL3VP FAST-CACHE mode -- banked certified demo values "
          "(sha-pinned, compare-gated).")
    print(f"  cache: {os.path.basename(FASTCACHE_PATH)}  banked "
          f"{prov['generated_utc']} by: {prov['generator_cmd']}")
    print(f"  banked live wall {prov['wall_s']:.0f}s; this fast start serves "
          f"the banked eps grid at dps <= {fc['cached_dps']} "
          f"(requested {DPS})")
    print("  cache = fast start ONLY; the live machinery in this file is "
          "the definition (--no-fastcache)\n")
    # --- [p1] integrity pins ---------------------------------------------
    for fname, want in fc['_pins'].items():
        if _fc_sha(fname) != want:
            _fc_fail(f"sha256 pin mismatch on {fname} (definition changed "
                     "since banking; cache is STALE)")
    npins = 0
    for sect in list(fc['points'].values()) + (
            [fc['pole_layers']] if 'pole_layers' in fc else []):
        for key, want in sect['_string_pins'].items():
            if _fc_str_sha(sect[key]) != want:
                _fc_fail(f"cached string '{key}' fails its sha pin (mutated)")
            npins += 1
    print(f"[p1] integrity: {len(fc['_pins'])} file pins + {npins} "
          f"cached-string pins OK")
    # --- [p2] held-out oracle gates recomputed NOW ------------------------
    print("[p2] held-out AMFlow-grid gates recomputed NOW from the banked "
          "strings vs the oracle JSON:")
    for ee in pts:
        lab = eps_label(ee)
        pt = fc['points'][lab]
        val = mp.mpf(pt['I_VP'])
        print(f"  eps = {lab}:")
        print(f"    I_VP    = {mp.nstr(val, 40)}   (banked; live diag at bank "
              f"time: closure {pt['closure']:.1e}, quad_agree "
              f"{pt['quad_agree']:.1f} d)")
        ostr = oracle_for(ee)
        if ostr is None:
            _fc_fail(f"banked eps {lab} has no held-out oracle (should not "
                     "have been banked)")
        _fc_gate(f"I_VP({lab}) vs held-out 3-loop AMFlow grid",
                 agree_d(val, mp.mpf(ostr)), pt['banked_d'], 30)
    if args.point is None:
        pl = fc['pole_layers']
        I2c, I1c = mp.mpf(pl['I_em2']), mp.mpf(pl['I_em1'])
        print("  Laurent pole layers (banked):")
        print(f"    eps^-2: -K_2^(0)             = {mp.nstr(I2c, 36)}")
        _fc_gate("eps^-2 layer vs shipped oracle",
                 agree_d(I2c, mp.mpf(ORC['laurent']['eps-2'])),
                 pl['banked_d_em2'], 30)
        print(f"    eps^-1: -K_2^(1)+B[s_-1]^(0) = {mp.nstr(I1c, 36)}")
        _fc_gate("eps^-1 layer vs shipped oracle",
                 agree_d(I1c, mp.mpf(ORC['laurent']['eps-1'])),
                 pl['banked_d_em1'], 30)
    # --- [p3] live cheap cross-check: exact 2-fold K_2^(1) ----------------
    print("[p3] live cheap cross-check (machinery exercised NOW):")
    t0 = time.time()
    kdps_live = 45
    _, K2_1_live = K2_threshold_eps01(kdps_live)
    d_bank = agree_d(K2_1_live, mp.mpf(fc['K2_1_bank']))
    d_held = agree_d(K2_1_live, mp.mpf(DATA['Kn']['2']['1']))
    bar = kdps_live - 5
    print(f"    K_2^(1) exact parametric 2-fold at dps {kdps_live} "
          f"({time.time()-t0:.1f}s): vs banked {d_bank:.1f} d, vs stored "
          f"held-out AMFlow string {d_held:.1f} d >= bar {bar} d")
    if not (d_bank >= bar and d_held >= bar):
        _fc_fail(f"live K_2^(1) cross-check: {d_bank:.1f}/{d_held:.1f} d "
                 f"< bar {bar} d")
    print(f"\ntotal wall time: {time.time()-t_start:.0f}s")
    print("PASS (fast-cache) -- banked values re-gated vs the held-out "
          "AMFlow grid + live K_2^(1) cross-check; for the full live run "
          "use --no-fastcache")
    raise SystemExit(0)


def _fc_bank(fc_points, pole, K2_1_str, DPS, wall_s, cmd):
    import datetime
    fc = {
        '_doc': ("fast-start cache for lbl3vp-evaluate.py: banked from ONE "
                 "full live run (command below).  The live machinery is the "
                 "definition; this file only fast-starts the banked demo "
                 "grid and is refused on any sha-pin or compare-gate "
                 "mismatch.  Regenerate: --fastcache-bank."),
        '_provenance': {
            'generated_utc': datetime.datetime.now(
                datetime.timezone.utc).isoformat(),
            'generator_cmd': cmd,
            'wall_s': wall_s,
            'source': '<archive>/axis3_wave/disp-fastcache/',
        },
        '_pins': {f: _fc_sha(f) for f in _FC_PIN_FILES},
        'cached_dps': DPS,
        'K2_1_bank': K2_1_str,
        'points': {},
    }
    for lab, pt in fc_points.items():
        pt = dict(pt)
        pt['_string_pins'] = {'I_VP': _fc_str_sha(pt['I_VP'])}
        fc['points'][lab] = pt
    if pole is not None:
        pole = dict(pole)
        pole['_string_pins'] = {k: _fc_str_sha(pole[k])
                                for k in ('I_em2', 'I_em1')}
        fc['pole_layers'] = pole
    with open(FASTCACHE_PATH, 'w') as fh:
        json.dump(fc, fh, indent=1)
    print(f"[fast-cache] BANKED {FASTCACHE_PATH} (depth {DPS} d class, "
          f"{len(fc['points'])} eps points, {len(fc['_pins'])} file pins)")


# =====================================================================
# Taylor-series stepping of the exact rational connection (unchanged from
# round 1; this is the standard series evaluation of the DE-defined
# iterated integral, seeded now by the analytic vacuum boundary).
# =====================================================================
class RatDE:
    def __init__(self, de_json, eps_rat, sings):
        self.masters = [tuple(m) for m in de_json['masters']]
        self.n = len(self.masters)
        W, D = sp.symbols('w d')
        d_sym = sp.Rational(4) - 2 * eps_rat
        self.funcs = [[None] * self.n for _ in range(self.n)]
        for i in range(self.n):
            for j in range(self.n):
                e = sp.sympify(de_json['A'][i][j], locals={'w': W, 'd': D})
                if e == 0:
                    continue
                e = sp.cancel(e.subs(D, d_sym))
                num, den = sp.fraction(e)
                nc = [mp.mpf(sp.Rational(c).p) / mp.mpf(sp.Rational(c).q)
                      for c in sp.Poly(num, W).all_coeffs()]
                dc = [mp.mpf(sp.Rational(c).p) / mp.mpf(sp.Rational(c).q)
                      for c in sp.Poly(den, W).all_coeffs()]
                self.funcs[i][j] = (nc, dc)
        self.sings = [mp.mpc(s) for s in sings]

    def nearest(self, w):
        w = mp.mpc(w)
        return min(abs(w - s) for s in self.sings)

    def _shift(self, coefs, w0):
        p = [mp.mpc(coefs[0])]
        for c in coefs[1:]:
            pn = [mp.mpc(0)] * (len(p) + 1)
            for m, a in enumerate(p):
                pn[m] += a * w0
                pn[m + 1] += a
            pn[0] += c
            p = pn
        return p

    def step(self, M, w0, h, N):
        n = self.n
        ser = [[] for _ in range(N + 1)]
        for i in range(n):
            for j in range(n):
                cd = self.funcs[i][j]
                if cd is None:
                    continue
                nc, dc = cd
                ph = self._shift(nc, w0)
                qh = self._shift(dc, w0)
                q0 = qh[0]
                r = [mp.mpc(1) / q0]
                for m in range(1, N + 1):
                    s = mp.mpc(0)
                    for l in range(1, min(m, len(qh) - 1) + 1):
                        s += qh[l] * r[m - l]
                    r.append(-s / q0)
                for m in range(N + 1):
                    s = mp.mpc(0)
                    for l in range(min(m, len(ph) - 1) + 1):
                        s += ph[l] * r[m - l]
                    if s != 0:
                        ser[m].append((i, j, s))
        C = [list(M)]
        for m in range(N):
            s = [mp.mpc(0)] * n
            for l in range(m + 1):
                cv = C[m - l]
                for (ri, ci, a) in ser[l]:
                    s[ri] += a * cv[ci]
            inv = mp.mpf(1) / (m + 1)
            C.append([x * inv for x in s])
        v = [mp.mpc(0)] * n
        hp = mp.mpc(1)
        for m in range(N + 1):
            cm = C[m]
            for i in range(n):
                v[i] += cm[i] * hp
            hp *= h
        return v, C

    def path(self, M0, w_from, targets, N, sf, deflect=None):
        n = self.n
        M = list(M0)
        w = mp.mpc(w_from)
        out = {}
        ti = 0
        nT = len(targets)
        end = mp.mpc(targets[-1])
        sgn = 1 if end.real > w.real else -1
        wps = []
        for s in sorted(deflect or [], key=lambda x: sgn * x):
            if min(w.real, end.real) < s < max(w.real, end.real):
                r = mp.mpf('0.35')
                wps += [mp.mpc(s - sgn * r, 0), mp.mpc(s, r), mp.mpc(s + sgn * r, 0)]
        wps.append(end)
        for wp in wps:
            while abs(w - wp) > mp.mpf('1e-200'):
                d = self.nearest(w)
                hmax = sf * d
                dirn = wp - w
                stp = dirn if abs(dirn) <= hmax else dirn / abs(dirn) * hmax
                Mnew, C = self.step(M, w, stp, N)
                w_next = w + stp
                if abs(w.imag) < mp.mpf('1e-50') and abs(w_next.imag) < mp.mpf('1e-50') and stp != 0:
                    while ti < nT:
                        t = mp.mpc(targets[ti])
                        frac = ((t - w) / stp).real
                        if frac < -mp.mpf('1e-50') or frac > 1 + mp.mpf('1e-50'):
                            break
                        hh = t - w
                        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 *= hh
                        out[targets[ti]] = v
                        ti += 1
                M = Mnew
                w = w_next
            w = wp
        while ti < nT:
            out[targets[ti]] = list(M)
            ti += 1
        return out, M


def vs_for(DPS):
    """Half-width of the threshold Taylor patch |w'-1| < vs. The Taylor
    coefficients are computed LIVE (arc-VoP Cauchy circle, radius 0.35), so
    unlike round 1 there is no shipped-series truncation cap; vs only sets
    where the series (convergence radius 1) takes over from the stepped
    kernel. VP_VS overrides."""
    env = os.environ.get('VP_VS')
    if env:
        return mp.mpf(env)
    return mp.mpf('0.02') if DPS <= 90 else mp.mpf('0.008')


def load_laur(block):
    """Shipped boundary Laurent data {idx: {order: [re,im]}} -> mpc dict.
    (Round 2: cross-check data only.)"""
    out = {}
    for idx, orders in block.items():
        key = tuple(int(x) for x in idx.split(','))
        out[key] = {int(o): mp.mpc(mp.mpf(v[0]), mp.mpf(v[1])) for o, v in orders.items()}
    return out


def laur_at(l, eps, kmin, kmax):
    return sum(l.get(k, mp.mpc(0)) * eps ** k for k in range(kmin, kmax + 1))


def ts_nodes(a, b, level, dps):
    old = mp.mp.dps
    mp.mp.dps = dps + 20
    ts = mp.calculus.quadrature.TanhSinh(mp.mp)
    prec = mp.mp.prec
    pts = set()
    for deg in range(1, level + 1):
        for x, w in ts.get_nodes(mp.mpf(a), mp.mpf(b), deg, prec):
            pts.add(+mp.mpf(x))
    mp.mp.dps = old
    return sorted(pts)


def quad_full(f, a, b, level, dps):
    old = mp.mp.dps
    mp.mp.dps = dps + 20
    ts = mp.calculus.quadrature.TanhSinh(mp.mp)
    prec = mp.mp.prec
    results = []
    for deg in range(1, level + 1):
        nodes = ts.get_nodes(mp.mpf(a), mp.mpf(b), deg, prec)
        results.append(ts.sum_next(f, nodes, deg, prec, results, False))
    mp.mp.dps = old
    return results[-1]


# =====================================================================
# Closed-form subsector masters of the box family, EXACT in eps.
# box1eq family (lbl3vp-data.json): D1 = l^2 - w', D2 = (l+p1)^2 - 1,
# D3 = (l+p1+p2)^2 - 1, D4 = (l-p4)^2 - 1;  s = -1, t = -1/3, u = 4/3.
# Tadpoles are Gamma closed forms; bubbles/triangles are elementary
# one-folds (one Feynman integration done analytically, exact in eps),
# with a log substitution removing the 1/w' boundary layer at large w'.
# Verified 39-53 d against archived 1-loop AMFlow values at w'=5
# (<archive>/lbl3vp/r2-eichler/v2_kernel.py); the shipped box1eq_bnd_w5 block
# is used below ONLY as a live cross-check.
# =====================================================================
_GL = {}


def _gl(n):
    key = (n, mp.mp.dps)
    if key not in _GL:
        old = mp.mp.dps
        mp.mp.dps = old + 10
        _GL[key] = mp.gauss_quadrature(n, 'legendre')
        mp.mp.dps = old
    return _GL[key]


def glint(f, a, b, n):
    xs, ws = _gl(n)
    h = (b - a) / 2
    c = (b + a) / 2
    return h * mp.fsum(w * f(c + h * x) for x, w in zip(xs, ws))


def _inner(A, C, X, e):
    """int_0^X (A + C t)^{-1-eps} dt, stable for small eps and small C."""
    z = C * X / A
    if abs(z) < mp.mpf(10) ** (-mp.mp.dps // 3):
        return X * A ** (-1 - e) * (1 - (1 + e) * z / 2 + (1 + e) * (2 + e) * z * z / 6)
    return -A ** (-e) * mp.expm1(-e * mp.log1p(z)) / (e * C)


def _logpanels(L, n):
    """Panels for GL-n integration of an ~e^z integrand on [0, L]: panel
    length ~ 0.3*n keeps the per-panel exponential within spectral reach."""
    plen = max(2, n * 3 // 10)
    K = max(1, int(mp.ceil(L / plen)))
    return [(L * i / K, L * (i + 1) / K) for i in range(K)]


U_RAT = (4, 3)


def tad(m2, e):
    return mp.gamma(e) / (1 - e) * m2 ** (1 - e)


def bub_u(e, n):
    u = mp.mpf(U_RAT[0]) / U_RAT[1]
    f = lambda x: (1 - u * x * (1 - x)) ** (-e)
    return mp.gamma(e) * glint(f, mp.mpf(0), mp.mpf(1), n)


def tri_u1(e, n):
    u = mp.mpf(U_RAT[0]) / U_RAT[1]
    f = lambda x2: _inner(mp.mpf(1), -u * x2, 1 - x2, e)
    return -mp.gamma(1 + e) * glint(f, mp.mpf(0), mp.mpf(1), n)


def bub_s(wp, e, n):
    s = kin()[0]
    wm = wp - 1
    if abs(wm) < mp.mpf('0.5'):
        f = lambda x: (1 + wm * x - s * x * (1 - x)) ** (-e)
        return mp.gamma(e) * glint(f, mp.mpf(0), mp.mpf(1), n)
    lw = mp.log(wp)

    def g(z):
        ez = mp.e ** z
        x = (ez - 1) / wm
        return (ez - s * x * (1 - x)) ** (-e) * ez / wm
    return mp.gamma(e) * mp.fsum(glint(g, a, b, n) for a, b in _logpanels(lw, n))


def tri_s(wp, e, n):
    s = kin()[0]
    wm = wp - 1
    if abs(wm) < mp.mpf('0.5'):
        f = lambda x1: _inner(1 + wm * x1, -s * x1, 1 - x1, e)
        return -mp.gamma(1 + e) * glint(f, mp.mpf(0), mp.mpf(1), n)
    lw = mp.log(wp)

    def g(z):
        ez = mp.e ** z
        x1 = (ez - 1) / wm
        return _inner(ez, -s * x1, 1 - x1, e) * ez / wm
    return -mp.gamma(1 + e) * mp.fsum(glint(g, a, b, n) for a, b in _logpanels(lw, n))


def tri_u(wp, e, n):
    u = mp.mpf(U_RAT[0]) / U_RAT[1]
    wm = wp - 1
    if abs(wm) < mp.mpf('0.5'):
        f = lambda x2: _inner(wp + (1 - wp) * x2, (1 - wp) - u * x2, 1 - x2, e)
        return -mp.gamma(1 + e) * glint(f, mp.mpf(0), mp.mpf(1), n)
    lw = mp.log(wp)

    def g(z):
        ez = mp.e ** z
        y = (ez - 1) / wm
        return _inner(ez, -wm - u * (1 - y), y, e) * ez / wm
    # complex singularity of the endpoint factor 1-u*y*(1-y) sits at
    # |Im z| ~ 0.95, so tri_u needs short panels regardless of n
    K = max(1, int(mp.ceil(lw / mp.mpf('2.5'))))
    return -mp.gamma(1 + e) * mp.fsum(
        glint(g, lw * i / K, lw * (i + 1) / K, n) for i in range(K))


# =====================================================================
# Kernel: seedless VoP from the vacuum boundary w' -> infinity.
# =====================================================================
class Kernel:
    def __init__(self, eps_rat, DPS, N_TAY, SF):
        self.er = sp.Rational(eps_rat)
        self.e = mp.mpf(self.er.p) / self.er.q
        self.DPS = DPS
        self.N = N_TAY
        self.SF = SF
        self.n_inner = DPS // 2 + 14
        self.w0 = mp.mpf(200)
        W, D = sp.symbols('w d')
        dv = sp.Rational(4) - 2 * self.er
        self.A0f = [sp.lambdify(W, sp.cancel(sp.sympify(
            DATA['box1eq_de']['A'][0][j], locals={'w': W, 'd': D}).subs(D, dv)),
            'mpmath') for j in range(8)]
        b_sings = [0, 1, mp.mpc(0, 1), mp.mpc(0, -1), mp.mpc(0, 2), mp.mpc(0, -2),
                   mp.mpc(mp.mpf(1) / 3, mp.sqrt(mp.mpf(8)) / 3),
                   mp.mpc(mp.mpf(1) / 3, -mp.sqrt(mp.mpf(8)) / 3)]
        self.de = RatDE(DATA['box1eq_de'], self.er, b_sings)
        e, n = self.e, self.n_inner
        self.cst = {'tad1': tad(mp.mpf(1), e), 'bubu': bub_u(e, n),
                    'triu1': tri_u1(e, n)}
        self.ap = (1 + 2 * e) / 2

    def logH(self, v):
        return -self.ap * mp.log(1 + v * v)

    def masters(self, v, n=None):
        e = self.e
        n = n or self.n_inner
        return [None, tad(v, e), self.cst['tad1'], bub_s(v, e, n),
                self.cst['bubu'], tri_s(v, e, n), tri_u(v, e, n),
                self.cst['triu1']]

    def R(self, v, n=None):
        M = self.masters(v, n)
        return mp.fsum(self.A0f[j](v) * M[j] for j in range(1, 8))

    def K_w0(self):
        """K(200) by the VoP one-fold from the w'->inf vacuum boundary,
        seed K(W) = -Tri_u1(eps)/W, W = 200*10^dps (relative seed error
        O(1/W) sits below the working precision). The inner one-fold order
        shrinks with v: the [v, e*v] band contributes O(200/v) of K(200),
        so far bands need proportionally fewer digits. Used as the LIVE
        independent cross-check of the parametric anchor (its quadrature
        carries a small ~1/eps-scaling residual, ~1e-30 at the defaults)."""
        e = self.e
        old = mp.mp.dps
        exp10 = old
        mp.mp.dps = old + int(2 * float(e) * (exp10 + 3)) + 10
        try:
            Winf = self.w0 * mp.mpf(10) ** exp10
            ngl0 = int(os.environ.get('VP_NGL_INF', 0)) or max(16, self.DPS // 5 + 6)
            lw, lW = mp.log(self.w0), mp.log(Winf)
            nseg = int(mp.ceil((lW - lw) / mp.log(10) * 2))  # 2 segments/decade
            seg = (lW - lw) / nseg
            acc = -self.cst['triu1'] / Winf * mp.e ** (-self.logH(Winf))
            for i in range(nseg):
                xa = lW - i * seg
                dec_above = float((xa - lw) / mp.log(10))   # decades above w0
                nloc = max(10, self.n_inner - int(3.0 * dec_above))
                ngl = max(10, ngl0 - int(dec_above) // 3)
                f = lambda x: (lambda v: self.R(v, nloc) * mp.e ** (-self.logH(v))
                               * v)(mp.e ** x)
                acc += glint(f, xa, xa - seg, ngl)
            K0 = mp.e ** self.logH(self.w0) * acc
        finally:
            mp.mp.dps = old
        return +K0

    def K_anchor(self):
        """K(w0) from the exact 2-fold Feynman-parametric representation
        (one inner integration done analytically, exact in eps):
          K = Gamma(2+eps) Int dx1 dx2 [A^{-1-eps} - (A+BX)^{-1-eps}]/((1+eps)B),
          A = w0 x1 + (1-x1) - u x2(1-x1-x2), B = -s x1 + u x2, X = 1-x1-x2.
        Elementary integrand, arbitrary precision (tanh-sinh)."""
        e = self.e
        s = kin()[0]
        u = mp.mpf(U_RAT[0]) / U_RAT[1]
        W = self.w0
        old = mp.mp.dps
        mp.mp.dps = old + 10
        try:
            def f(x1, x2):
                A = W * x1 + (1 - x1) - u * x2 * (1 - x1 - x2)
                B = -s * x1 + u * x2
                X = 1 - x1 - x2
                return (A ** (-1 - e) - (A + B * X) ** (-1 - e)) / ((1 + e) * B)
            K0 = mp.gamma(2 + e) * mp.quad(
                lambda x1: mp.quad(lambda x2: f(x1, x2), [0, 1 - x1]), [0, 1])
        finally:
            mp.mp.dps = old
        return +K0

    def transport(self, tg_desc, tg_asc):
        """Taylor-step the 8-vector, seeded ANALYTICALLY at w0=200 by the
        parametric anchor, to all real-axis targets. Returns {w': K}."""
        V0 = [self.K_anchor()] + self.masters(self.w0)[1:]
        K_at = {}
        for tg in (tg_desc, tg_asc):
            if tg:
                res, _ = self.de.path(V0, self.w0, tg, self.N, self.SF)
                for v in tg:
                    K_at[v] = res[v][0]
        self._V0 = V0
        return K_at

    def circle_cn(self, K_start, r, nmax):
        """Threshold Taylor coefficients c_n of K at w'=1: arc-VoP around the
        Cauchy circle |w'-1| = r starting from K(1+r), trapezoid/FFT sum.
        Returns (coeffs c_0..c_nmax, closure defect of the loop)."""
        M = int(os.environ.get('VP_CIRCLE_M', 0)) or max(128, int(2.4 * self.DPS))
        acc = K_start * mp.e ** (-self.logH(1 + r))
        vals = [K_start]
        th = [2 * mp.pi * k / M for k in range(M + 1)]
        I = mp.mpc(0, 1)
        pt = lambda t: 1 + r * mp.e ** (I * t)
        for k in range(1, M + 1):
            a, b = th[k - 1], th[k]
            f = lambda t: (lambda v: self.R(v) * mp.e ** (-self.logH(v))
                           * r * I * mp.e ** (I * t))(pt(t))
            acc += glint(f, a, b, 8)
            vals.append(mp.e ** self.logH(pt(th[k])) * acc)
        closure = abs(vals[M] - vals[0]) / abs(vals[0])
        cs = []
        for nn in range(nmax + 1):
            ssum = mp.fsum(vals[k] * mp.e ** (-I * nn * th[k]) for k in range(M)) / M
            cs.append(mp.re(ssum / r ** nn))
        return cs, closure


# =====================================================================
# Density: closed-form 2-body + Gamma_1(6) period one-fold 3-body.
# =====================================================================
def rho_2b(w, e):
    G = mp.gamma(1 + e) * mp.gamma(1 - e) / (e * mp.gamma(2 - 2 * e))
    return -mp.pi * G * (-(1 - 2 * e) * (w + 1) / (w - 1) + e * (w - 1) / 6) \
        * (w - 1) ** (-2 * e) * w ** (e - 1)


def rho_3b(w, e, qdps):
    """Exact-in-d three-body density: sunrise-elliptic period one-fold (see
    module docstring). tanh-sinh with log-split interior points (handles the
    (q-4)^{1/2-eps}, (hi-q)^{1/2-eps} endpoint branches at any w)."""
    old = mp.mp.dps
    mp.mp.dps = qdps + 15
    try:
        w = mp.mpf(w)
        sw = mp.sqrt(w)
        hi = (sw - 1) ** 2
        if hi <= 4:
            return mp.mpf(0)
        hb = (sw + 1) ** 2
        c1 = mp.gamma(1 - e) / mp.gamma(2 - 2 * e)
        he = mp.mpf(1) / 2 - e

        def f(q):
            return (q - 4) ** he * q ** mp.mpf('-2.5') * (hi - q) ** he * (hb - q) ** he
        pts = [mp.mpf(4)]
        c = mp.mpf(40)
        while c < hi / 4:
            pts.append(c)
            c *= 100
        pts.append(hi)
        try:
            I = mp.quad(f, pts)
        except ZeroDivisionError:
            # mpmath TanhSinh.estimate_error divides by an exactly-zero level
            # difference when two successive levels agree to the bit (far-field
            # nodes at very large w'; first seen at dps 100, 2026-07-05).
            # Retry off the tie at +7 dps, else fixed-level panel sums (no
            # error estimator). Exception path only: non-degenerate runs are
            # bit-identical to the previous behavior.
            try:
                mp.mp.dps += 7
                I = mp.quad(f, pts)
            except ZeroDivisionError:
                I = mp.fsum(quad_full(f, a, b, 6, mp.mp.dps)
                            for a, b in zip(pts[:-1], pts[1:]))
            finally:
                mp.mp.dps -= 7
        out = mp.pi * c1 * c1 * w ** (e - 1) * I
    finally:
        mp.mp.dps = old
    return +out


# =====================================================================
# Per-eps assembly (round-1 panel algebra; new kernel + density)
# =====================================================================
def I_disp_at_eps(eps_rat, DPS, N_TAY, SF, LMAX, work_dps=None, verbose=True):
    """I_VP(eps) at s=-1, t=-1/3, m^2=1, any rational 0 < eps <= 1/16.
    Fully analytic: no numeric seeds; everything below RUNS now."""
    t_e = time.time()
    if work_dps is None:
        work_dps = max(45, DPS - 25)
    eps_rat = sp.Rational(eps_rat)
    mp.mp.dps = DPS
    e = mp.mpf(eps_rat.p) / eps_rat.q
    G = mp.gamma(1 + e) * mp.gamma(1 - e) / (e * mp.gamma(2 - 2 * e))

    # kernel pipeline precision: subsector masters carry Gamma(eps) ~ 1/eps
    # while K is O(1); the assembly cancellation costs ~log10(1/eps) digits.
    kdps = DPS + int(mp.ceil(-mp.log10(mp.mpf(eps_rat.p) / eps_rat.q))) + 2
    mp.mp.dps = kdps
    ker = Kernel(eps_rat, DPS, N_TAY, SF)
    mp.mp.dps = DPS

    # quadrature nodes / transport targets (levels LMAX and LMAX+1)
    mp.mp.dps = work_dps + 20
    LN = LMAX + 1
    bb = mp.mpf('1.4')
    Lb = bb - 1
    WP_MAX = mp.mpf(10) ** max(35, work_dps + 5)
    vs = vs_for(DPS)
    rc = mp.mpf('0.375')            # dyadic: exact at every precision
    n_p1a = ts_nodes(mp.mpf(0), Lb, LN, work_dps)
    n_p1b = ts_nodes(bb, mp.mpf(9), LN, work_dps)
    n_p2 = ts_nodes(mp.mpf(9), mp.mpf(200), LN, work_dps)
    n_tt = ts_nodes(mp.mpf(-1), mp.mpf(1), LN, work_dps)
    n_tail = [mp.mpf(200) + mp.mpf(200) * (1 + t) / (1 - t) for t in n_tt if t < 1]
    tg_desc = sorted(set([mp.mpf(1) + rc]
                         + [v for v in (n_p1b + n_p2) if v < 200]
                         + [mp.mpf(1) + v for v in n_p1a if v >= vs]
                         + ts_nodes(mp.mpf(1) + vs, bb, LN, work_dps)), reverse=True)
    tg_asc = sorted([v for v in (n_p2 + n_tail) if 200 <= v < WP_MAX])
    mp.mp.dps = DPS

    # === kernel: seedless transport + threshold circle (at kdps) ===
    t0 = time.time()
    mp.mp.dps = kdps
    K_at = ker.transport(tg_desc, tg_asc)
    tK = time.time() - t0
    t0 = time.time()
    nmax = int(mp.ceil(DPS / -mp.log10(vs))) + 8
    cs, closure = ker.circle_cn(K_at[mp.mpf(1) + rc], rc, nmax)
    mp.mp.dps = DPS
    tC = time.time() - t0
    K0e, K1e, K2e = cs[0], cs[1], cs[2]

    # live cross-checks (archived 1-loop AMFlow data; NEVER used as input)
    box_b = load_laur(DATA['box1eq_bnd_w5'])
    Knj = {int(n): {int(o): mp.mpf(v) for o, v in d.items()}
           for n, d in DATA['Kn'].items()}
    K5 = None
    if mp.mpf(5) in K_at:
        K5 = K_at[mp.mpf(5)]
    else:
        res5, _ = ker.de.path(ker._V0, ker.w0, [mp.mpf(5)], N_TAY, SF)
        K5 = res5[mp.mpf(5)][0]
    chk_K5 = agree_d(mp.re(K5), mp.re(laur_at(box_b[(1, 1, 1, 1)], e, -1, 20)))
    chk_K2 = agree_d(K2e, sum(Knj[2].get(j, mp.mpf(0)) * e ** j for j in range(21)))
    # independent SEEDLESS check: VoP one-fold from the w'->inf vacuum
    # boundary vs the parametric anchor (skip with VP_VAC_CHECK=0)
    chk_vac = None
    if os.environ.get('VP_VAC_CHECK', '1') != '0':
        mp.mp.dps = kdps
        chk_vac = agree_d(ker.K_w0(), ker._V0[0])
        mp.mp.dps = DPS
    diag = {'nK': len(K_at), 'tK': tK, 'tC': tC, 'closure': float(closure),
            'chk_K5': chk_K5, 'chk_K2': chk_K2, 'chk_vac': chk_vac}
    if verbose:
        print(f"    kernel: transport {len(K_at)} pts {tK:.0f}s | circle {tC:.0f}s"
              f" (closure {closure:.1e})")
        print(f"    [check] K(5) vs archived 1-loop AMFlow: {chk_K5:.1f} d |"
              f" K_2(eps) circle vs archived series: {chk_K2:.1f} d"
              + (f" | vacuum-boundary VoP vs anchor: {diag['chk_vac']:.1f} d"
                 if diag['chk_vac'] is not None else ""))

    K_keys = sorted(K_at.keys())
    miss = []

    def _near(keys, wp):
        i = bisect.bisect_left(keys, wp)
        for j in (i - 1, i, i + 1):
            if 0 <= j < len(keys) and abs(keys[j] - wp) < abs(wp) * mp.mpf('1e-30'):
                return keys[j]
        return None

    def K_P(wp):
        wp = mp.mpf(wp)
        u = wp - 1
        if abs(u) < vs:
            acc = mp.mpf(0)
            for c in reversed(cs[2:]):
                acc = acc * u + c
            return acc
        if wp >= WP_MAX:
            return mp.mpf(0)
        Kv = K_at.get(wp)
        if Kv is None:
            kk = _near(K_keys, wp)
            if kk is not None:
                Kv = K_at[kk]
        if Kv is None:
            # far-field: K_P ~ -K1/w' ~ 1e-20 there; a zeroed stray node costs
            # ~1e-40 relative (counted and reported in the diagnostics)
            miss.append(wp)
            if wp > mp.mpf('1e20'):
                return mp.mpf(0)
            raise KeyError(mp.nstr(wp, 20))
        return (mp.re(Kv) - K0e - K1e * u) / (u * u)

    # === density (closed forms; 3-body one-fold memoized per node) ===
    rho3_cache = {}

    def rho(wp):
        wp = mp.mpf(wp)
        if wp <= 1 or wp >= WP_MAX:
            return mp.mpf(0)
        r = rho_2b(wp, e)
        if wp > 9:
            key = mp.nstr(wp, work_dps)
            r3 = rho3_cache.get(key)
            if r3 is None:
                # far nodes carry ~(200/w') weight in the dispersion integral
                qd = max(20, work_dps + 10 - (int(mp.log10(wp / 200)) if wp > 200 else 0))
                r3 = rho_3b(wp, e, qd)
                rho3_cache[key] = r3
            r += r3
        return r

    # live check: closed-form density vs nothing transported — instead check
    # the 3-body fold at doubled quadrature precision at one point
    w_chk = mp.mpf(12)
    chk_r3 = agree_d(rho_3b(w_chk, e, work_dps + 10), rho_3b(w_chk, e, 2 * work_dps + 10))
    diag['chk_r3'] = chk_r3
    if verbose:
        print(f"    [check] rho_3b(12) quadrature doubling: {chk_r3:.1f} d")

    # === threshold panel + main panels (round-1 algebra, unchanged) ===
    def f1(u):
        u = mp.mpf(u)
        return (2 + u) / (1 + u) ** (1 - e) * K_P(1 + u)

    def f2(u):
        u = mp.mpf(u)
        return K_P(1 + u) / (1 + u) ** (1 - e)

    f10 = 2 * K2e

    def assem(level):
        old = mp.mp.dps
        mp.mp.dps = work_dps + 20
        A = f10 * Lb ** (-2 * e) / (-2 * e)
        A += quad_full(lambda u: (mp.mpf(u)) ** (-2 * e) * (f1(u) - f10) / mp.mpf(u),
                       mp.mpf(0), Lb, level, work_dps)
        A *= G * (1 - 2 * e)
        B = -G * e / 6 * quad_full(lambda u: mp.mpf(u) ** (1 - 2 * e) * f2(u),
                                   mp.mpf(0), Lb, level, work_dps)
        P1 = A + B
        P2 = quad_full(lambda wp: rho(wp) * K_P(wp) / mp.pi, bb, mp.mpf(9),
                       level, work_dps)
        P3 = quad_full(lambda wp: rho(wp) * K_P(wp) / mp.pi, mp.mpf(9),
                       mp.mpf(200), level, work_dps)

        def ft(t):
            wp = mp.mpf(200) + mp.mpf(200) * (1 + t) / (1 - t)
            return rho(wp) * K_P(wp) / mp.pi * mp.mpf(200) * 2 / (1 - t) ** 2

        P4 = quad_full(ft, mp.mpf(-1), mp.mpf(1), level, work_dps)
        mp.mp.dps = old
        return P1 + P2 + P3 + P4

    I5 = assem(LMAX)
    I6 = assem(LMAX + 1)
    R56 = 2 * I6 - I5
    diag['quad_agree'] = agree_d(I5, I6)
    diag['nmiss'] = len(miss)
    diag['t_total'] = time.time() - t_e
    if verbose:
        print(f"    quad levels L{LMAX}/L{LMAX+1} agree {diag['quad_agree']:.1f} d"
              f" -> Richardson  ({diag['t_total']:.0f}s total)")
    return R56, diag


# =====================================================================
# Laurent pole layers eps^-2, eps^-1 (round-1 machinery): closed-form dilog
# kernel; K_2^{(1)} runtime-derived (K2_threshold_eps01) since 2026-07-05 —
# no shipped numeric inputs; depth knobs scale with --dps (SPEC_ROW33).
# =====================================================================
_GL_CACHE = {}


def _gl_work(n, work):
    key = (n, work)
    xw = _GL_CACHE.get(key)
    if xw is None:
        old = mp.mp.dps
        mp.mp.dps = work
        xw = mp.gauss_quadrature(n, 'legendre')
        mp.mp.dps = old
        _GL_CACHE[key] = xw
    return xw


def _g_of_v(v, s, u, m2, M5, a):
    vp1 = v + 1
    p1 = a * v + (M5 + m2)
    r1 = m2 * vp1 * vp1
    p2 = (M5 + m2) * vp1
    r2 = r1 - u * v
    sd1 = mp.sqrt(p1 * p1 - 4 * M5 * r1)
    sd2 = mp.sqrt(p2 * p2 - 4 * M5 * r2)
    twoM5 = 2 * M5
    mal = (p1 - sd1) / twoM5
    malp = (p1 + sd1) / twoM5
    mbe = (p2 - sd2) / twoM5
    mbep = (p2 + sd2) / twoM5
    Ral = -1 / ((u + s * mal) * sd1)
    Ralp = 1 / ((u + s * malp) * sd1)
    Rbe = 1 / ((u + s * mbe) * sd2)
    Rbep = -1 / ((u + s * mbep) * sd2)
    return -(Ral * mp.log(mal) + Ralp * mp.log(malp)
             + Rbe * mp.log(mbe) + Rbep * mp.log(mbep))


def box_closed_dilog(M5, dps=40):
    work = dps + 18
    old_dps = mp.mp.dps
    mp.mp.dps = work
    s, t, m2 = kin()
    M5 = mp.mpf(M5)
    u = -s - t
    a = M5 + m2 - s
    n = work + 12
    nodes, weights = _gl_work(n, work)
    L = mp.mpf(2)
    one = mp.mpf(1)
    tot = mp.mpf(0)
    for tk, wk in zip(nodes, weights):
        v = L * (one + tk) / (one - tk)
        dv = L * 2 / (one - tk) ** 2
        tot += wk * _g_of_v(v, s, u, m2, M5, a) * dv
    mp.mp.dps = dps
    return +tot


# Runtime derivation of K_2^{(0)}, K_2^{(1)} — retires the shipped AMFlow
# string Kn['2']['1'] (2026-07-05 wiring; the LBL3VP build spec for the derivation +
# verification: two-precision 161.4 d, vs stored 120.6 d = full stored length).
def K2_threshold_eps01(dps, guard=25):
    r"""eps^0 and eps^1 of the 2nd threshold Taylor coefficient K_2(eps) of the
    box kernel at w'=1, computed from the exact parametric 2-fold (no seeds).

    K analytic at w'=1 (A >= 2/3 on the simplex); differentiating the exact
    representation K = Gamma(2+eps) II [A^{-1-eps}-C^{-1-eps}]/((1+eps)B)
    twice in w' (dA/dw' = dC/dw' = x1, C = A + BX) gives, at w'=1,

        K_2(eps) = [Gamma(2+eps)(2+eps)/2] II x1^2 [A^{-3-eps}-C^{-3-eps}]/B
        K_2^{(0)} = Q0,   K_2^{(1)} = (3-2*euler)/2 * Q0 + Q1

    with the cancellation-free divided-difference integrands (C-A = BX):
        Q0: x1^2 X (A^2+AC+C^2)/(A^3 C^3)
        Q1: x1^2 X [A^{-3} log1p(BX/A)/(BX) - (A^2+AC+C^2)/(A^3 C^3) ln C]
    Tensor Gauss-Legendre on (x1,x2)=(xi,eta(1-xi)); integrand analytic on the
    closed square -> geometric convergence (n = 2.2*dps + 40 nodes).
    """
    work = dps + guard
    n = int(mp.mpf('2.2') * dps) + 40
    old = mp.mp.dps
    mp.mp.dps = work
    try:
        s, t, m2 = kin()
        u = -s - t
        nodes, weights = _gl_work(n, work)
        one = mp.mpf(1)
        half = one / 2
        xs = [half * (one + tk) for tk in nodes]
        ws = [half * wk for wk in weights]
        Q0 = mp.mpf(0)
        Q1 = mp.mpf(0)
        for xi, wxi in zip(xs, ws):            # xi = x1
            om = one - xi
            x1sq_j = xi * xi * om * wxi        # x1^2 * Jacobian * weight
            acc0 = mp.mpf(0)
            acc1 = mp.mpf(0)
            for et, wet in zip(xs, ws):
                x2 = et * om
                X = om - x2
                A = one - u * x2 * X
                B = -s * xi + u * x2
                BX = B * X
                C = A + BX
                A3 = A ** -3
                R = (A * A + A * C + C * C) * A3 * C ** -3
                ddlog = mp.log1p(BX / A) / BX if BX != 0 else 1 / A
                acc0 += wet * X * R
                acc1 += wet * X * (A3 * ddlog - R * mp.log(C))
            Q0 += x1sq_j * acc0
            Q1 += x1sq_j * acc1
        K2_0 = Q0
        K2_1 = (3 - 2 * mp.euler) / 2 * Q0 + Q1
        mp.mp.dps = dps
        K2_0 = +K2_0            # round at dps BEFORE restoring caller precision
        K2_1 = +K2_1            # (module-level dps=15 would truncate otherwise)
    finally:
        mp.mp.dps = old
    return K2_0, K2_1


# Depth-knob scaling for the pole-layer demo (2026-07-05, SPEC_ROW33 measured
# rates): at --dps <= LEGACY_DPS every depth knob freezes to the shipped value,
# so the default invocation's output is unchanged; above it the fit degree
# scales so the delivered layer digits track --dps instead of the frozen
# nfit=28 truncation cap (measured 45.03 d at any dps).
LEGACY_DPS = 60


def nfit_for(DPS):
    """Chebyshev-fit degree for the pole layers.  Measured rates (SPEC_ROW33):
    the r=0.05 fit caps c2 at 1.61 d/order (nfit 28/44/60 -> 45.03/70.83/96.58 d,
    flat in dps); the |u|<0.04 Horner window carries log10(1/0.04)=1.398 d/term,
    the binding (slower) rate.  Target: DPS+8 digits at the window rate."""
    if DPS <= LEGACY_DPS:
        return 28
    return max(28, math.ceil((DPS + 8) / 1.398) + 1)


def kernel_taylor_live(nfit=28, r='0.05', dps=100):
    """Taylor coefficients c_n of K(1+u)|_{eps^0} at u=0 from the CLOSED-FORM
    dilog kernel: Chebyshev-Vandermonde fit (round-1 route, kept for the
    pole-layer demo; the fixed-eps gate uses the arc-VoP circle instead).
    The Vandermonde solve guard grows with nfit (conditioning headroom) for
    nfit > 28; at the legacy depth it is the shipped +40."""
    sg = 40 + (0 if nfit <= 28 else nfit // 2)
    old = mp.mp.dps
    mp.mp.dps = dps + sg
    r = mp.mpf(r)
    us = [r * mp.cos(mp.pi * j / nfit) for j in range(nfit + 1)]
    vals = mp.matrix([[box_closed_dilog(1 + u, dps)] for u in us])
    mp.mp.dps = dps + sg
    V = mp.matrix(nfit + 1, nfit + 1)
    for i, u in enumerate(us):
        p = mp.mpf(1)
        for j in range(nfit + 1):
            V[i, j] = p
            p *= u
    c = mp.lu_solve(V, vals)
    mp.mp.dps = old
    return [c[j] for j in range(nfit + 1)]


def pole_layers(Knj, qdps=55, kdps=90, nfit=28):
    """eps^-2 and eps^-1 Laurent coefficients of I_VP, live (round-1 route).
    Seed-free since 2026-07-05: K_2^{(1)} is computed at runtime by
    K2_threshold_eps01; the stored AMFlow string is a held-out gate only.
    Depth knobs (qdps, kdps, nfit) are wired to --dps by main() (frozen at the
    shipped values for dps <= LEGACY_DPS)."""
    mp.mp.dps = kdps + 40
    c = kernel_taylor_live(nfit=nfit, r='0.05', dps=kdps)
    K2_live = c[2]
    K2_ship = Knj[2][0]
    fit_chk = agree_d(K2_live, K2_ship)

    K0, K1 = c[0], c[1]

    def K_P0(wp):
        wp = mp.mpf(wp)
        u = wp - 1
        if abs(u) < mp.mpf('0.04'):
            acc = mp.mpf(0)
            for cc in reversed(c[2:]):
                acc = acc * u + cc
            return acc
        old = mp.mp.dps
        v = (box_closed_dilog(wp, kdps - 20) - K0 - K1 * u) / (u * u)
        mp.mp.dps = old
        return +v

    def f1(wp):
        wp = mp.mpf(wp)
        u = wp - 1
        if abs(u) < mp.mpf('0.04'):
            acc = mp.mpf(0)
            for cc in reversed(c[3:]):
                acc = acc * u + cc
            return acc
        return (K_P0(wp) - K2_live) / u

    def f2(wp):
        return K_P0(wp) / (mp.mpf(wp) - 1)

    def f3(wp):
        return K_P0(wp) / mp.mpf(wp)

    b = mp.mpf('1.4')
    mp.mp.dps = qdps
    I1 = mp.quad(f1, [1, b], method='tanh-sinh')
    I2 = mp.quad(f2, [b, 9], method='tanh-sinh') + mp.quad(f2, [9, 200], method='tanh-sinh')
    I2 += mp.quad(lambda t: f2(200 + 200 * (1 + t) / (1 - t)) * 200 * 2 / (1 - t) ** 2,
                  [-1, 1], method='tanh-sinh')
    I3 = mp.quad(f3, [1, 9], method='tanh-sinh') + mp.quad(f3, [9, 200], method='tanh-sinh')
    I3 += mp.quad(lambda t: f3(200 + 200 * (1 + t) / (1 - t)) * 200 * 2 / (1 - t) ** 2,
                  [-1, 1], method='tanh-sinh')
    Bsm1 = 2 * K2_live * (mp.euler + mp.log(b - 1)) + 2 * I1 + 2 * I2 - I3
    K2_0_live2, K2_1_live = K2_threshold_eps01(kdps)   # runtime, no seeds
    chk_K21 = agree_d(K2_1_live, Knj[2][1])   # stored AMFlow string -> held-out gate
    K2_1 = K2_1_live
    I_em2 = -K2_live
    I_em1 = -K2_1 + Bsm1
    return I_em2, I_em1, fit_chk, chk_K21


# =====================================================================
# CLI / gate runner
# =====================================================================
def parse_eps(txt):
    txt = txt.strip()
    if txt.startswith('2^-'):
        return sp.Rational(1, 2 ** int(txt[3:]))
    return sp.Rational(txt)


def eps_label(e):
    if e.p == 1 and (e.q & (e.q - 1)) == 0:
        return f"2^-{e.q.bit_length() - 1}"
    return f"{e.p}/{e.q}"


def oracle_for(e):
    if e.p == 1 and (e.q & (e.q - 1)) == 0:
        k = e.q.bit_length() - 1
        return ORC['I_AMF_grid'].get(str(k))
    return None


def run_gate(pts, DPS, N_TAY, SF, LMAX, n_jobs):
    def run_pt(ee):
        val, diag = I_disp_at_eps(ee, DPS, N_TAY, SF, LMAX, verbose=False)
        return mp.nstr(val, max(60, DPS)), diag

    results = {}
    if n_jobs > 1 and len(pts) > 1:
        import multiprocessing
        ctx = multiprocessing.get_context('fork')

        def worker(i, conn):
            try:
                conn.send((i,) + run_pt(pts[i]))
            except Exception as ex:          # pragma: no cover
                conn.send((i, None, {'error': repr(ex)}))
            conn.close()

        pipes, procs = {}, {}
        for i in range(len(pts)):
            pr, pw = ctx.Pipe(False)
            p = ctx.Process(target=worker, args=(i, pw))
            p.start()
            pipes[i], procs[i] = pr, p
        for i in range(len(pts)):
            ii, vstr, diag = pipes[i].recv()
            procs[i].join()
            results[ii] = (vstr, diag)
    else:
        for i, ee in enumerate(pts):
            results[i] = run_pt(ee)

    mind = mp.inf
    bank = {}
    for i, ee in enumerate(pts):
        vstr, diag = results[i]
        if vstr is None:
            print(f"  eps = {eps_label(ee)}: FAILED: {diag.get('error')}")
            continue
        mp.mp.dps = max(80, DPS + 20)
        val = mp.mpf(vstr)
        print(f"  eps = {eps_label(ee)}:")
        print(f"    kernel: transport {diag['nK']} pts {diag['tK']:.0f}s |"
              f" circle {diag['tC']:.0f}s (closure {diag['closure']:.1e})")
        print(f"    [check] K(5) vs archived 1-loop AMFlow {diag['chk_K5']:.1f} d |"
              f" K_2(eps) circle vs archived series {diag['chk_K2']:.1f} d |"
              f" rho_3b dps-doubling {diag['chk_r3']:.1f} d"
              + (f" | vacuum VoP vs anchor {diag['chk_vac']:.1f} d"
                 if diag.get('chk_vac') is not None else ""))
        print(f"    quad levels L{LMAX}/L{LMAX+1} agree {diag['quad_agree']:.1f} d"
              f" -> Richardson  ({diag['t_total']:.0f}s"
              + (f"; {diag['nmiss']} far-field node(s) zeroed" if diag.get('nmiss') else "")
              + ")")
        print(f"    I_VP    = {mp.nstr(val, 40)}   (computed now)")
        ostr = oracle_for(ee)
        if ostr is None:
            print("    no held-out oracle at this eps (grid: eps=2^-k, k=4..18);"
                  " value computed, ungated\n")
            continue
        orac = mp.mpf(ostr)
        d = agree_d(val, orac)
        mind = min(mind, d)
        print(f"    I_AMF   = {mp.nstr(orac, 40)}   (held-out 3-loop oracle)")
        print(f"    agreement = {d:.2f} d\n")
        bank[eps_label(ee)] = {
            'I_VP': vstr, 'banked_d': float(d),
            'closure': float(diag['closure']),
            'quad_agree': float(diag['quad_agree']),
            'chk_K5': float(diag['chk_K5']), 'chk_K2': float(diag['chk_K2']),
            'chk_r3': float(diag['chk_r3'])}
    return mind, bank


def settings_for(DPS):
    N_TAY = int(os.environ.get('VP_TAYLOR_N', 0)) or max(80, int(0.95 * DPS))
    LMAX = int(os.environ.get('VP_LMAX', 0)) or (5 if DPS <= 90 else 6 if DPS <= 200 else 7)
    return N_TAY, LMAX


def main():
    import argparse
    ap = argparse.ArgumentParser(
        description="LBL3VP fully-analytic (seed-free) evaluation at s=-1, t=-1/3, "
                    "m^2=1; runtime inputs: eps (rational, 0 < eps <= 1/16), dps.")
    ap.add_argument('--point', metavar='EPS', default=None,
                    help="eps as '2^-k', 'p/q' or exact decimal "
                         "(default: demo grid 2^-6,2^-9,2^-13)")
    ap.add_argument('--dps', type=int, default=None,
                    help="working precision in decimal digits (default 70)")
    ap.add_argument('--double', action='store_true',
                    help="dps-doubling demo: rerun one gate point at dps and 2*dps")
    ap.add_argument('--no-fastcache', action='store_true',
                    help="skip the sha-pinned fast-start cache "
                         "(lbl3vp-fastcache.json) and run the live machinery "
                         "even at banked eps points; the cache is a "
                         "fast-start convenience ONLY -- the live machinery "
                         "in this file remains the definition")
    ap.add_argument('--fastcache-bank', action='store_true',
                    help="maintainer mode: run fully live and, on a gate run "
                         "meeting the >=30d bars, write/refresh "
                         "lbl3vp-fastcache.json (banked values + measured "
                         "agreements + sha pins)")
    args = ap.parse_args()

    t_start = time.time()
    DPS = args.dps if args.dps else int(os.environ.get('VP_DPS', 60))
    SF = mp.mpf(os.environ.get('VP_STEP_FRAC', '0.28'))
    N_TAY, LMAX = settings_for(DPS)

    if args.point is not None:
        pts = [parse_eps(args.point)]
    else:
        pts = [sp.Rational(1, 2 ** int(x))
               for x in os.environ.get('VP_EPS_KS', '6,9,13').split(',')]
    for ee in pts:
        if not (0 < ee <= sp.Rational(1, 16)):
            raise SystemExit(f"eps={ee} outside the supported domain 0 < eps <= 1/16"
                             " (see module docstring)")

    n_jobs = int(os.environ.get('VP_JOBS', min(len(pts), os.cpu_count() or 1)))

    _fc, _fc_note = _fc_eligible(args, DPS, pts)
    if _fc_note:
        print(f"[fast-cache] {_fc_note}")
    if _fc is not None:
        _fc_fast_path(_fc, args, DPS, pts, t_start)  # SystemExit(0) on success

    print("LBL3VP fully-analytic evaluation:  s=-1, t=-1/3, m^2=1  (kinematics fixed"
          " by the shipped exact rational connection; runtime inputs: eps, dps)")
    print("I_VP(eps) = (1/pi) int_1^inf rho_VP(w';eps) K_P(w';eps) dw'   with")
    print("  rho = closed-form 2-body (exact in d) + Gamma_1(6) sunrise-period"
          " one-fold 3-body (exact in d),")
    print("  K   = exact 2-fold parametric anchor + exact-rational-connection"
          " Taylor stepping (vacuum-boundary VoP fold cross-checked live).")
    print(f"settings: dps={DPS}, Taylor order={N_TAY}, quad levels {LMAX}/{LMAX+1},"
          f" {n_jobs} worker(s)\n")

    if args.double:
        ee = pts[0] if args.point is not None else sp.Rational(1, 1024)
        print(f"dps-doubling demo at eps = {eps_label(ee)}: same analytic pipeline,"
              f" working precision {DPS} then {2*DPS}\n")
        legs = []
        for dcur in (DPS, 2 * DPS):
            ncur, lcur = settings_for(dcur)
            print(f"-- pass at dps={dcur} (Taylor {ncur}, quad levels {lcur}/{lcur+1},"
                  f" quad work-dps {max(45, dcur - 25)}):")
            t0 = time.time()
            d, _ = run_gate([ee], dcur, ncur, SF, lcur, 1)
            legs.append((dcur, d, time.time() - t0))
        (d1c, d1, t1), (d2c, d2, t2) = legs
        print(f"dps {d1c} -> {d2c}: live agreement {d1:.2f} d -> {d2:.2f} d"
              f"  (walls {t1:.0f}s / {t2:.0f}s)")
        print("remaining caps: oracle grid strings ~320 d (AMFlow working_pre 340);"
              " all round-1 analytic-input caps (eps^20 seeds, shipped threshold"
              " series) are gone — inputs are exact formulas.")
        print(f"\ntotal wall time: {time.time()-t_start:.0f}s")
        return

    mind, _fc_pts = run_gate(pts, DPS, N_TAY, SF, LMAX, n_jobs)
    if mp.isfinite(mind):
        print(f"  min gate over {{{', '.join(eps_label(ee) for ee in pts)}}}:"
              f" {float(mind):.2f} d")
        print("  (round-1 archived dispersion gate at these settings: 31.70-32.79 d;"
              " that form is superseded by this seed-free representation)\n")

    if args.point is not None:
        print("(--point mode: Laurent pole-layer demo skipped; run without flags"
              " to see it)")
        if args.fastcache_bank:
            print("[fast-cache] bank REFUSED: banking requires the full "
                  "default demo (no --point)")
        print(f"\ntotal wall time: {time.time()-t_start:.0f}s")
        return

    print("Laurent pole layers (computed now from the closed-form dilog kernel):")
    mp.mp.dps = 140
    Knj = {int(n): {int(o): mp.mpf(v) for o, v in d.items()}
           for n, d in DATA['Kn'].items()}
    I_em2, I_em1, fit_chk, chk_K21 = pole_layers(Knj, qdps=max(55, DPS - 15),
                                                 kdps=max(90, DPS + 20),
                                                 nfit=nfit_for(DPS))
    mp.mp.dps = 60
    o2 = mp.mpf(ORC['laurent']['eps-2'])
    o1 = mp.mpf(ORC['laurent']['eps-1'])
    print(f"  [check] live threshold-Taylor K_2^(0) vs shipped 1-loop AMFlow: {fit_chk:.1f} d")
    print(f"  [check] live K_2^(1) (runtime 2-fold) vs stored AMFlow string"
          f" (held-out): {chk_K21:.1f} d")
    print(f"  eps^-2: -K_2^(0)            = {mp.nstr(I_em2, 36)}   (computed)")
    print(f"          oracle              = {mp.nstr(o2, 36)}")
    print(f"          agreement           = {agree_d(I_em2, o2):.2f} d   (archived: 45.76 d)")
    print(f"  eps^-1: -K_2^(1)+B[s_-1]^(0) = {mp.nstr(I_em1, 36)}   (computed;"
          f" K_2^(1) runtime-derived from the exact 2-fold)")
    print(f"          oracle              = {mp.nstr(o1, 36)}")
    print(f"          agreement           = {agree_d(I_em1, o1):.2f} d   (archived: 39.57 d)")
    print("  (the fixed-eps gate above certifies the full function; the layers are"
          " demos of the closed pole structure)")
    if args.fastcache_bank:
        d_em2, d_em1 = agree_d(I_em2, o2), agree_d(I_em1, o1)
        ok_bank = (mp.isfinite(mind) and float(mind) >= 30
                   and d_em2 >= 30 and d_em1 >= 30
                   and len(_fc_pts) == len(pts))
        if ok_bank:
            _, K2_1b = K2_threshold_eps01(max(90, DPS + 20))
            # banked_d derived from the STORED strings (what the fast path
            # re-checks) at the fast path's parse precision max(80, DPS+20)
            # -- in-memory agreements shift past the 0.05 d tight bar once
            # the string truncation enters (caught in a scratch run).
            s2, s1 = mp.nstr(I_em2, 52), mp.nstr(I_em1, 52)
            with mp.workdps(max(80, DPS + 20)):
                bd2 = float(agree_d(mp.mpf(s2), mp.mpf(ORC['laurent']['eps-2'])))
                bd1 = float(agree_d(mp.mpf(s1), mp.mpf(ORC['laurent']['eps-1'])))
            _fc_bank(_fc_pts,
                     {'I_em2': s2, 'I_em1': s1,
                      'banked_d_em2': bd2, 'banked_d_em1': bd1,
                      'fit_chk': float(fit_chk), 'chk_K21': float(chk_K21)},
                     mp.nstr(K2_1b, 80), DPS, time.time() - t_start,
                     f'python3 lbl3vp-evaluate.py --dps {DPS} --fastcache-bank')
        else:
            print("[fast-cache] bank REFUSED: >=30 d bars not met on every "
                  "gated eps + both pole layers")
    print(f"\ntotal wall time: {time.time()-t_start:.0f}s")


if __name__ == '__main__':
    main()