#!/usr/bin/env python3
"""c5-ladder-evaluate.py -- the three-loop C5 massless ladder of the mpl-suite page: STANDALONE.

Closed form (Usyukina-Davydychev / Isaev L=3 tower, as printed in the paper and on the page):

    f_3(x) = 120 Li_6(-x) - 60 ln x Li_5(-x) + 12 ln^2 x Li_4(-x) - ln^3 x Li_3(-x)

with x = p_1^2/(p_3^2 - p_1^2) the single ratio of the off-shell three-point ladder, X = x/(1+x)
(the paper's Sec 2.4); coefficients c^(3) = {120,-60,12,-1} are the
exact UD c^(L)_k = (-1)^k (2L-k)!/(k! (L-k)!) at L=3. All coefficients are exact integers --
NO interim constants anywhere in this row.

DOMAIN: x real, x > 0 (physical Euclidean ratio; f_3 is real there). Singular locus
x in {0, -1, infinity}. Complex x off the cut (-inf,-1] also evaluates (mpmath polylog),
but only x>0 is the documented/gated domain.

REFERENCES (held-out literals printed in the paper's "Independent check"
paragraph -- the paper's own held-out points x in {11/20, 7/25, 1/50},
never in the original fit set):
    f_3(11/20) : 40-digit literal printed in the paper
    f_3(7/25), f_3(1/50) : 13-digit literals printed in the paper
The paper's full ~440-digit held-out residual files are not vendored here; beyond the
printed 40/13 digits this script SELF-CHECKS via (a) --check two-precision rule and
(b) an independent integral-representation control (mp.quad of
Li_n(z) = z/Gamma(n) * int_0^1 ln(1/t)^{n-1}/(1-z t) dt), run at every gate point.

Requirements: python3 + mpmath ONLY. No AMFlow/Kira/network/file reads. mp.dps is set
INSIDE main() after argparse (module-level mpf would silently truncate).

EXIT CODE: 0 only when f_3(11/20) agrees with the paper literal to >= 40
digits, both 13-digit literals to >= 12, and the quadrature control to >= 40
digits (and, with --check, the two-precision rerun is stable); nonzero
otherwise. --mutate perturbs the 40-digit reference literal and must exit
nonzero -- a control that cannot fail is vacuous.

CLI:
  python3 c5-ladder-evaluate.py                          # reference table + control
  python3 c5-ladder-evaluate.py --dps 100 --point 3/7    # f_3 at arbitrary x>0
  python3 c5-ladder-evaluate.py --check                  # rerun at dps+60, diff
  python3 c5-ladder-evaluate.py --mutate                 # mutation control (rc != 0)
"""

import argparse
import sys
import time
from fractions import Fraction

import mpmath as mp

sys.stdout.reconfigure(line_buffering=True)


def mprat(s):
    fr = Fraction(s)
    return mp.mpf(fr.numerator) / mp.mpf(fr.denominator)

# Held-out reference literals -- transcribed from the paper (2026-07-03).
GATE_REFS = {
    "11/20": "-87.26509162505445234451646574652768104466",
    "7/25":  "-60.57473726123",
    "1/50":  "-11.95353137881",
}

C3 = [120, -60, 12, -1]  # exact UD coefficients c^(3)_k, k=0..3


def f3(x):
    """f_3(x) via mpmath.polylog (primary evaluator)."""
    x = mp.mpf(x) if mp.im(mp.mpc(x)) == 0 else mp.mpc(x)
    L = mp.log(x)
    return sum(C3[k] * L**k * mp.polylog(6 - k, -x) for k in range(4))


def li_quad(n, z):
    """Independent control: Li_n(z) = z/Gamma(n) int_0^1 ln(1/t)^{n-1}/(1-z t) dt."""
    return z / mp.gamma(n) * mp.quad(lambda t: mp.log(1 / t)**(n - 1) / (1 - z * t), [0, 1])


def f3_quad(x):
    """f_3 built entirely from the quadrature Li_n control stack."""
    x = mp.mpf(x)
    L = mp.log(x)
    return sum(C3[k] * L**k * li_quad(6 - k, -x) for k in range(4))


def digits_agree(a, b):
    a, b = mp.mpf(a), mp.mpf(b)
    if a == b:
        return float(mp.mp.dps)
    return float(-mp.log10(abs((a - b) / b)))


def run(points):
    out = {}
    for p in points:
        x = mprat((p))
        assert x > 0, "documented domain is x > 0"
        out[p] = f3(x)
    return out


def main():
    ap = argparse.ArgumentParser(description="Row 04: three-loop C5 ladder f_3(x), UD closed form")
    ap.add_argument("--dps", type=int, default=60)
    ap.add_argument("--point", action="append", default=[],
                    help="extra x > 0 (rational string, e.g. 3/7); repeatable")
    ap.add_argument("--check", action="store_true", help="two-precision rule: rerun at dps+60")
    ap.add_argument("--mutate", action="store_true",
                    help="perturb the 40-digit reference literal; run must exit nonzero")
    ap.add_argument("--control-dps", type=int, default=40,
                    help="dps for the independent quadrature control (default 40)")
    args = ap.parse_args()

    mp.mp.dps = args.dps + 10  # guard digits
    t0 = time.time()
    fails = []
    if args.mutate:
        s = GATE_REFS["11/20"]
        GATE_REFS["11/20"] = s[:20] + str((int(s[20]) + 1) % 10) + s[21:]
        print("[mutate] 40-digit reference literal perturbed at digit 20 "
              "(this run MUST exit nonzero)")

    pts = list(GATE_REFS.keys()) + [p for p in args.point if p not in GATE_REFS]
    print("c5-ladder  three-loop C5 ladder  f_3(x) = 120 Li6(-x) - 60 lnx Li5(-x) "
          "+ 12 ln^2x Li4(-x) - ln^3x Li3(-x)", flush=True)
    print(f"dps={args.dps} (+10 guard); domain x>0; references = the paper's held-out "
          "literals (40d / 13d printed digits)\n")

    vals = run(pts)
    print(f"{'x':>8} | {'f_3(x) (computed now)':>44} | reference")
    for p, v in vals.items():
        if p in GATE_REFS:
            d = digits_agree(v, mp.mpf(GATE_REFS[p]))
            nref = len(GATE_REFS[p].replace("-", "").replace(".", ""))
            tgt = 40.0 if nref >= 40 else 12.0
            if d < tgt:
                fails.append(f"f_3({p}) {d:.1f}d < {tgt:.0f}d")
            print(f"{p:>8} | {mp.nstr(v, 40):>44} | {d:6.1f} d vs {nref}-digit paper literal"
                  f" (need >= {tgt:.0f})")
        else:
            print(f"{p:>8} | {mp.nstr(v, 40):>44} | (new point, no reference)")

    # Independent integral-representation control (self-gate component)
    t1 = time.time()
    saved = mp.mp.dps
    mp.mp.dps = args.control_dps + 10
    print(f"\nindependent control (Li_n via mp.quad integral rep, dps={args.control_dps}):")
    for p in GATE_REFS:
        a = f3(mprat((p)))
        b = f3_quad(mprat((p)))
        dq = digits_agree(a, b)
        if dq < 40.0:
            fails.append(f"quadrature control x={p} {dq:.1f}d < 40d")
        print(f"  x={p:>6}: polylog vs quadrature agree to {dq:.1f} d "
              f"(need >= 40)")
    mp.mp.dps = saved
    t_ctrl = time.time() - t1

    if args.check:
        print(f"\n--check: rerunning at dps {args.dps + 60} ...")
        mp.mp.dps = args.dps + 60 + 10
        vals2 = run(pts)
        ok = True
        for p in pts:
            d = digits_agree(vals[p], vals2[p])
            tgt = args.dps - 5
            ok = ok and d >= tgt
            print(f"    x={p}: stable to {d:.1f} d (target >= {tgt})")
        print(f"--check verdict: {'PASS' if ok else 'FAIL'}")
        if not ok:
            fails.append("--check two-precision rerun unstable")

    print(f"\ntimings [s]: control {t_ctrl:.2f}, TOTAL {time.time() - t0:.2f}")
    if fails:
        print(f"OVERALL FAIL: {'; '.join(fails)}")
        return 1
    print("OVERALL PASS")
    return 0


if __name__ == "__main__":
    sys.exit(main())