#!/usr/bin/env python3
"""C3 central-rung-mass double box: evaluate the explicit closed form at runtime
and reproduce the held-out gate vs AMFlow.

Self-contained (python3 + mpmath only). What it does, all live at runtime:
  1. Evaluates the explicit weight<=2 closed form g0,g1,g2 (valid in full 2D)
     at the three verification points (s=-1, t in {-1/2,-1,-2}).
  2. Evaluates the explicit s=-1 slice closed form for weights 3-4: the 15-word
     (w3) and 44-word (w4) Goncharov-polylog expressions of c3-dbox-expression.md,
     G(a_vec; x) with letters {0,-1,-1/2} in x=-t, computed here by a small
     GPL evaluator (Taylor-stepped iterated-integral ODE d/dx G(a,w;x) =
     G(w;x)/(x-a) transported from x=0, where every word with a nonzero trailing
     letter vanishes; poles at x=0,-1/2,-1 are off the Euclidean path x>0).
  3. Assembles the Laurent coefficients eps^-4..eps^0 of the top master
        J = e^{-2 eps gamma_E} eps^{-4} g(eps) / R,   R = s^2 (t-1),
     and compares each order to the independent held-out AMFlow oracle,
     with the agreement RECOMPUTED as -log10(|J - oracle|/|oracle|).
  4. NEW (2026-07-04 back-port of the closed 2-variable forms): evaluates the
     full-2D explicit closed forms for weights 3-4 at any Euclidean (s, t),
     s < 0, t < 0 — w3 as 39 exact-rational columns and w4 as 123 (t-sector)
     + 29 (pure-y block) exact-rational columns over products of GPL words
     G(a_vec; x = -t), x-letters {0, -1, s, s/(1-s)} (the tables' fifth
     x-letter -s^2 is used by no column and is not transported), and
     G(b_vec; y = -s), y-letters {0, 1, -1}, times log powers and the
     classical constants 1, zeta2, zeta3, zeta4 — and gates them live
     against 21 held-out AMFlow points (110-digit strings) on five slices
     s in {-4, -7/2, -3, -5/2, -3/2} never used in any fit. Data vendored in
     c3-dbox-2var.json (exact rationals; provenance inside).
  5. HARDENING (2026-07-08, "certified transport + raising Im-gate"):
     (a) every Taylor transport step (gpl_transport AND gpl_transport_2var)
     is CERTIFIED at runtime by the trailing-term geometric tail bound
     (r = 1/2 by the half-distance step rule; tail <= max trailing
     |c_n h^n| * r/(1-r) = max*1), growing N (x1.5) from the unchanged
     starting guess until the bound < 10^-(dps + TAIL_GUARD=12),
     RuntimeError at the cap 64*N0 — a step is never silently
     under-resolved; (b) the |Im| branch-cancellation residual is PROMOTED
     from printed diagnostic to a RAISING gate in _assemble_2var:
     |Im g_k|/max(|Re g_k|,1) > 10^-(dps - IM_MARGIN=8) for k = 2,3,4
     raises (the w4 t-sector is never assembled standalone here; on its own
     it carries a t-independent imaginary part that only cancels in the
     full g4, so it is not checked in isolation); (c) the
     default gate demo is FAIL-CLOSED: any gate miss exits nonzero.

The only literals kept below (and in c3-dbox-2var.json) are the AMFlow ORACLE
strings (provenance in the comment above the dict / the json _provenance
block) and the exact rational word coefficients of the closed forms
themselves (the analytic results being tested; derivation record:
c3-dbox-expression.md and the closure records of this work). The s=-1 slice
oracle strings carry ~45 significant digits, so ~45 d is the measurable
ceiling there; the 2-var oracle strings carry ~110 digits, so the 2-var gate
measures ~110 d live at --dps 130.

Interface (deps: python3 + mpmath, nothing else):
  python3 c3-dbox-evaluate.py                     # slice gate demo + dps-doubling
                                                  # + 2-var gate (5-pt subset, dps 130)
                                                  # FAIL-CLOSED: rc != 0 on any gate miss
  python3 c3-dbox-evaluate.py --gate-full         # 2-var gate on all 21 held-out pts
  python3 c3-dbox-evaluate.py --skip-2var         # historical slice-only demo
  python3 c3-dbox-evaluate.py --dps 120           # same, at doubled working precision
  python3 c3-dbox-evaluate.py --point -3/7        # J Laurent coeffs at t = -3/7, s = -1
  python3 c3-dbox-evaluate.py --point -2 --s -3 --dps 130   # off-slice: 2-var forms
or, from python, the callables are  evaluate(t, dps)  (s = -1 slice) and
evaluate2(s, t, dps)  (full 2D) below (pass arguments as strings 'p/q'/'-0.83'
or exact mpf; they set mp.dps BEFORE converting, avoiding the float/mpf
import-precision footgun).

Domain: s < 0, t < 0 (Euclidean). On the s = -1 slice the historical
slice-only path is used (x = -t > 0 keeps the GPL path clear of the letters
at x = -1/2, -1; any t < 0 works, cost grows only logarithmically in |t|).
Off the slice the 2-var closed forms are used; for s < -1 the y = -s
transport continues along a fixed upper-half-plane detour around the y = 1
letter, and the assembled weights are real (max |Im| printed AND gated: a
RAISING branch-cancellation gate fires if |Im g_k|/max(|Re g_k|,1) >
10^-(dps-8) for k=2,3,4). Precision is arbitrary: the Taylor length
N (starting guess; every step tail-certified, see gpl_transport), every
polylog and every constant recompute at the requested dps; no
cached numeric nodes anywhere.
"""
import itertools, json, math, os
from fractions import Fraction as Fr
from mpmath import (mp, mpf, mpc, log, pi, zeta, polylog, euler, factorial,
                    fabs, log10, im, re as mpre)

mp.dps = 60

# ---------------------------------------------------------------------------
# Weights 0-2: explicit 2D closed form (c3-dbox-expression.md, w<=2 section)
# ---------------------------------------------------------------------------
# Euclidean region s < 0: log(s) in the closed form means log(-s).
def g1(s, t):
    return -2*log(-s) - 4*log(1 - t)

def g2(s, t):
    # 2D weight-2 closed form. At s = -1 the term log(-s)*log(s+1) is the
    # limit 0 * log(0) = 0, handled explicitly.
    out = 2*log(-s)**2 + 8*log(-s)*log(1 - t) + 8*log(1 - t)**2
    if s != -1:
        out += -2*log(-s)*log(s + 1)
    out += -2*polylog(2, -s) + 10*polylog(2, t) + 2*polylog(2, s + t - s*t)
    out += -pi**2/6
    return out

# ---------------------------------------------------------------------------
# Weights 3-4 on the s = -1 slice: runtime GPL evaluator.
# Letters in x = -t:  t -> 0,  1-t -> -1,  2t-1 -> -1/2  (indices 0,1,2).
# G(a1..an; x) = int_0^x dt/(t-a1) G(a2..an; t),  G(a;x) = log(1 - x/a).
# ---------------------------------------------------------------------------
LETTERS = (mpf(0), mpf(-1), mpf(-1)/2)

def _all_words(maxw):
    words = []
    def rec(prefix, w):
        if len(prefix) == w:
            if prefix[-1] != 0:      # nonzero trailing letter => G(w;0) = 0
                words.append(tuple(prefix))
            return
        for i in range(3):
            rec(prefix + [i], w)
    for w in range(1, maxw + 1):
        rec([], w)
    return words

WORDS = _all_words(4)                # 80 words, suffix-closed
ORDER = sorted(WORDS, key=len)

# 2026-07-08 hardening: every
# transport step is CERTIFIED at runtime by the trailing-term geometric tail
# bound (half-distance step rule gives series ratio r <= 1/2 for every word,
# so  tail <= max trailing |c_n h^n| * r/(1-r) = max*1); N GROWS (x1.5) from
# the unchanged starting guess until the bound < 10^-(dps + TAIL_GUARD),
# RuntimeError at the cap -- a step is never silently under-resolved.
TAIL_GUARD = 12    # certified per-step Taylor tail < 10^-(dps + TAIL_GUARD)
TAIL_WINDOW = 8    # trailing terms scanned for the max (guards lone zeros)

def gpl_transport(targets):
    """G(word; x) for every word in WORDS at each target x > 0.

    Taylor-steps the iterated-integral ODE d/dx G(a,w;x) = G(w;x)/(x-a) from
    x = 0 (all words vanish there). Because T[1/(x-a)] at a regular center is
    a geometric series, each word's Taylor coefficients follow from a first-
    order recurrence, so a step costs O(N) per word. Step size = half the
    distance to the nearest pole seen from the center x0 (x = 0 for x0 > 0,
    x = -1/2 for the first step from x0 = 0, where the x = 0 letter is instead
    integrated term-wise), so the series ratio is <= 1/2 for EVERY letter and
    the step count to reach x is only logarithmic in x. N is the STARTING
    guess only: every step is certified by the trailing-term tail gate
    (see TAIL_GUARD above), growing N until the bound passes, RuntimeError
    at the cap 64*N0.
    """
    N0 = int((mp.dps + 12) / 0.301) + 2      # starting guess (formula kept)
    NCAP = 64 * N0
    tail_tol = mpf(10) ** (-(mp.dps + TAIL_GUARD))
    targets = sorted(mpf(t) for t in targets)
    vals = {w: mpf(0) for w in WORDS}
    x0 = mpf(0)
    out, ti = {}, 0
    while ti < len(targets):
        h = (x0 + mpf(1)/2) / 2 if x0 == 0 else x0 / 2
        hit = None
        if targets[ti] <= x0 + h:
            h = targets[ti] - x0
            hit = targets[ti]
        habs = fabs(h)
        N = N0
        while True:
            C = {(): [mpf(1)] + [mpf(0)] * N}
            for w in ORDER:
                a = LETTERS[w[0]]
                rest = C[w[1:]]
                c = [mpf(0)] * (N + 1)
                c[0] = vals[w]
                if a == 0 and x0 == 0:
                    # T[rest]/h is analytic at 0 (rest[0] = 0): integrate term-wise
                    for n in range(1, N + 1):
                        c[n] = rest[n] / n
                else:
                    d = x0 - a               # 1/(h+d) = sum_k p0 r^k h^k
                    p0, r = 1 / d, -1 / d
                    s = rest[0] * p0
                    c[1] = s
                    for n in range(2, N + 1):
                        s = r * s + p0 * rest[n - 1]
                        c[n] = s / n         # c[n] = s_{n-1}/n
                C[w] = c
            # certified trailing-term tail gate: half-distance step rule
            # gives series ratio r <= 1/2 for every word, so the truncated
            # tail obeys  tail <= max trailing |c_n h^n| * r/(1-r) = max*1.
            nlo = max(N - TAIL_WINDOW + 1, 1)
            hp0 = habs ** nlo
            bound = mpf(0)
            for w in ORDER:
                cc = C[w]
                hp = hp0
                for n in range(nlo, N + 1):
                    tt = fabs(cc[n]) * hp
                    if tt > bound:
                        bound = tt
                    hp *= habs
            if bound < tail_tol:
                break                     # step CERTIFIED at this N
            if N >= NCAP:
                raise RuntimeError(
                    f"c3-dbox transport tail gate FAILED: certified tail "
                    f"bound {mp.nstr(bound, 3)} >= tol "
                    f"{mp.nstr(tail_tol, 3)} (= 10^-(dps+{TAIL_GUARD}), "
                    f"dps={mp.dps}) at step x0={mp.nstr(x0, 8)}, "
                    f"h={mp.nstr(h, 8)} with N={N} (cap {NCAP}); "
                    f"raise the cap or investigate the pole geometry")
            N = N * 3 // 2 + 8            # grow until the bound passes
        for w in ORDER:
            c = C[w]
            v = mpf(0)
            for n in range(N, -1, -1):
                v = v * h + c[n]
            vals[w] = v
        x0 = x0 + h
        if hit is not None:
            out[hit] = dict(vals)
            ti += 1
    return out

# Exact rational word coefficients of the s=-1 slice closed form, transcribed
# from c3-dbox-expression.md (closure record of this work, 2026-06-26).
# Word letters: 0 = t, 1 = 1-t, 2 = 2t-1 (i.e. GPL letters 0, -1, -1/2).
W3_PURE = [(-8, (0, 0, 1)), (40, (0, 1, 1)), (-12, (0, 2, 1)),
           (40, (1, 0, 1)), (-64, (1, 1, 1)), (4, (1, 2, 1)),
           (-16, (2, 0, 1)), (8, (2, 1, 1)), (8, (2, 2, 1))]
W3_LOG2 = [(-12, (0, 2)), (4, (1, 2)), (8, (2, 2))]

W4_PURE = [(32, (0, 0, 1, 1)), (-20, (0, 0, 2, 1)), (164, (0, 1, 0, 1)),
           (-160, (0, 1, 1, 1)), (-28, (0, 1, 2, 1)), (-96, (0, 2, 0, 1)),
           (48, (0, 2, 1, 1)), (48, (0, 2, 2, 1)), (52, (1, 0, 0, 1)),
           (-160, (1, 0, 1, 1)), (24, (1, 0, 2, 1)), (-152, (1, 1, 0, 1)),
           (256, (1, 1, 1, 1)), (-8, (1, 1, 2, 1)), (32, (1, 2, 0, 1)),
           (-16, (1, 2, 1, 1)), (-16, (1, 2, 2, 1)), (-36, (2, 0, 0, 1)),
           (64, (2, 0, 1, 1)), (12, (2, 0, 2, 1)), (-28, (2, 1, 0, 1)),
           (-32, (2, 1, 1, 1)), (20, (2, 1, 2, 1)), (64, (2, 2, 0, 1)),
           (-32, (2, 2, 1, 1)), (-32, (2, 2, 2, 1))]
W4_LOG2 = [(-20, (0, 0, 2)), (-28, (0, 1, 2)), (48, (0, 2, 2)),
           (24, (1, 0, 2)), (-8, (1, 1, 2)), (-16, (1, 2, 2)),
           (12, (2, 0, 2)), (20, (2, 1, 2)), (-32, (2, 2, 2))]

def g3_slice(G):
    """g^(3) at s=-1 from the 15-word GPL closed form; G = word values at x=-t."""
    L2, Z3 = log(2), zeta(3)
    v = sum(mpf(c)*G[w] for c, w in W3_PURE)
    v += L2*sum(mpf(c)*G[w] for c, w in W3_LOG2)
    v += mpf(5)/3*pi**2*G[(1,)] + (-mpf(2)/3*pi**2 + 4*L2**2)*G[(2,)]
    v -= mpf(50)/3*Z3
    return v

def g4_slice(G):
    """g^(4) at s=-1 from the 44-word GPL closed form; G = word values at x=-t."""
    L2, Z3 = log(2), zeta(3)
    v = sum(mpf(c)*G[w] for c, w in W4_PURE)
    v += L2*sum(mpf(c)*G[w] for c, w in W4_LOG2)
    v += mpf(16)/3*pi**2*G[(0, 1)] + (-4*pi**2 + 24*L2**2)*G[(0, 2)]
    v += -mpf(10)/3*pi**2*G[(1, 1)] + (mpf(4)/3*pi**2 - 8*L2**2)*G[(1, 2)]
    v += -mpf(8)/3*pi**2*G[(2, 1)] + (mpf(8)/3*pi**2 - 16*L2**2)*G[(2, 2)]
    v += mpf(140)/3*Z3*G[(1,)] + (-8*Z3 - mpf(16)/3*L2**3)*G[(2,)]
    v += -pi**4/5 - 14*Z3*L2 + mpf(2)/3*pi**2*L2**2 - mpf(2)/3*L2**4 \
         - 16*polylog(4, mpf(1)/2)
    return v

def J_coeffs(s, t, G):
    """Laurent coefficients of J at eps^-4..eps^0 from g0..g4 (s = -1 slice)."""
    R = s**2*(t - 1)
    g = [mpf(1), g1(s, t), g2(s, t), g3_slice(G), g4_slice(G)]
    E = [(-2*euler)**k/factorial(k) for k in range(5)]   # e^{-2 eps gamma_E}
    return [sum(g[w]*E[n - w] for w in range(n + 1))/R for n in range(5)]

# ---------------------------------------------------------------------------
# ORACLE (legitimate held-out literal): independent AMFlow evaluation, dps 55,
# never used in the transport, fit, or slice-word closure. Provenance
# (repo artifacts, <archive>/Physics/Bootstrap/stage4/polylog-series/
#  C3-dbox-1internalmass-multiroot/):
#   t=-1/2: samples/cr_s-1_1_t-1_2.json, W2-detransport/graded_farm/cr_tm1_2_hi_out.json
#   t=-1  : W2-detransport/graded_farm/cr_tm1_hi_out.json  (gate: W2-detransport/CLOSE_t.json)
#   t=-2  : samples/cr_s-1_1_t-2_1.json, W2-detransport/graded_farm/cr_tm2_hi_out.json
# ---------------------------------------------------------------------------
S = mpf(-1)
AMFLOW = {  # t : {eps_order: amflow string}
    '-1/2': {-4: "-0.666666666666666666666666666666666666666666667",
             -3: "1.85086117482381549941671776134813427224818691",
             -2: "5.62569268601302292659205792861938881359426457",
             -1: "8.1712721018592524625748614890468985273128844",
              0: "4.45707528946417406900366064065150578639592329"},
    '-1':   {-4: "-0.5",
             -3: "1.9635100260214234794409763329987555671931596",
             -2: "4.66374006587296358961814961014709185907076684",
             -1: "4.66231282712134873339635966741485120473095613",
              0: "-2.4296731128324343743023059593891353272779447"},
    '-2':   {-4: "-0.333333333333333333333333333333333333333333333",
             -3: "1.84962682815850149559800170928496922689142697",
             -2: "3.13526226866885021938623257633193642923936119",
             -1: "0.86926666727117618178584722491150724604846346",
              0: "-7.6940541602659281372205091282438020380272746"}}
TVALS = {'-1/2': mpf(-1)/2, '-1': mpf(-1), '-2': mpf(-2)}

def digits(a, b):
    a, b = mpf(a), mpf(b)
    if a == b:
        return float(mp.dps)
    return float(-log10(fabs(a - b)/max(fabs(b), mpf('1e-30'))))

LABS = ['eps^-4', 'eps^-3', 'eps^-2', 'eps^-1', 'eps^0 ']

def _parse_mpf(sstr):
    """Parse 'p/q' or a decimal string at the CURRENT mp.dps."""
    sstr = str(sstr).strip()
    if '/' in sstr:
        num, den = sstr.split('/')
        return mpf(num.strip()) / mpf(den.strip())
    return mpf(sstr)

def _sig(sstr):
    """Significant digits carried by an oracle literal."""
    return len(sstr.replace('-', '').replace('.', '').lstrip('0'))

def evaluate(t, dps=60):
    """The callable f: Laurent coefficients [eps^-4..eps^0] of J at s = -1.

    Domain: t < 0 (Euclidean s = -1 slice; see module docstring). Sets
    mp.dps = dps FIRST, then converts t, so string/rational inputs keep full
    precision. Everything downstream (GPL transport length N, polylogs,
    constants) recomputes at this dps — arbitrary precision, no caches.
    """
    mp.dps = int(dps)
    t = _parse_mpf(t) if isinstance(t, str) else mpf(t)
    if not t < 0:
        raise ValueError("domain: t < 0 (Euclidean s = -1 slice)")
    G = gpl_transport([-t])[-t]
    return J_coeffs(mpf(-1), t, G)

# ===========================================================================
# FULL-2D SECTION (added 2026-07-04; closed forms of 2026-07-03): explicit
# 2-variable closed forms
# for w3 (39 exact-Q columns) and w4 (123-column t-sector + 29-term pure-y
# block, both exact; dual-decode confirmed). Constants: 1, zeta2, zeta3,
# zeta4 ONLY (classical). Data: c3-dbox-2var.json in this directory (exact
# rationals + the 21 held-out AMFlow gate strings; provenance inside).
# The s = -1 slice machinery above is untouched (historical page content);
# the two representations agree on the slice: the 123-term table restricted
# to s = -1 reproduces the 44-word slice form exactly (26 weight-four words,
# 9 log 2 words, 9 pi^2 / log^2 2 coefficients; exact rational arithmetic;
# c3-dbox-identities.py).
# ===========================================================================
_HERE = os.path.dirname(os.path.abspath(__file__))
with open(os.path.join(_HERE, 'c3-dbox-2var.json')) as _f:
    _D2 = json.load(_f)
# 6-tuple column schema: [coeff(Q), kt, xword, ks, yword, const];
# column value = coeff * log(x)^kt/kt! * G(xword;x) * log(y)^ks/ks! *
#                G(yword;y) * const,  x = -t, y = -s.
W3_2VAR = [(Fr(c), kt, tuple(wx), ks, tuple(wy), const)
           for c, kt, wx, ks, wy, const in _D2['w3_2var']]
W4T_2VAR = [(Fr(c), kt, tuple(wx), ks, tuple(wy), const)
            for c, kt, wx, ks, wy, const in _D2['w4t_2var']]
W4_PUREY_EXACT = [(Fr(c), kt, tuple(wx), ks, tuple(wy), const)
                  for c, kt, wx, ks, wy, const in _D2['w4_purey_exact']]
GATE_FEED = _D2['gate_feed']   # 21 held-out AMFlow points, NEVER used in any fit
# default 5-point gate subset (one per held-out slice of the first four + one; --gate-full = all 21)
GATE_SUBSET = {('-3', '-2'), ('-7/2', '-3'), ('-5/2', '-4'),
               ('-3/2', '-2/3'), ('-3', '-1/3')}

def gpl_transport_2var(letters, maxw, targets, waypoints=None):
    """G(word; z) for ALL words over `letters` (letters[0] must be 0) with
    nonzero trailing letter, length <= maxw, at each target, by the same
    Taylor-stepped ODE d G(a,w;z) = G(w;z)/(z-a) as gpl_transport above,
    generalized to arbitrary letter sets and complex waypoints (needed for
    the upper-half-plane detour when y = -s > 1). Step = half the distance
    to the nearest pole => series ratio <= 1/2; N = 3.33*(dps+15) terms is
    the STARTING guess only: every step is certified by the trailing-term
    tail gate (see TAIL_GUARD above), growing N (x1.5) until the bound
    < 10^-(dps + TAIL_GUARD), RuntimeError at the cap 64*N0.
    Targets are visited in the given order (callers pass ascending paths)."""
    letters = [x if isinstance(x, mpc) else mp.mpf(x) for x in letters]
    words = []
    for w in range(1, maxw + 1):
        for pre in itertools.product(range(len(letters)), repeat=w - 1):
            for last in range(1, len(letters)):
                words.append(pre + (last,))
    ORDER2 = sorted(words, key=len)
    N0 = int((mp.dps + 15) / 0.301) + 4      # starting guess (formula kept)
    NCAP = 64 * N0
    tail_tol = mpf(10) ** (-(mp.dps + TAIL_GUARD))
    vals = {w: mp.mpf(0) for w in words}
    z0 = mp.mpf(0)
    pts = list(waypoints or []) + list(targets)
    tstart = len(waypoints or [])
    out = {}
    poles = list(letters[1:])
    for pi_, ztgt in enumerate(pts):
        arrived = (z0 == ztgt)
        while not arrived:
            direction = ztgt - z0
            dist = abs(direction)
            if z0 == 0:
                h = min(abs(p) for p in poles) / 2
            else:
                ds = [abs(z0 - p) for p in poles] + [abs(z0)]
                h = min(ds) / 2
            if dist <= h:
                h_step = direction; arrived = True
            else:
                h_step = direction / dist * h
            habs = fabs(h_step)
            N = N0
            while True:
                C = {(): [mp.mpf(1)] + [mp.mpf(0)] * N}
                for w in ORDER2:
                    a = letters[w[0]]
                    rest = C[w[1:]]
                    c = [mp.mpf(0)] * (N + 1)
                    c[0] = vals[w]
                    if a == 0 and z0 == 0:
                        for n in range(1, N + 1):
                            c[n] = rest[n] / n
                    else:
                        d = z0 - a
                        p0, r = 1 / d, -1 / d
                        ssum = rest[0] * p0
                        c[1] = ssum
                        for n in range(2, N + 1):
                            ssum = r * ssum + p0 * rest[n - 1]
                            c[n] = ssum / n
                    C[w] = c
                # certified trailing-term tail gate: half-distance step rule
                # gives series ratio r <= 1/2 for every word, so the truncated
                # tail obeys  tail <= max trailing |c_n h^n| * r/(1-r) = max*1.
                nlo = max(N - TAIL_WINDOW + 1, 1)
                hp0 = habs ** nlo
                bound = mp.mpf(0)
                for w in ORDER2:
                    cc = C[w]
                    hp = hp0
                    for n in range(nlo, N + 1):
                        tt = fabs(cc[n]) * hp
                        if tt > bound:
                            bound = tt
                        hp *= habs
                if bound < tail_tol:
                    break                     # step CERTIFIED at this N
                if N >= NCAP:
                    raise RuntimeError(
                        f"c3-dbox transport tail gate FAILED: certified tail "
                        f"bound {mp.nstr(bound, 3)} >= tol "
                        f"{mp.nstr(tail_tol, 3)} (= 10^-(dps+{TAIL_GUARD}), "
                        f"dps={mp.dps}) at step z0={mp.nstr(z0, 8)}, "
                        f"h={mp.nstr(h_step, 8)} with N={N} (cap {NCAP}); "
                        f"raise the cap or investigate the pole geometry")
                N = N * 3 // 2 + 8            # grow until the bound passes
            for w in ORDER2:
                cc = C[w]
                v = mp.mpf(0)
                for n in range(N, -1, -1):
                    v = v * h_step + cc[n]
                vals[w] = v
            z0 = ztgt if arrived else z0 + h_step
        if pi_ >= tstart:
            out[pi_ - tstart] = dict(vals)
    return out

def _xylet(s):
    """x-letters (dlog images of alphabet letters t, 1-t, s+t, s+t-st in
    x = -t) and y-letters (letters s, s+1, 1-s in y = -s). The tables encode a
    fifth x-letter -s^2 (index 4) that no column uses; it is not transported."""
    xlet = [mp.mpf(0), mp.mpf(-1), s, s / (1 - s)]
    ylet = [mp.mpf(0), mp.mpf(1), mp.mpf(-1)]
    return xlet, ylet

def _gy_transport(ylet, y0, maxw=4):   # maxw 4: pure-y w4 block has len-4 words
    if y0 < 1:
        return gpl_transport_2var(ylet, maxw, [y0])[0]
    # s < -1: fixed upper-half-plane detour around the y = 1 letter
    return gpl_transport_2var(ylet, maxw, [y0],
                              waypoints=[mpc(mpf('0.5'), mpf('0.6')),
                                         mpc(y0, mpf('0.6'))])[0]

def _term_sum(table, gx, gy, logx, logy, cv):
    tot = mp.mpf(0)
    for c, kt, wx, ks, wy, const in table:
        v = cv[const] * mp.mpf(c.numerator) / c.denominator
        if kt: v *= logx**kt / math.factorial(kt)
        if wx: v *= gx[wx]
        if ks: v *= logy**ks / math.factorial(ks)
        if wy: v *= gy[wy]
        tot += v
    return tot

IM_MARGIN = 8      # raise if |Im g_k|/max(|Re g_k|,1) > 10^-(dps-IM_MARGIN);
                   # set from a measured dps-60 run: worst residual 3.9e-60
                   # = 10^-(dps-0.6) (g3 at (-3,-2)); see gate in _assemble_2var

def _assemble_2var(gx, gy, s, t):
    """(g1, g2, g3, g4, imax) from transported word dicts at (s, t).
    g4 = t-sector (123 exact columns) + pure-y block C4(y) (29 exact terms);
    all four weights are explicit closed forms over Q[zeta2, zeta3, zeta4].
    imax = max |Im| over g2..g4: the branch-cancellation certificate for the
    upper-half-plane detour (s < -1); the assembled weights are real —
    printed AND gated (RAISING, see IM_MARGIN above)."""
    x0, y0 = -t, -s
    logx, logy = log(x0), log(y0)
    cv = {'1': mp.mpf(1), 'z2': pi**2 / 6, 'z3': zeta(3), 'z4': pi**4 / 90}
    gg1 = -2 * log(-s) - 4 * log(1 - t)
    # w2: the 8-term 2D closed form (c3-dbox-expression.md) rewritten in
    # G-column language so
    # every piece is real-by-transport for s < -1 (log(s+1) = G(1;y) cont.;
    # Li2(-s) = -G(0,1;y); Li2(s) = -G(0,-1;y); Li2(t) = -G(0,-1;x);
    # Li2(s+t-st) = -log(1-s) G(p3;x) - G(p3,-1;x) + Li2(s)):
    gg2 = (2 * logy**2 - 2 * logy * gy[(1,)]
           + 8 * logy * gx[(1,)] + 16 * gx[(1, 1)]
           + 2 * gy[(0, 1)] - 2 * gy[(0, 2)]
           - 10 * gx[(0, 1)]
           - 2 * gy[(2,)] * gx[(3,)] - 2 * gx[(3, 1)]
           - pi**2 / 6)
    gg3 = _term_sum(W3_2VAR, gx, gy, logx, logy, cv)
    gg4 = (_term_sum(W4T_2VAR, gx, gy, logx, logy, cv)
           + _term_sum(W4_PUREY_EXACT, gx, gy, logx, logy, cv))
    imax = max(fabs(im(gg2)), fabs(im(gg3)), fabs(im(gg4)))
    # RAISING branch-cancellation gate (2026-07-08): the assembled
    # g2/g3/g4 must be real (s<-1 detour
    # Im parts cancel by construction); the residual is pure roundoff,
    # measured at dps 60 over all 20 gate points: |Im|/|Re| in
    # [1.9e-62, 3.9e-60] ~ 10^-(dps-1) worst.  Threshold 10^-(dps-IM_MARGIN)
    # leaves ~10^7 headroom yet catches any genuine cancellation failure
    # (which contaminates at O(1)).  |Re| floored at 1 so an accidental real
    # zero cannot false-fire the ratio; below the floor the gate binds the
    # ABSOLUTE residual, which is the roundoff scale anyway.  NOTE: the w4
    # t-sector alone is NOT gated (it legitimately
    # carries a t-independent Im part that only cancels in fixed-s
    # t-differences and in full g4 = t-sector + pure-y); this file never
    # assembles it standalone, so gating g2/g3/g4 preserves that scoping.
    im_tol = mpf(10) ** (-(mp.dps - IM_MARGIN))
    for _name, _gv in (('g2', gg2), ('g3', gg3), ('g4', gg4)):
        _rel = fabs(im(_gv)) / max(fabs(mpre(_gv)), mp.mpf(1))
        if _rel > im_tol:
            raise RuntimeError(
                f"c3-dbox branch-cancellation gate FAILED at (s,t)=("
                f"{mp.nstr(s, 8)},{mp.nstr(t, 8)}), dps={mp.dps}: "
                f"|Im {_name}| = {mp.nstr(fabs(im(_gv)), 3)}, "
                f"|Im|/max(|Re|,1) = {mp.nstr(_rel, 3)} > "
                f"10^-(dps-{IM_MARGIN}) = {mp.nstr(im_tol, 3)}; the "
                f"transported y-words did not cancel to a real value")
    return gg1, gg2, gg3, gg4, imax

def evaluate2(s, t, dps=60):
    """The 2D callable f: Laurent coefficients [eps^-4..eps^0] of J at any
    Euclidean (s, t), s < 0, t < 0, plus the |Im| branch certificate.
    At s = -1 exactly it delegates to the historical slice evaluator (the
    2-var y-words degenerate at y = 1). Sets mp.dps FIRST, then converts."""
    mp.dps = int(dps)
    s = _parse_mpf(s) if isinstance(s, str) else mpf(s)
    t = _parse_mpf(t) if isinstance(t, str) else mpf(t)
    if not (s < 0 and t < 0):
        raise ValueError("domain: s < 0, t < 0 (Euclidean)")
    if s == -1:
        return evaluate(t, dps), mpf(0)
    xlet, ylet = _xylet(s)
    gx = gpl_transport_2var(xlet, 4, [-t])[0]
    gy = _gy_transport(ylet, -s)
    gg1, gg2, gg3, gg4, imax = _assemble_2var(gx, gy, s, t)
    g = [mpf(1), mpre(gg1), mpre(gg2), mpre(gg3), mpre(gg4)]
    R = s**2 * (t - 1)
    E = [(-2 * euler)**k / factorial(k) for k in range(5)]
    return [sum(g[w] * E[n - w] for w in range(n + 1)) / R
            for n in range(5)], imax

def run_gate_2var(full=False):
    """Live 2-var gate: recompute w3 and full w4 at the held-out AMFlow points
    (110-digit strings; slices s in {-4,-7/2,-3,-5/2,-3/2} NEVER used in any fit
    — the oracle points that fit the pure-y block are NOT in this set) and
    return (s, t, w3 digits, w4 digits, |Im| certificate) rows. One transport
    pair per s-slice; x-transport is multi-target (ascending in x = -t)."""
    pts = [p for p in GATE_FEED
           if full or (p['s'], p['t']) in GATE_SUBSET]
    groups = {}
    for p in pts:
        groups.setdefault(p['s'], []).append(p)
    rows = []
    for sstr, grp in groups.items():
        s = _parse_mpf(sstr)
        grp = sorted(grp, key=lambda p: -_parse_mpf(p['t']))  # ascend in x=-t
        xlet, ylet = _xylet(s)
        gxs = gpl_transport_2var(xlet, 4, [-_parse_mpf(p['t']) for p in grp])
        gy = _gy_transport(ylet, -s)
        for i, p in enumerate(grp):
            t = _parse_mpf(p['t'])
            _g1, _g2, gg3, gg4, imax = _assemble_2var(gxs[i], gy, s, t)
            rows.append((p['s'], p['t'],
                         digits(mpre(gg3), p['g3']),
                         digits(mpre(gg4), p['g4']),
                         mp.nstr(max(imax, fabs(im(gg4))), 3)))
    rows.sort(key=lambda r: (_parse_mpf(r[0]), _parse_mpf(r[1])))
    return rows

if __name__ == '__main__':
    import argparse, time
    ap = argparse.ArgumentParser(
        description="C3 double box closed form at s = -1: gate demo, or "
                    "evaluate at any t < 0 and any dps (see module docstring).")
    ap.add_argument('--point', metavar='T', default=None,
                    help="evaluate at t = T ('p/q' or decimal), t < 0")
    ap.add_argument('--s', metavar='S', default='-1',
                    help="s value for --point (default -1 = historical slice; "
                         "any s < 0 uses the 2-var closed forms)")
    ap.add_argument('--dps', type=int, default=60,
                    help="working precision, decimal digits (default 60)")
    ap.add_argument('--skip-doubling', action='store_true',
                    help="skip the dps-doubling demo in gate mode")
    ap.add_argument('--gate-full', action='store_true',
                    help="2-var gate on all 21 held-out points (default: "
                         "5-point subset, one per held-out slice of the first four + one)")
    ap.add_argument('--skip-2var', action='store_true',
                    help="skip the 2-var gate section (historical slice-only demo)")
    # argparse mistakes '-3/7' for an option flag (the slash defeats its
    # negative-number heuristic); join '--point -3/7' into '--point=-3/7'
    # (same for '--s -3').
    import sys
    argv = sys.argv[1:]
    for flag in ('--point', '--s'):
        for i, a in enumerate(argv):
            if a == flag and i + 1 < len(argv):
                argv[i:i + 2] = [flag + '=' + argv[i + 1]]
                break
    args = ap.parse_args(argv)
    wall0 = time.perf_counter()
    mp.dps = args.dps

    if args.point is not None and _parse_mpf(args.s) != -1:
        # off-slice point: the 2-var closed forms (one transport pair,
        # reused for both the Laurent coefficients and the gate comparison)
        s, t = _parse_mpf(args.s), _parse_mpf(args.point)
        if not (s < 0 and t < 0):
            raise SystemExit("domain: s < 0, t < 0 (Euclidean)")
        xlet, ylet = _xylet(s)
        gx = gpl_transport_2var(xlet, 4, [-t])[0]
        gy = _gy_transport(ylet, -s)
        gg1, gg2, gg3, gg4, imax = _assemble_2var(gx, gy, s, t)
        g = [mpf(1), mpre(gg1), mpre(gg2), mpre(gg3), mpre(gg4)]
        R = s**2 * (t - 1)
        E = [(-2*euler)**k/factorial(k) for k in range(5)]
        c = [sum(g[w]*E[n - w] for w in range(n + 1))/R for n in range(5)]
        print(f"C3 double box J at s = {args.s}, t = {args.point}  "
              f"(dps {args.dps}; 2-var closed forms, w3+w4 ALL EXACT)")
        for n in range(5):
            print(f"  {LABS[n]} = {mp.nstr(c[n], args.dps)}")
        print(f"  |Im| branch-cancellation residual: {mp.nstr(imax, 3)}")
        ref = next((p for p in GATE_FEED
                    if _parse_mpf(p['s']) == s and _parse_mpf(p['t']) == t), None)
        if ref is not None:
            cap = min(_sig(ref['g3']), _sig(ref['g4']))
            print(f"\n  held-out gate point (oracle strings carry ~{cap} sig "
                  f"digits; never used in any fit):")
            print(f"    w3  agree {digits(mpre(gg3), ref['g3']):.2f} d")
            print(f"    w4  agree {digits(mpre(gg4), ref['g4']):.2f} d")
        else:
            print("\n  (no stored oracle at this point — this value is the "
                  "closed form's prediction)")
        print(f"\nwall time: {time.perf_counter() - wall0:.2f} s")
        raise SystemExit(0)

    if args.point is not None:
        c = evaluate(args.point, args.dps)
        print(f"C3 double box J at s = -1, t = {args.point}  (dps {args.dps})")
        for n in range(5):
            print(f"  {LABS[n]} = {mp.nstr(c[n], args.dps)}")
        key = next((k for k, tv in TVALS.items() if tv == _parse_mpf(args.point)), None)
        if key is not None:
            cap = max(_sig(v) for v in AMFLOW[key].values())
            print(f"\n  gate point: live agreement vs held-out AMFlow oracle "
                  f"(strings carry {cap} sig digits — that is the cap):")
            for n in range(5):
                print(f"    {LABS[n]}  agree {digits(c[n], AMFLOW[key][n - 4]):.2f} d")
        else:
            print("\n  (no stored oracle at this point — oracle points are "
                  "t = -1/2, -1, -2; this value is the closed form's prediction)")
        print(f"\nwall time: {time.perf_counter() - wall0:.2f} s")
        raise SystemExit(0)

    # fail-closed gate demo (2026-07-08 hardening): every printed agreement is
    # also CHECKED; any miss is collected and the demo exits nonzero.
    # Thresholds are slack-below-measured but strictly binding (a 1e-30
    # mutation of any stored reference lands ~29-30 d, far below either bar):
    #   slice: measured live worst 44.62 d at dps 60 (45-46 d oracle cap)
    #   2-var: measured live worst 109.37 d at dps 130 (110-digit oracle cap)
    fails = []
    SLICE_GATE_MIN = min(40.0, args.dps - 20)
    GATE2VAR_MIN = 100.0

    print("C3 double box, s = -1: closed form (w<=2 2D + w3/w4 slice GPL words)")
    print(f"evaluated at runtime (dps {mp.dps}) vs independent held-out AMFlow oracle\n")
    GX = gpl_transport([-t for t in TVALS.values()])   # x = -t > 0
    for key, t in TVALS.items():
        c = J_coeffs(S, t, GX[-t])
        print(f"  t = {key}:")
        for n in range(5):
            ref = AMFLOW[key][n - 4]
            d = digits(c[n], ref)
            print(f"    {LABS[n]}  closed-form {mp.nstr(c[n], 22)}...  "
                  f"AMFlow {ref[:22]}...  agree {d:.2f} d")
            if not d > SLICE_GATE_MIN:
                fails.append(f"slice gate t={key} {LABS[n].strip()}: "
                             f"{d:.2f} d <= {SLICE_GATE_MIN:.1f} d")
        print()
    print("All five orders are computed here from the analytic closed form;")
    print("agreement is recomputed, capped ~45-46 d by the stored oracle strings")
    print("(the archived dps-55 gate: w2 40.4-54.9 d, w3 39.1-54.8 d, w4 39.4-54.9 d;")
    print(" slice words vs 4 fresh AMFlow points: w3 >= 108.7 d, w4 >= 108.8 d).")
    print(f"\ngate demo wall time: {time.perf_counter() - wall0:.2f} s")

    if not args.skip_doubling:
        key = '-1'
        cap = max(_sig(v) for v in AMFLOW[key].values())
        print(f"\ndps-doubling demo at gate point t = {key} "
              f"(oracle strings carry {cap} sig digits):")
        cs = {}
        for d in (max(args.dps // 2, 20), args.dps, 2 * args.dps):
            cs[d] = evaluate(TVALS[key], d)
            agr = min(digits(cs[d][n], AMFLOW[key][n - 4]) for n in range(5))
            note = (f"  <- capped: stored oracle carries only {cap} digits"
                    if agr > cap - 3 else "")
            print(f"  dps {d:4d}: min live agreement vs oracle = {agr:.2f} d{note}")
        mp.dps = 2 * args.dps
        cross = min(float(-log10(fabs(cs[args.dps][n] - cs[2 * args.dps][n])
                                 / fabs(cs[2 * args.dps][n]))) for n in range(1, 5))
        print(f"  self-check J(dps={args.dps}) vs J(dps={2 * args.dps}): "
              f"min agreement {cross:.2f} d")
        if not cross > args.dps - 15:
            fails.append(f"dps-doubling self-check: {cross:.2f} d <= "
                         f"{args.dps - 15} d (= dps - 15)")
        print(f"  (agreement vs oracle grows with dps until the {cap}-digit oracle-"
              "string cap;\n   the evaluation itself keeps sharpening — the "
              "self-check digits track dps)")

    if not args.skip_2var:
        gdps = max(args.dps, 130)   # 2-var oracle strings carry ~110 digits
        mp.dps = gdps
        npts = len(GATE_FEED) if args.gate_full else len(GATE_SUBSET)
        print(f"\nfull-2D gate: explicit 2-var closed forms (w3: 39 exact "
              f"columns; w4: 123-column\nt-sector + 29-term pure-y block, all "
              f"exact, constants zeta2/zeta3/zeta4 only)\nrecomputed live at "
              f"dps {gdps} vs {npts} held-out AMFlow points (110-digit "
              f"strings,\nslices never used in any fit"
              + ("" if args.gate_full else f"; --gate-full for all {len(GATE_FEED)}") + "):")
        t2 = time.perf_counter()
        rows2 = run_gate_2var(full=args.gate_full)
        for sstr, tstr, d3, d4, imstr in rows2:
            print(f"  s = {sstr:4s}  t = {tstr:4s}   w3 agree {d3:7.2f} d   "
                  f"w4 agree {d4:7.2f} d   |Im| {imstr}")
            for wname, dd in (('w3', d3), ('w4', d4)):
                if not dd > GATE2VAR_MIN:
                    fails.append(f"2-var gate s={sstr} t={tstr} {wname}: "
                                 f"{dd:.2f} d <= {GATE2VAR_MIN:.1f} d")
        d3s = [r[2] for r in rows2]; d4s = [r[3] for r in rows2]
        print(f"  {len(rows2)} pts: w3 {min(d3s):.2f}-{max(d3s):.2f} d, "
              f"w4 {min(d4s):.2f}-{max(d4s):.2f} d "
              f"(oracle-string-capped; agreement recomputed live)")
        print(f"2-var gate wall time: {time.perf_counter() - t2:.2f} s")

    print(f"\ntotal wall time: {time.perf_counter() - wall0:.2f} s")
    if fails:
        print(f"\nGATE DEMO: FAIL (fail-closed) -- {len(fails)} gate miss(es):")
        for f in fails:
            print("  " + f)
        raise SystemExit(1)
    print("\nGATE DEMO: PASS (fail-closed: any gate miss exits nonzero)")
