#!/usr/bin/env python3
"""LBL3E — sunrise-dressed box (3-loop light-by-light, 6 lines, all mass m):
independent evaluation of the homogeneous period
psi1(w)/pi = (2/sqrt 3) * varpi0(w), the anchor of the closed form.

varpi0 is the holomorphic Frobenius solution of the equal-mass-sunrise
Picard-Fuchs operator (t = w/m^2, m^2 = 1):

    L_sun = t(t-1)(t-9) y'' + (3t^2 - 20t + 9) y' + (t - 3) y = 0,  varpi0(0) = 1
    =>  a_{m+1} = [ (10m^2 + 10m + 3) a_m  -  m^2 a_{m-1} ] / ( 9 (m+1)^2 ),  a_0 = 1.

This script recomputes varpi0 at the Euclidean anchors via the PF series plus
Taylor-step analytic continuation (singularities only at t in {0, 1, 9}), and
compares against INDEPENDENT Gamma_1(6) q-series reference values computed in
this work (embedded below; the original gate record cross-checks the two
routes at 35.7-116.4 digits). It also checks the leading boundary constant
B^(0) = 2*Cl2(pi/3) against its recorded value (PSLQ residual 110.0 digits).

Evaluation interface (arbitrary point / arbitrary precision):
    python3 lbl3e-evaluate.py                        # gate demo + dps-doubling demo
    python3 lbl3e-evaluate.py --point=-7/3 --dps 200 # any Euclidean w, any dps
    python3 lbl3e-evaluate.py --dps 120              # demo at a different base dps
or, in code, psi1_over_pi(w, dps).  Series/step truncation orders scale with
dps (PF-operator period series: analytic integrand, no numeric node cache), so
precision is limited only by the requested dps.

DOMAIN: real Euclidean w <= 0 (m^2 = 1, t = w/m^2); the continuation path is
the negative real axis, and the singular fibres sit at w in {0, 1, 9} (cusps
of the Gamma_1(6) fibration). Minkowski w > 0 (or complex w) needs an i0+
path prescription and is NOT implemented here.
SCOPE: this script evaluates the COMPUTED pieces of the LBL3E closed form —
the homogeneous period psi1(w)/pi and the boundary constant B^(0). The
Eichler integral Eichler_f3(tau(w)), psi2, and the per-master (c1, c2)
assembly are recorded in lbl3e-expression.md (anchor data), not evaluated
here.

Runtime: seconds at default precision (measured wall time is printed).
Standalone: standard library + mpmath only (pip install mpmath).
NOTE: the gate-table agreement is capped near ~81 d by the 81-digit reference
literals embedded below; the original full-precision gate record compares the
two routes at 115.9-116.4 d on these four points. The dps-doubling demo uses
a 621-digit certified reference at w = -2, so it is NOT literal-capped at the
demo precisions (the cap is stated in the output either way).

CHANGELOG 2026-07-05 (axis3 wave2): fail-closed certified truncation on the
value path (mirror of tools/rowscripts/eval_row27.py hardening).  The
Frobenius seed series and every Taylor continuation step now certify their
truncation at runtime with the trailing-8-term-window geometric tail bound
(max |c_n h^n| x r/(1-r), r certified by the series radius / step rule); the
old nmax/nterm = f(digits) formulas survive as STARTING SEEDS only;
escalation continues the SAME recurrences x1.7 to a cap, RuntimeError at the
cap naming step/bound/tol/N/cap; per-step bounds accumulate to a value-level
certified bound printed as [certified] lines.  Default outputs are unchanged
except those new lines (healthy runs never escalate: measured headrooms
>= 1e4, calibration in <archive>/axis3_wave/row27/).
"""
import argparse
import math
import time

from mpmath import mp, mpf, log, pi, sqrt, im, polylog, exp, mpc, fabs

mp.dps = 150

# --- reference values (independent Gamma_1(6) q-series route) computed in this
# --- work, copied verbatim from its gate record (81 digits per point)
REF = {
    "-1/4": ("1.06990420417018736030436154346086918375296946817309071249078650275034810661105109", 116.4),
    "-1/2": ("1.00231667829912791054260671349038686405093521275267783741379929414027922898843651", 116.3),
    "-1":   ("0.899026529977333134189954924667982179257205514830453633179078524829944688247218676", 116.2),
    "-2":   ("0.761166273198895035051999080472538833664342316901239650920230189658510232833387943", 115.9),
}
# leading boundary constant computed in this work (PSLQ-closed as 2*Cl2(pi/3))
REF_B0 = "2.029883212819307250042405108549040571883378615060599584034978213553194952516488044272940708456513389891723655062719770803"

# psi1(-2)/pi -- INDEPENDENT certified reference for the dps-doubling demo:
# Arb ball-certified midpoint (radius ~ 2e-770) from this campaign's stage-1
# Gamma_1(6) dictionary grid (Eichler.jl, dual-route certified q-series vs
# eta-quotient; stage1/dictionaries/psi1-grid.tsv, row t = -2; the same string
# seeds the banked icecream-cone period CERT). 621 digits stored verbatim.
REF_LONG_W2 = "7.61166273198895035051999080472538833664342316901239650920230189658510232833387943389361537444200806126014140498020923008307426039511680245368880977336718158457303936680162623695006598331419827024579353294072194910881008662185051674381723382153038704057837767628769184753472831952017920667031011678942155311629692653947028811470602414686304303447033526768026055151278676051552481820050704775292609778615931047605232017455939702951279595133749716140782358379847880335875357153351779579224834855463652324033807906004225386391762193582008259491535031886204663515480211939127187899553219903248061590680035426896625844556593740e-01"

DEMO_DPS = 90  # demo target digits (sized to saturate the 81-digit gate literals)

# certified-truncation guards (axis3 wave2; calibration measured in
# <archive>/axis3_wave/row27/scratch/probe_calib.log: healthy
# per-step window bound ~1e-(digits+10), seed-series bound ~1e-(digits+24))
TAILWIN = 8          # trailing-window width (accepted axis3 pilot pattern)
TAIL_GUARD = 6       # value-series tail gate: bound < 10^-(digits+6)
TAIL_GUARD_DER = 2   # derivative-series tail gate: 10^-(digits+2)
ESCALATE = 1.7       # exact-continuation escalation factor
SEED_CAP = 16        # RuntimeError if bound not reached at 16 x seed order


def series_at_zero(t, nmax=400, digits=None):
    """varpi0 and derivative from the CERTIFIED Frobenius series at t=0.

    axis3 wave2: nmax is a STARTING SEED only.  The truncation is certified
    by the trailing-TAILWIN-window bound max|a_n t^n| * r/(1-r) with
    r = |t| <= 1/4 certified by the series radius (L_sun is regular-singular
    with singularities only at {0,1,9}, and varpi0 is the analytic Frobenius
    solution at 0, so sum a_n t^n converges on |t| < 1).  Escalation = EXACT
    continuation of the same recurrence and running sums x1.7 up to
    SEED_CAP x seed; RuntimeError at the cap.  Returns (y, yp, err_y,
    err_yp), errors gated at 10^-(digits + TAIL_GUARD[_DER]).  Caveat
    (accepted axis3 pilot pattern): the window max is the geometric envelope
    constant; subexponential factors are absorbed by the guard digits +
    escalation."""
    if digits is None:
        digits = mp.dps
    r = fabs(t)
    if r >= 1:
        # 2026-07-06 (row27-r1guard blog flag): the trailing-window bound
        # r/(1-r) is only a tail envelope for r < 1 (series radius; callers
        # clamp |t| <= 1/4).  Raise loud on direct calls outside the radius
        # rather than return a garbage bound.
        raise ValueError(
            f"series_at_zero: |t| = {mp.nstr(r, 8)} >= 1 outside the "
            f"certified series radius; transport instead of direct call")
    ncap = SEED_CAP * nmax
    fac = r / (1 - r)
    tol_v = mpf(10) ** (-(digits + TAIL_GUARD))
    tol_d = mpf(10) ** (-(digits + TAIL_GUARD_DER))
    a_prev, a_cur = mpf(0), mpf(1)  # a_{-1}, a_0
    y = a_cur
    yp = mpf(0)
    tp = mpf(1)
    m = 0
    win = [abs(a_cur)]              # |a_n t^n| trailing window
    target = nmax
    while True:
        while m < target:           # EXACT continuation of sums + recurrence
            a_next = ((10 * m * m + 10 * m + 3) * a_cur - m * m * a_prev) / (9 * (m + 1) ** 2)
            tp_next = tp * t  # t^{m+1}
            y += a_next * tp_next
            yp += (m + 1) * a_next * tp
            tp = tp_next
            a_prev, a_cur = a_cur, a_next
            m += 1
            win.append(abs(a_cur * tp))
            if len(win) > TAILWIN:
                win.pop(0)
        e_val = max(win) * fac
        e_der = (max(win) / r) * (m * fac + r / (1 - r) ** 2) if r > 0 else mpf(0)
        if e_val < tol_v and e_der < tol_d:
            break
        if target >= ncap:
            raise RuntimeError(
                f"series_at_zero: certified tail bound not reached at cap: "
                f"t = {mp.nstr(t, 10)}, bound = {mp.nstr(e_val, 3)} (tol "
                f"{mp.nstr(tol_v, 3)}), N = {m} (cap {ncap})")
        target = min(ncap, int(target * ESCALATE) + 1)
    return y, yp, e_val, e_der


def taylor_step(t0, y, yp, h, nterm, rloc, digits=None):
    """One CERTIFIED Taylor step of L_sun from t0 to t0+h, given (y, y')(t0).

    axis3 wave2: nterm is a STARTING SEED only; the truncation is certified
    by the trailing-TAILWIN-window bound max|c_n h^n| * r/(1-r) with
    r = |h|/rloc certified by the step rule (|h| <= 0.45 * rloc, rloc = the
    solution's local radius = distance to the nearest singularity in {1,9};
    varpi0 is analytic at t=0).  The derivative tail carries the exact
    n-weighted factor sum_{n>M} n r^n <= M r/(1-r) + r/(1-r)^2 on the same
    envelope.  Escalation = EXACT continuation of the same recurrence (the
    c-list is extended in place) x1.7 up to SEED_CAP x seed; RuntimeError at
    the cap.  Returns (ynew, ypnew, bound_y, bound_yp)."""
    # local polynomial coefficients of p, q, r in x = t - t0
    p = [t0**3 - 10 * t0**2 + 9 * t0, 3 * t0**2 - 20 * t0 + 9, 3 * t0 - 10, mpf(1)]
    q = [3 * t0**2 - 20 * t0 + 9, 6 * t0 - 20, mpf(3)]
    r = [t0 - 3, mpf(1)]
    if digits is None:
        digits = mp.dps
    rr = fabs(h) / rloc
    assert rr < 1
    fac = rr / (1 - rr)
    ha = fabs(h)
    tol_v = mpf(10) ** (-(digits + TAIL_GUARD))
    tol_d = mpf(10) ** (-(digits + TAIL_GUARD_DER))
    ncap = SEED_CAP * nterm
    target = nterm
    c = [y, yp]
    while True:
        for n in range(len(c) - 2, target - 2):   # EXACT continuation
            s = mpf(0)
            for j in range(1, 4):  # p_j terms, j=0 is the solved-for term
                k = n - j + 2
                if 0 <= k < len(c):
                    s += p[j] * k * (k - 1) * c[k]
            for j in range(3):
                k = n - j + 1
                if 0 <= k < len(c):
                    s += q[j] * k * c[k]
            for j in range(2):
                k = n - j
                if 0 <= k < len(c):
                    s += r[j] * c[k]
            c.append(-s / (p[0] * (n + 2) * (n + 1)))
        M = len(c) - 1
        wmax = max(abs(c[n]) * ha ** n
                   for n in range(max(0, M - (TAILWIN - 1)), M + 1))
        bound_y = wmax * fac
        bound_yp = (wmax / ha) * (M * fac + rr / (1 - rr) ** 2)
        if bound_y < tol_v and bound_yp < tol_d:
            break
        if target >= ncap:
            raise RuntimeError(
                f"taylor_step: certified tail bound not reached at cap: "
                f"t0 = {mp.nstr(t0, 10)}, h = {mp.nstr(h, 6)}, bound_y = "
                f"{mp.nstr(bound_y, 3)} (tol {mp.nstr(tol_v, 3)}), bound_yp = "
                f"{mp.nstr(bound_yp, 3)} (tol {mp.nstr(tol_d, 3)}), "
                f"N = {M + 1} (cap {ncap})")
        target = min(ncap, int(target * ESCALATE) + 1)
    ynew = mpf(0)
    ypnew = mpf(0)
    for n in range(len(c) - 1, -1, -1):
        ynew = ynew * h + c[n]
    for n in range(len(c) - 1, 0, -1):
        ypnew = ypnew * h + n * c[n]
    return ynew, ypnew, bound_y, bound_yp


def varpi0_cert(t_target, digits=None):
    """Analytic continuation along the negative real axis (sing. at 0,1,9 only).

    digits: target decimal accuracy; the series / Taylor-step truncation
    orders are sized from it (default mp.dps) as STARTING SEEDS; every
    truncation is certified at runtime (see series_at_zero / taylor_step) and
    RAISES on non-convergence.  Domain: real t_target <= 0 (Euclidean).
    Returns (varpi0, err): err = accumulated certified truncation bound
    (per-step bounds summed at value level with first-order (c0,c1) error
    transport E += |h|*Ep + bound_y, Ep += bound_yp; higher-order flow
    sensitivity is absorbed by the guard digits — the accepted axis3-pilot-
    class caveat).
    """
    if digits is None:
        digits = mp.dps
    t_target = mpf(t_target)
    if t_target > 0:
        raise ValueError(
            "Euclidean domain only: need real w <= 0 (m^2 = 1). Singular "
            "fibres at w in {0, 1, 9}; Minkowski w > 0 needs an i0+ path "
            "prescription not implemented in this demo script."
        )
    # seed with the Frobenius series at t=0 (radius 1: nearest true
    # singularity of varpi0 is t=1); keep |t0| <= 1/4 so the term ratio <= 1/4
    t0 = t_target if t_target > mpf(-1) / 4 else mpf(-1) / 4
    if t0 == 0:
        return mpf(1), mpf(0)
    nmax = int(digits / (-math.log10(abs(float(t0))))) + 40
    y, yp, E, Ep = series_at_zero(t0, nmax, digits=digits)
    eps_land = mpf(10) ** (-digits)
    while t0 > t_target + eps_land:
        dist = min(abs(t0), abs(t0 - 1), abs(t0 - 9))  # ODE-coefficient step control
        h = -min(mpf("0.45") * dist, t0 - t_target)
        # truncation order from the SOLUTION's local radius (distance to {1,9};
        # varpi0 is analytic at t=0, so t=0 limits the step size, not the tail)
        rloc = min(abs(t0 - 1), abs(t0 - 9))
        ratio = float(abs(h) / rloc)
        nterm = int(digits / (-math.log10(ratio))) + 30
        y, yp, by, byp = taylor_step(t0, y, yp, h, nterm, rloc, digits=digits)
        E, Ep = E + fabs(h) * Ep + by, Ep + byp
        t0 += h
    return y, E


def varpi0(t_target, digits=None):
    """Back-compat wrapper: certified value only (see varpi0_cert)."""
    return varpi0_cert(t_target, digits=digits)[0]


def _parse_point(w):
    """Parse a kinematic point: number, decimal string, or fraction 'p/q'."""
    if isinstance(w, str) and "/" in w:
        num, den = w.split("/")
        return mpf(num.strip()) / mpf(den.strip())
    return mpf(w)


def psi1_over_pi(w, dps=60, return_cert=False):
    """psi1(w)/pi = (2/sqrt 3) * varpi0(w/m^2) at Euclidean w <= 0, m^2 = 1.

    w:   number or string ("-2", "-7/3", "-0.37"); real, w <= 0 (see DOMAIN
         in the module docstring).
    dps: target decimal digits — arbitrary; work is done at dps+10 with
         truncation orders sized to dps+5 as certified STARTING SEEDS
         (refine-until-bound loops inside, RAISE on non-convergence).
    return_cert=True additionally returns the accumulated certified
    truncation bound (scaled by the exact prefactor 2/sqrt3).
    """
    with mp.workdps(dps + 10):
        t = _parse_point(w)
        v, err = varpi0_cert(t, digits=dps + 5)
        v = (2 / sqrt(3)) * v
        err = (2 / sqrt(3)) * err
    if return_cert:
        return v, err
    return v


def agree_digits(a, b):
    d = fabs(a - b)
    if d == 0:
        return mp.dps
    return float(-log(d / fabs(b), 10))


def run_gate_demo(base_dps=DEMO_DPS):
    print("LBL3E homogeneous-period gate reproduction (PF Frobenius vs stored q-series)")
    print("point      psi1(w)/pi (this script, PF route)        ref-agreement  recorded gate")
    worst_cert = mpf(0)
    for key, (ref, rec) in REF.items():
        v, err = psi1_over_pi(key, dps=base_dps, return_cert=True)
        worst_cert = max(worst_cert, err)
        with mp.workdps(base_dps + 10):
            a = agree_digits(v, mpf(ref))
        print(f"w={key:>5} : {mp.nstr(v, 40)}   {a:7.1f} d   ({rec} d)")
    print("(table agreement is capped by the 81-digit stored q-series literals, ~81 d cap)")
    print(f"[certified] PF-transport truncation bound <= {mp.nstr(worst_cert, 3)} "
          f"(worst gate point; per-step tol 1e-{base_dps + 5 + TAIL_GUARD}; "
          f"refine-until-bound loops, RAISE on non-convergence)")

    # --- boundary constant check: B^(0) = 2*Cl2(pi/3), Cl2(x) = Im Li2(e^{ix})
    b0 = 2 * im(polylog(2, exp(mpc(0, 1) * pi / 3)))
    print(f"\nB^(0) = 2*Cl2(pi/3) = {mp.nstr(b0, 40)}")
    print(f"recorded value (this work) = {REF_B0[:42]}...")
    print(f"agreement: {agree_digits(b0, mpf(REF_B0)):.1f} d  (recorded PSLQ residual 110.0 d)")


def run_doubling_demo(base_dps=DEMO_DPS):
    """Rerun the w=-2 gate point at 2x dps: live agreement digits must grow."""
    ref_ndig = len(REF_LONG_W2.split("e")[0].replace(".", ""))
    print(f"\ndps-doubling demo at w = -2 (vs {ref_ndig}-digit certified stage-1 reference):")
    vals = []
    certs = []
    for dps in (base_dps, 2 * base_dps):
        t_run = time.perf_counter()
        v, err = psi1_over_pi(-2, dps=dps, return_cert=True)
        dt = time.perf_counter() - t_run
        with mp.workdps(2 * base_dps + 40):
            a = agree_digits(v, mpf(REF_LONG_W2))
        vals.append(v)
        certs.append(err)
        print(f"  dps = {dps:4d} : live agreement {a:7.1f} d   ({dt:.2f} s)")
    with mp.workdps(2 * base_dps + 40):
        sc = agree_digits(vals[0], vals[1])
    print(f"  self-consistency of the two runs: {sc:.1f} d")
    print(f"  [certified] accumulated truncation bounds: {mp.nstr(certs[0], 3)} "
          f"(dps {base_dps}) / {mp.nstr(certs[1], 3)} (dps {2 * base_dps})")
    if 2 * base_dps + 5 >= ref_ndig:
        print(f"  NOTE: stored reference string is {ref_ndig} digits — agreement caps there.")
    else:
        print(f"  (stored reference string is {ref_ndig} digits; cap not reached here)")


if __name__ == "__main__":
    ap = argparse.ArgumentParser(
        description="LBL3E homogeneous period psi1(w)/pi evaluator (PF Frobenius "
        "route). Domain: real Euclidean w <= 0, m^2 = 1 (singular fibres at "
        "w in {0,1,9}; Minkowski w > 0 not implemented)."
    )
    ap.add_argument("--point", default=None,
                    help="evaluate psi1(w)/pi at this w (e.g. --point -2, "
                         "--point -0.37, or --point=-7/3 — use the '=' form "
                         "for fractions); default: gate demo + dps-doubling demo")
    ap.add_argument("--dps", type=int, default=None,
                    help=f"target decimal digits (demo default {DEMO_DPS}, "
                         "point-mode default 60)")
    args = ap.parse_args()

    t_wall = time.perf_counter()
    if args.point is None:
        base = args.dps if args.dps else DEMO_DPS
        run_gate_demo(base)
        run_doubling_demo(base)
    else:
        dps = args.dps if args.dps else 60
        try:
            v, cert = psi1_over_pi(args.point, dps=dps, return_cert=True)
        except ValueError as e:
            ap.exit(2, f"error: {e}\n")
        print(f"psi1(w)/pi at w = {args.point}  (target {dps} dps, PF Frobenius route):")
        print(f"  {mp.nstr(v, dps)}")
        print(f"  [certified] accumulated truncation bound <= {mp.nstr(cert, 3)} "
              f"(per-step tol 1e-{dps + 5 + TAIL_GUARD}; fail-closed RAISE at cap)")
        with mp.workdps(dps + 20):
            t_in = _parse_point(args.point)
            matched = False
            for key, (ref, rec) in REF.items():
                if t_in == _parse_point(key):
                    a = agree_digits(v, mpf(ref))
                    print(f"  stored 81-digit q-series ref: agreement {a:.1f} d "
                          f"(literal cap ~81 d; recorded gate {rec} d)")
                    matched = True
            if t_in == mpf(-2):
                a = agree_digits(v, mpf(REF_LONG_W2))
                print(f"  stored 621-digit certified ref: agreement {a:.1f} d "
                      f"(literal cap ~621 d)")
                matched = True
            if not matched:
                print("  (no stored reference at this point; value computed live "
                      "from the PF operator)")
    print(f"\nmeasured wall time: {time.perf_counter() - t_wall:.2f} s")