#!/usr/bin/env python3
"""AMFlow-free definition + arbitrary-precision evaluator for the equal-mass
sunrise eps^3/eps^4 boundary constants B3, B4 (and the whole tower B0..B4).

CONTEXT.  The blind bootstrap of the equal-mass two-loop sunrise
J111(2-2eps, t), m^2 = mu^2 = 1, t = p^2 (engine norm, no e^{gamma*eps};
J111 = -S111, AW sign) closed the function as a Gamma_1(6) Eichler-word tower

    e^{2*gamma*eps} J111(2-2eps, t) = (2/sqrt3) u(q) e^{-eps*wt(q)}
                                      * Sum_{j>=0} eps^j H_j(q),

with q the Gamma_1(6) nome (hauptmodul t(q) = 9 eta(6t)^8 eta(t)^4 /
(eta(2t)^8 eta(3t)^4)), u = (sqrt3/2)*psi1-hat the normalized holomorphic
period (u(0) = 1), wt = D^{-1}((k2+1)/2) the exact eps-factorizing rotation
(wt(0) = 0: (k2+1)/2 vanishes at q = 0), and H_j exact Q(sqrt3)-combinations
of iterated q-integrals I(l1,...,ln; q) over the derived alphabet
{one, k2h, k3, k4E, c4}, tangential base point q = 0.  Recorded structure:
  <archive>/Physics/Bootstrap/stage1/results/sunrise/exact/H_words.json
  <archive>/Physics/Bootstrap/stage1/results/sunrise/RESULT.md
In H_words.json the EMPTY word of H_j carries exactly the coefficient B_j.
B3, B4 were until now pinned only as 500-digit protocol numerics fitted
against the AMFlow sample runs (patched amflow-cpp; provenance:
  <archive>/Physics/Bootstrap/stage1/samples/sunrise/README.md).

AMFLOW-FREE DEFINITION (this file).  All regularized iterated q-integrals
vanish at the tangential base point q = 0, and u(0) = 1, wt(0) = 0,
t(q=0) = 0.  Hence the tower collapses at the MUM cusp q -> 0 to

    H_j(0) = B_j,   and   Sum_j eps^j B_j = (sqrt3/2) e^{2*gamma*eps}
                                            * J111(2-2eps, t=0).

J111 at t = p^2 = 0 is the EQUAL-MASS TWO-LOOP VACUUM SUNSET (all propagator
masses 1, below every threshold, analytic at t=0; the t^eps branch of the
{0, eps} indicials is excluded by analyticity — same argument the original
recursion used to fix the log freedom).  With S_vac = -J111(t=0) > 0
(Euclidean AW sign):

  DEFINITION:   B_j := -(sqrt3/2) * [eps^j] ( e^{2*gamma*eps} S_vac(eps) )

and S_vac(eps) has two explicit AMFlow-free forms, both implemented here:

 (I) convergent Bessel-moment integral (used as runtime cross-check):
       S_vac(eps) = 2^{2-eps}/Gamma(1-eps) * int_0^inf x^{1+eps} K_eps(x)^3 dx
     (position space: propagator = (2pi)^{-d/2} x^{1-d/2} K_{d/2-1}(x),
      d = 2-2eps; angular volume + engine normalization give the prefactor).

 (II) convergent Gamma-term series (used for the eps-jets; derivation:
      one-loop bubble B(p^2) = Gamma(2-d/2) * int_0^1 dx [1+x(1-x)p^2]^{d/2-2},
      Mellin-Barnes for the bracket, int_0^1 (x(1-x))^z dx = Gamma(z+1)^2/
      Gamma(2z+2), then the p-integral against 1/(p^2+1) in closed Gamma form;
      closing the 1-fold MB rightwards picks up z = n and z = n+eps, and the
      Gamma-reflection collapses each pair to:

      S_vac(eps) = pi/(sin(pi*eps) * Gamma(1-eps)) * Sum_{n>=0}
          [ n! Gamma(1+n+eps)/Gamma(2n+2)
            - Gamma(1+n+2eps) Gamma(1+n+eps)/Gamma(2n+2+2eps) ].

      The bracket vanishes at eps = 0 term by term, cancelling the
      pi/sin(pi*eps) pole EXACTLY in jet arithmetic (implemented as
      c_n*(exp(u_n)-exp(v_n)) with u_n, v_n log-Gamma jets that share a zero
      constant term — no numerical cancellation).  Terms decay like 4^{-n},
      so N ~= 1.67*dps terms give dps digits; all polygamma data at integer
      arguments are generated by exact recurrences from zeta values.

NO AMFlow-derived number enters the definition or the computation path.
The AMFlow-era 500-digit protocol numerics remain ONLY as held-out gate
literals below (byte-copied; grep them at
  <archive>/Physics/Bootstrap/stage1/results/sunrise/fits/fit_results.json
  keys results.B1/B2/B3/B4, prefix "NUMERIC: [").

POSITIVE CONTROLS (independent closed forms, run every time):
  B0 = -(3/2) sqrt3 L(chi-3,2)                       [exact, this work]
  B1 = -6 ImLi3(1-r3) - (1/2) pi log^2 3 - (5/54) pi^3
  B2 = 12 ImLi4(1-r3) + (3/2) sqrt3 L(chi-3,4) - (1/4) sqrt3 pi^2 L(chi-3,2)
       + (5/54) pi^3 log3 + (1/6) pi log^3 3
with r3 = e^{2*pi*i/3} (so 1-r3 = sqrt3 * e^{-i*pi/6}); L-values via Hurwitz
zeta.  These are the B1/B2 closed forms of this work (RESULT.md) and match the
series here to full working precision.

FOOTGUNS BAKED IN: mp.dps is set inside main() after argument parsing (never
at import); series length auto-sized from dps;
--check reruns everything at dps+60 and diffs (two-precision rule).

CLI:
    python3 sunrise-B34.py --dps 520            # evaluate + controls + gate
    python3 sunrise-B34.py --dps 520 --check    # + two-precision rerun/diff
    python3 sunrise-B34.py --dps 520 --quad     # + Bessel-integral cross-check
    python3 sunrise-B34.py --dps 300 --pslq     # + PSLQ naming attempt (honest
                                                #   null expected; controls first)
Dependency: pip mpmath only.
"""
import argparse
import time

import mpmath as mp

NJET = 8  # eps^0..eps^7 carried internally (B0..B4 needed; 3 guard orders)

# ---- axis-3 hardening knobs (2026-07-06; seeds/margins, fail-closed) --------
TAIL_GUARD = 8    # certified relative series tail must clear 10^-(dps+TAIL_GUARD)
CERT_CAP = 8      # max x1.5 series extensions in the certified loop (fail-closed)
GATE_MARGIN = 8   # held-out gate table RAISES below min(dps, digits) - GATE_MARGIN
CRANK_STEP = 40   # in-run D -> D+CRANK_STEP two-precision crank gate (raising)
CRANK_MARGIN = 4  # crank gate RAISES below dps - CRANK_MARGIN

# ---------------------------------------------------------------------------
# Held-out gate literals (AMFlow protocol numerics; NOT used in any
# computation below — comparison only).  Byte-copied from
# stage1/results/sunrise/fits/fit_results.json : results.B1/B2/B3/B4
# ("NUMERIC: [...]", staggered-goal PSLQ-anchor fit at t=-1, 1/2; engine =
# patched amflow-cpp per stage1/samples/sunrise/README.md).
# ---------------------------------------------------------------------------
GATE = {
    1: "3.4966434998229078606310448758241536228847941038927438085663018320609100789928608147103125043901883847137163043440340933851817466539440620607576810208877424457611180784573464183332544071667382883599725992725537894091511452625634721254096562560043432194844467373857445637334321796132724981312570086401452571995666599435954946246284669381312569727915081234350968902591458901206389033399374402230823399861483685530360152531837364464667033988903922780218616483105980675545321822703402953391160260867039453863450626366146161",
    2: "-10.30408706529854878384416029021417491888216163689755231223899338688874806721883707716091215080118185087812619295057740941135342677940996716179276091094300374566382745191400687469487029897962930799224226295919299039188718830029973940947082755083298206215999304782914080925288967198843156506128611367460713873347760650281675063031379011071548190660730371812542113839282304895934453636205201381609315095358906374875284042140327860665677895475446373568650490144661781832910220552353342713817135813583003711932782278633444",
    3: "21.257175978757173084657206013858551747273653522272268718848407280046917795440032078265667581622465661731460954905494235154629580522729952777561837878760711899916985849411044288955067112007143842338841663398734741749185162912570854434690112744734778574346162882763713138731446509959777234055742709783310058707426668986818285545965361437647189112082892677703368066594464093842021580774754804782071721211158674279535661472724106924465273592848916630123248798941101204554701671946892976132185833438442393972942119889809243",
    4: "-45.83139704914937206518906003650497171490939802849120602897162307010260150147736501741430193371762858554227617308042444237143226510796367765714593585760757040473742595015698620029534655115969936297606392988220581339208286525239004236931521007077350741313920603680633487157004688018826180752816353057134994669937543810203878773848426827718920474972590401059250755953811780853708440460644704165219321997020636106467178518933285988396174760532432584924873492289649393554064579697913352878617882108177368684943896851295417",
}


# ------------------------------ jet arithmetic ------------------------------
def jmul(a, b):
    r = [mp.mpf(0)] * NJET
    for i, ai in enumerate(a):
        if ai:
            for k in range(NJET - i):
                if b[k]:
                    r[i + k] += ai * b[k]
    return r


def jexp(u):
    """exp of a jet with u[0] == 0 (returns jet with [0] == 1, exact)."""
    assert u[0] == 0
    r = [mp.mpf(1)] + [mp.mpf(0)] * (NJET - 1)
    t = list(r)
    for k in range(1, NJET):
        t = [x / k for x in jmul(t, u)]
        r = [x + y for x, y in zip(r, t)]
    return r


def jinv(a):
    r = [mp.mpf(0)] * NJET
    r[0] = 1 / a[0]
    for k in range(1, NJET):
        r[k] = -sum(a[i] * r[k - i] for i in range(1, k + 1)) / a[0]
    return r


# --------------------------- the defining series ----------------------------
def svac_jets(nterms=None, tail_out=None):
    """Jet (eps^0..eps^{NJET-1}) of e^{2*gamma*eps} * S_vac(eps), form (II).

    nterms auto-sized from mp.dps (terms decay 4^{-n}: log10(4) digits/term).
    All polygamma values at integer arguments by exact recurrence:
        psi^{(k)}(x+1) = psi^{(k)}(x) + (-1)^k k! / x^{k+1},
    seeded at psi^{(k)}(1) = -gamma, (-1)^{k+1} k! zeta(k+1).

    tail_out (axis-3, 2026-07-06): optional dict, filled with a CERTIFIED
    truncation-tail bound (a BOUND, not an estimate):
      * c_n <= c_N * 4^-(n-N) for n >= N, exactly (term ratio
        (n+1)^2/((2n+2)(2n+3)) < 1/4; post-loop c holds c_N).
      * |jexp(a)_k| <= exp(sum_i |a_i|) (multinomial-exp domination).
      * post-loop state ps1/ps2 hold psi^{(k)}(N+1) / psi^{(k)}(2N+2);
        |psi^{(k>=1)}| decreases in the argument, psi increases with
        psi(x + j) <= psi(x) + j/x, so for n = N + j every jet entry sum is
        <= its n=N value + 5j/(N+1)  =>  per-jet-order
        |term_n| <= c_N 4^-j (e^Au + e^Av) e^{5j/(N+1)},  j = n - N >= 0,
        summing to the geometric bound below (g = e^{5/(N+1)}/4 < 1/3).
    The bound propagates to the returned jets via the l1-norms of the
    (1/Gamma, sin, e^{2 gamma eps}) jet factors (convolution triangle
    inequality); caller gates it fail-closed."""
    N = nterms if nterms else int(mp.ceil(mp.mpf(mp.mp.dps + 12) / mp.log10(4))) + 30
    kfact = [mp.gamma(k + 1) for k in range(NJET)]
    # psi^{(k)}(n+1) and psi^{(k)}(2n+2), running in n (k = 0..NJET-2)
    ps1 = [-mp.euler] + [(-1) ** (k + 1) * kfact[k] * mp.zeta(k + 1)
                         for k in range(1, NJET - 1)]
    ps2 = [ps1[k] + (-1) ** k * kfact[k] for k in range(NJET - 1)]  # psi^{(k)}(2)
    c = mp.mpf(1)                       # c_n = n!^2/(2n+1)!;  c_0 = 1
    tot = [mp.mpf(0)] * NJET
    for n in range(N):
        u = [mp.mpf(0)] * NJET          # lnGamma(1+n+eps) - lnGamma(1+n)
        v = [mp.mpf(0)] * NJET          # ln[ G(1+n+2e)G(1+n+e)/G(1+n)^2 ] -
        for k in range(1, NJET):        #   ln[ G(2n+2+2e)/G(2n+2) ]
            u[k] = ps1[k - 1] / kfact[k]
            v[k] = (ps1[k - 1] * (2 ** k + 1) - ps2[k - 1] * 2 ** k) / kfact[k]
        eu, ev = jexp(u), jexp(v)       # both have constant term exactly 1
        for k in range(1, NJET):        # bracket eps^0 coefficient is EXACT 0
            tot[k] += c * (eu[k] - ev[k])
        for k in range(NJET - 1):       # advance n -> n+1
            s = (-1) ** k * kfact[k]
            ps1[k] += s / mp.mpf(n + 1) ** (k + 1)
            ps2[k] += s / mp.mpf(2 * n + 2) ** (k + 1) \
                + s / mp.mpf(2 * n + 3) ** (k + 1)
        c *= mp.mpf(n + 1) ** 2 / ((2 * n + 2) * (2 * n + 3))
    # certified tail of tot (see docstring; c = c_N, ps1/ps2 = state at N)
    Au = abs(ps1[0]) + mp.fsum(abs(ps1[k - 1]) / kfact[k]
                               for k in range(2, NJET))
    Av = 3 * abs(ps1[0]) + 2 * abs(ps2[0]) + mp.fsum(
        ((2 ** k + 1) * abs(ps1[k - 1]) + 2 ** k * abs(ps2[k - 1])) / kfact[k]
        for k in range(2, NJET))
    g = mp.exp(mp.mpf(5) / (N + 1)) / 4
    tail_tot = c * (mp.exp(Au) + mp.exp(Av)) / (1 - g)
    bracket = tot[1:] + [mp.mpf(0)]     # divide by eps (exact shift)
    # pi*eps/sin(pi*eps): invert the sin(pi eps)/(pi eps) jet
    sc = [mp.mpf(0)] * NJET
    for m in range((NJET + 1) // 2):
        sc[2 * m] = (-1) ** m * mp.pi ** (2 * m) / mp.gamma(2 * m + 2)
    # 1/Gamma(1-eps) = exp(-gamma eps - sum_{k>=2} zeta(k) eps^k / k)
    lg = [mp.mpf(0)] * NJET
    lg[1] = -mp.euler
    for k in range(2, NJET):
        lg[k] = -mp.zeta(k) / k
    j_inv_sc, j_lg = jinv(sc), jexp(lg)
    s = jmul(jmul(bracket, j_inv_sc), j_lg)
    eg = [mp.mpf(0)] * NJET
    eg[1] = 2 * mp.euler                # e^{2 gamma eps}
    j_eg = jexp(eg)
    if tail_out is not None:            # convolution triangle inequality
        l1 = lambda a: mp.fsum(abs(x) for x in a)
        tail_out["jet_tail"] = tail_tot * l1(j_inv_sc) * l1(j_lg) * l1(j_eg)
        tail_out["N"] = N
    return jmul(s, j_eg)


def boundary_constants(nterms=None, tail_out=None):
    """B_j = -(sqrt3/2) [eps^j] e^{2 gamma eps} S_vac(eps),  j = 0..NJET-1."""
    jets = svac_jets(nterms, tail_out)
    if tail_out is not None:
        tail_out["B_tail"] = mp.sqrt(3) / 2 * tail_out["jet_tail"]
    return [-mp.sqrt(3) / 2 * x for x in jets]


# ------------------------- independent cross-checks -------------------------
def Lchi3(s):
    """L(chi_{-3}, s) via Hurwitz zeta (independent of the series above)."""
    return (mp.zeta(s, mp.mpf(1) / 3) - mp.zeta(s, mp.mpf(2) / 3)) / 3 ** s


def closed_controls():
    """Closed forms of B0, B1, B2 from this work (RESULT.md; r3 = e^{2 pi i/3})."""
    s3, l3 = mp.sqrt(3), mp.log(3)
    y0 = 1 - mp.exp(mp.mpc(0, 2 * mp.pi / 3))
    IL3 = mp.im(mp.polylog(3, y0))
    IL4 = mp.im(mp.polylog(4, y0))
    B0 = -mp.mpf(3) / 2 * s3 * Lchi3(2)
    B1 = -6 * IL3 - mp.pi * l3 ** 2 / 2 - mp.mpf(5) / 54 * mp.pi ** 3
    B2 = (12 * IL4 + mp.mpf(3) / 2 * s3 * Lchi3(4)
          - mp.mpf(1) / 4 * s3 * mp.pi ** 2 * Lchi3(2)
          + mp.mpf(5) / 54 * mp.pi ** 3 * l3 + mp.pi * l3 ** 3 / 6)
    return B0, B1, B2


def svac_quad(eps):
    """Form (I): S_vac(eps) by direct Bessel-moment quadrature (independent
    of the Mellin-Barnes residue derivation; used as a runtime identity check
    at numeric eps)."""
    f = lambda x: x ** (1 + eps) * mp.besselk(eps, x) ** 3
    upper = max(40, int(3.5 * mp.mp.dps))  # e^{-3x} tail < 10^{-1.3*3.5*dps}
    I = mp.quad(f, [0, 1, 5, 30, upper])
    return 2 ** (2 - eps) / mp.gamma(1 - eps) * I


def svac_scalar(eps, nterms):
    """Form (II) at a NUMERIC eps (scalar, no jets) — cross-check partner."""
    s = mp.mpf(0)
    for n in range(nterms):
        a = mp.gamma(n + 1) * mp.gamma(1 + n + eps) / mp.gamma(2 * n + 2)
        b = (mp.gamma(1 + n + 2 * eps) * mp.gamma(1 + n + eps)
             / mp.gamma(2 * n + 2 + 2 * eps))
        s += a - b
    return mp.pi / (mp.sin(mp.pi * eps) * mp.gamma(1 - eps)) * s


def agreed(a, b):
    d = abs(a - b)
    if d == 0:
        return mp.mp.dps
    return int(mp.floor(-mp.log10(d / max(abs(b), mp.mpf(1)))))


# --------------------------------- PSLQ -------------------------------------
def pslq_attempt(B, dps_hi):
    """Naming attempt in the weight-graded Gamma_1(6)/zeta/Cl2 ring
    (sixth-root-orbit dictionary: ImLi_k(1-r3), sqrt3*L(chi-3,even), pi,
    log3, odd zetas) plus generalized log-sines Ls_j^{(k)}(2pi/3).
    Protocol: positive controls (B1 at w3, B2 at w4) MUST land before the
    B3 (w5) / B4 (w6) scans are reported; two-precision stability; honest
    null otherwise.  An earlier PSLQ hunt already showed the committed dictionary
    insufficient (fits/bhunt_log.txt: n=131, 482 digits, null) — this re-runs
    the question against the measured-precision series values."""
    s3, l3, pi = mp.sqrt(3), mp.log(3), mp.pi
    y0 = 1 - mp.exp(mp.mpc(0, 2 * pi / 3))
    IL = {k: mp.im(mp.polylog(k, y0)) for k in (3, 4, 5, 6)}
    L2, L4, L6 = Lchi3(2), Lchi3(4), Lchi3(6)
    z3, z5 = mp.zeta(3), mp.zeta(5)

    def lsjk(j, k):  # Ls_j^{(k)}(2pi/3) = -int_0^{2pi/3} x^k log^{j-1-k}|2 sin(x/2)|
        f = lambda x: -x ** k * mp.log(2 * mp.sin(x / 2)) ** (j - 1 - k)
        return mp.quad(f, [0, mp.pi / 3, 2 * mp.pi / 3])

    def monos(w):  # pi^a * log3^b of total weight w
        return {f"pi^{a}l3^{b}": pi ** a * l3 ** b
                for a in range(w + 1) for b in [w - a]}

    baskets = {
        "B1(w3 control)": dict(monos(3), IL3=IL[3], z3=z3,
                               s3L2pi=s3 * L2 * pi, s3L2l3=s3 * L2 * l3),
        "B2(w4 control)": dict(monos(4), IL4=IL[4], IL3pi=IL[3] * pi,
                               IL3l3=IL[3] * l3, s3L4=s3 * L4, z3pi=z3 * pi,
                               z3l3=z3 * l3, s3L2pi2=s3 * L2 * pi ** 2,
                               s3L2pil3=s3 * L2 * pi * l3,
                               s3L2l32=s3 * L2 * l3 ** 2, L2sq=L2 ** 2),
        "B3(w5)": dict(monos(5), IL5=IL[5], IL4pi=IL[4] * pi, IL4l3=IL[4] * l3,
                       IL3pi2=IL[3] * pi ** 2, IL3pil3=IL[3] * pi * l3,
                       IL3l32=IL[3] * l3 ** 2, z5=z5, z3pi2=z3 * pi ** 2,
                       z3pil3=z3 * pi * l3, z3l32=z3 * l3 ** 2,
                       s3L4pi=s3 * L4 * pi, s3L4l3=s3 * L4 * l3,
                       s3L2pi3=s3 * L2 * pi ** 3, s3L2pi2l3=s3 * L2 * pi ** 2 * l3,
                       s3L2pil32=s3 * L2 * pi * l3 ** 2,
                       s3L2l33=s3 * L2 * l3 ** 3, s3L2z3=s3 * L2 * z3,
                       s3L2IL3=s3 * L2 * IL[3], L2sqpi=L2 ** 2 * pi,
                       L2sql3=L2 ** 2 * l3,
                       Ls5=lsjk(5, 0), Ls4pi=lsjk(4, 0) * pi,
                       Ls4_1pi=lsjk(4, 1) * pi, Ls5_1=lsjk(5, 1),
                       Ls5_2=lsjk(5, 2), Ls3pi2=lsjk(3, 0) * pi ** 2),
        "B4(w6)": dict(monos(6), IL6=IL[6], IL5pi=IL[5] * pi, IL5l3=IL[5] * l3,
                       IL4pi2=IL[4] * pi ** 2, IL4pil3=IL[4] * pi * l3,
                       IL4l32=IL[4] * l3 ** 2, IL3z3=IL[3] * z3,
                       IL3pi3=IL[3] * pi ** 3, IL3pi2l3=IL[3] * pi ** 2 * l3,
                       IL3pil32=IL[3] * pi * l3 ** 2, IL3l33=IL[3] * l3 ** 3,
                       IL3sq=IL[3] ** 2, s3L6=s3 * L6, s3L4pi2=s3 * L4 * pi ** 2,
                       s3L4pil3=s3 * L4 * pi * l3, s3L4l32=s3 * L4 * l3 ** 2,
                       L2L4=L2 * L4, s3L2z3l0=s3 * L2 * z3,  # w2+w3=w5? no: guard
                       z3sq=z3 ** 2, z5pi=z5 * pi, z5l3=z5 * l3,
                       z3pi3=z3 * pi ** 3, z3pi2l3=z3 * pi ** 2 * l3,
                       z3pil32=z3 * pi * l3 ** 2, z3l33=z3 * l3 ** 3,
                       s3L2pi4=s3 * L2 * pi ** 4, s3L2pi3l3=s3 * L2 * pi ** 3 * l3,
                       s3L2pi2l32=s3 * L2 * pi ** 2 * l3 ** 2,
                       s3L2IL4=s3 * L2 * IL[4], L2sqpi2=L2 ** 2 * pi ** 2,
                       Ls6=lsjk(6, 0), Ls5pi=lsjk(5, 0) * pi,
                       Ls6_1=lsjk(6, 1), Ls6_2=lsjk(6, 2),
                       Ls5_1pi=lsjk(5, 1) * pi, Ls4pi2=lsjk(4, 0) * pi ** 2),
    }
    # drop the mistaken-weight guard entry
    baskets["B4(w6)"].pop("s3L2z3l0")

    targets = {"B1(w3 control)": B[1], "B2(w4 control)": B[2],
               "B3(w5)": B[3], "B4(w6)": B[4]}
    dps_lo = dps_hi - 60
    controls_ok = True
    for name, basket in baskets.items():
        tgt = targets[name]
        names = sorted(basket)
        rels = {}
        for dps in (dps_lo, dps_hi):
            with mp.workdps(dps):
                vec = [mp.mpf(tgt)] + [mp.mpf(basket[n]) for n in names]
                r = mp.pslq(vec, maxcoeff=10 ** 6, maxsteps=200000)
            rels[dps] = tuple(r) if r else None
        if rels[dps_lo] and rels[dps_lo] == rels[dps_hi] and rels[dps_hi][0]:
            r = rels[dps_hi]
            combo = " + ".join(f"({-c}/{r[0]})*{n}" for c, n in
                               zip(r[1:], names) if c)
            resid = abs(sum(c * v for c, v in zip(
                r, [tgt] + [basket[n] for n in names])))
            print(f"  {name}: RELATION (two-precision stable) "
                  f"residual={mp.nstr(resid, 3)}")
            print(f"    {name.split('(')[0]} = {combo}")
        else:
            stable = "unstable-across-dps" if (rels[dps_lo] or rels[dps_hi]) \
                else "null"
            print(f"  {name}: NO relation ({stable}; n={len(names)+1}, "
                  f"dps legs {dps_lo}/{dps_hi}, maxcoeff 1e6)")
            if "control" in name:
                controls_ok = False
    if not controls_ok:
        print("  WARNING: a positive control failed — the null reports above "
              "are NOT meaningful (protocol violation).")
    return controls_ok


# --------------------------------- driver -----------------------------------
def run(dps, quad_check=False):
    mp.mp.dps = dps + 15  # guard digits; results quoted at dps
    t0 = time.perf_counter()
    # certified refine-until-bound loop (axis-3): the default nterms formula
    # is a SEED; the certified relative tail must clear 10^-(dps+TAIL_GUARD)
    # or the series is extended x1.5, RAISING (fail-closed) at CERT_CAP.
    tb, esc, nt = {}, 0, None
    tol = mp.mpf(10) ** (-(dps + TAIL_GUARD))
    while True:
        B = boundary_constants(nterms=nt, tail_out=tb)
        rel = max(tb["B_tail"] / abs(b) for b in B[:5])
        if rel < tol:
            break
        esc += 1
        if esc > CERT_CAP:
            raise RuntimeError(
                f"FAIL-CLOSED: certified series tail {mp.nstr(rel, 3)} "
                f"(relative) >= tol {mp.nstr(tol, 3)} at N={tb['N']} after "
                f"{CERT_CAP} x1.5 extensions")
        nt = int(1.5 * tb["N"]) + 20
    t_series = time.perf_counter() - t0
    print(f"series jets: {t_series:.2f} s wall (measured), "
          f"dps={dps} (+15 guard)")
    print(f"[cert] certified series-truncation tail: |dB_j|/|B_j| <= "
          f"{mp.nstr(rel, 3)} < 1e-{dps + TAIL_GUARD} -- BOUND, not estimate "
          f"(N = {tb['N']}, escalations = {esc}, fail-closed at cap "
          f"{CERT_CAP})")

    print("\npositive controls (independent closed forms, recomputed now):")
    for j, c in enumerate(closed_controls()):
        d = agreed(B[j], c)
        print(f"  B{j}: series vs closed form -> {d} digits "
              f"(floor: working precision)")
        if d < dps - 3:                 # raising, not assert (fail-closed)
            raise RuntimeError(f"CONTROL FAILURE at B{j}: {d} < {dps - 3}")

    print("\nheld-out gate (AMFlow protocol literals, comparison ONLY):")
    gate_report = {}
    gate_fails = []
    for j, lit in GATE.items():
        nd = sum(ch.isdigit() for ch in lit)
        with mp.workdps(max(mp.mp.dps, nd + 10)):
            d = agreed(B[j], mp.mpf(lit))
        cap = min(dps, nd)
        gate_report[j] = (d, nd)
        print(f"  B{j}: {d} digits vs {nd}-digit stored literal "
              f"[cap = min(dps, {nd}) = {cap}]")
        if d < cap - GATE_MARGIN:
            gate_fails.append(f"B{j}: {d} digits < bar {cap - GATE_MARGIN}")
    if gate_fails:
        print("[gate] FAIL (fail-closed): " + "; ".join(gate_fails))
        raise SystemExit(1)
    print(f"[gate] held-out table PASS: all >= min(dps, digits) - "
          f"{GATE_MARGIN}; a 1e-30 mutation of any GATE literal exits "
          "nonzero (mutation rc-gate)")

    if quad_check:
        t1 = time.perf_counter()
        with mp.workdps(50):
            e = mp.mpf(1) / 137
            q = svac_quad(e)
            s = svac_scalar(e, 120)
            dq = agreed(s, q)
        print(f"\nBessel-moment identity check (form I vs II at eps=1/137, "
              f"dps=50): {dq} digits  [{time.perf_counter()-t1:.1f} s]")
        if dq < 40:                     # raising, not assert (fail-closed)
            raise RuntimeError(
                f"form (I) vs (II) cross-check failed: {dq} < 40")

    print("\nB3 =", mp.nstr(B[3], dps))
    print("B4 =", mp.nstr(B[4], dps))
    return B, gate_report


def main():
    ap = argparse.ArgumentParser(description=__doc__.splitlines()[0])
    ap.add_argument("--dps", type=int, default=520)
    ap.add_argument("--check", action="store_true",
                    help="rerun at dps+60 and diff (two-precision rule)")
    ap.add_argument("--quad", action="store_true",
                    help="runtime Bessel-integral cross-check of the series")
    ap.add_argument("--pslq", action="store_true",
                    help="PSLQ naming attempt (controls first; honest null ok)")
    args = ap.parse_args()

    B, _ = run(args.dps, quad_check=args.quad)

    # in-run two-precision crank gate D -> D+CRANK_STEP (axis-3, raising).
    # The recompute is itself certified (refine-until-bound, fail-closed).
    mp.mp.dps = args.dps + CRANK_STEP + 15
    tb40, nt40 = {}, None
    tol40 = mp.mpf(10) ** (-(args.dps + CRANK_STEP + TAIL_GUARD))
    for _ in range(CERT_CAP + 1):
        B40 = boundary_constants(nterms=nt40, tail_out=tb40)
        if max(tb40["B_tail"] / abs(b) for b in B40[:5]) < tol40:
            break
        nt40 = int(1.5 * tb40["N"]) + 20
    else:
        raise RuntimeError(
            "FAIL-CLOSED: crank recompute series tail did not certify at "
            f"N={tb40['N']} after {CERT_CAP} x1.5 extensions")
    with mp.workdps(args.dps + CRANK_STEP + 30):
        dmin = min(agreed(B[j], B40[j]) for j in range(5))
    bar = args.dps - CRANK_MARGIN
    ok = dmin >= bar
    print(f"\n[gate] crank D/D+{CRANK_STEP} (dps {args.dps} -> "
          f"{args.dps + CRANK_STEP}): min agreement over B0..B4 = {dmin} "
          f"digits (raising bar {bar}) {'PASS' if ok else 'FAIL'}")
    if not ok:
        raise SystemExit(1)

    if args.check:
        print(f"\n--check: rerun at dps={args.dps + 60} and diff")
        B2, _ = run(args.dps + 60)
        with mp.workdps(args.dps + 75):
            for j in range(5):
                d = agreed(mp.mpf(str(B[j])), B2[j])
                print(f"  B{j}: two-precision agreement = {d} digits "
                      f"(expect ~{args.dps})")

    if args.pslq:
        print(f"\nPSLQ naming attempt (measured precision, controls first):")
        mp.mp.dps = args.dps + 15
        pslq_attempt(B, dps_hi=min(args.dps - 20, 460))


if __name__ == "__main__":
    main()