#!/usr/bin/env python3 """sudakov-evaluate.py -- the two-loop Sudakov vertex of the mpl-suite page, the planar massless ladder I[1,0,1,1,1,1,1](s): STANDALONE evaluator of its exact Gamma-function closed form and of every digit-agreement claim in the page's Sudakov section. FAMILY (two loops k1, k2; p1^2 = p2^2 = 0, (p1+p2)^2 = s; d = 4-2 eps) D1 = k1^2 D2 = (k1+p1)^2 D3 = (k1+p1+p2)^2 D4 = k2^2 D5 = (k2+p1)^2 D6 = (k2+p1+p2)^2 D7 = (k1-k2)^2 I[a_1..a_7](s) = int prod_{l=1,2} d^d k_l/(i pi^{d/2}) prod_j 1/D_j^{a_j}, Minkowski propagators with +i0, no e^{eps gamma_E} or (4 pi)^eps factor. The ladder drops D2 (six lines: the off-shell vertex (D1,D3) at the apex, the on-shell vertices (D4,D5), (D5,D6), the internal vertices (D1,D4,D7), (D3,D6,D7)); seven scalar products make D2 the family's one irreducible numerator. Values of record at s = -1; the held-out point is s = -2. CLOSED FORM (exact in d): I[1,0,1,1,1,1,1](s) = (-s)^(-2-2 eps) * [ -3/eps * G(1,2) * T(1,1+eps,1) + (1-2 eps)^2/eps^2 * G(1,1)^2 + 3(1-3 eps)/eps^2 * G(1,2) * G(1+eps,1) ] G(a,b) = Gamma(a+b-d/2) Gamma(d/2-a) Gamma(d/2-b) / (Gamma(a) Gamma(b) Gamma(d-a-b)) T(a,b,c) = Gamma(a+b+c-d/2) Gamma(d/2-a-b) Gamma(d/2-b-c) / (Gamma(a) Gamma(c) Gamma(d-a-b-c)) G is the one-loop massless bubble, T the one-loop triangle with two on-shell legs. The form is the family's own reduction onto five Gamma-product masters (Kira 3.1, at s = -1; three genuine masters -- the sunrise, the bubble-inserted triangle and the square of the one-loop bubble -- plus two dotted companions that are exact rational multiples of them). The three rational prefactors are the reduction coefficients -6/(d-4), 4(d-3)^2/(d-4)^2, 6(3d-10)/(d-4)^2 of I[1,0,1,0,1,0,2], I[1,0,1,1,0,1,0], I[0,0,1,1,0,0,2]; the masters' Minkowski signs (-1)^{sum a} are folded into the bracket. Production receipt of that reduction: sha256 f567a5b0a8d894da... Masters: Gehrmann, Huber, Maitre, hep-ph/0507061 (A_2^2, A_3, A_4); the ladder itself is not in their paper, so the agreement below is not by construction. HISTORY: the closed form shown for this row earlier in 2026 described a different object and was withdrawn; the row was re-derived from the family's own reduction, and this evaluator ships on the re-derived form above. Nothing from the withdrawn form is used. eps-EXPANSION at s = -1 (Laurent series from eps^-4). Two independent routes: (A) analytic: each Gamma(1+x eps) = exp(-gamma_E x eps + sum_{k>=2} (-1)^k zeta(k) (x eps)^k/k), so eps^4 * bracket = (3/4) E1 + E2 - (3/2) E3 with E1 = G(1+e)G(1-e)^2 G(1+2e)G(1-2e)/G(1-3e), E2 = G(1+e)^2 G(1-e)^4/G(1-2e)^2, E3 = G(1-e)^3 G(1+2e)/G(1-3e) (G = Gamma; e = eps); (B) numerical: Cauchy-circle Taylor coefficients (256-point trapezoidal rule on |eps| = 0.1, nearest Gamma pole at |eps| = 1/2) of eps^4 * bracket(eps) evaluated directly through mp.gamma. Multiplied by e^{2 eps gamma_E} the tower is pure zeta values: 1/(4 eps^4) + 0/eps^3 + 5 pi^2/(24 eps^2) + 29 zeta3/(6 eps) + 3 pi^4/32 + eps (329 zeta5/10 - 107 pi^2 zeta3/36) + ... (weight = order + 4). CHECKS (the default run; exit 0 only if ALL pass; every threshold is >= 30 agreeing digits): 1. route (A) vs route (B), and route (A) at dps vs dps+40 (two-precision rule); 2. s = -1 reference: coefficients vs the vendored auxiliary-mass-flow run (goal 70 digits); 3. HELD-OUT s = -2 reference: the vendored s = -2 run vs the closed form re-expanded through (-s)^(-2-2eps) = 2^-2 * exp(-2 eps ln 2) -- a cross-point test of the exact form; 4. zeta-value control: the e^{2 eps gamma_E}-rescaled coefficients vs the pure-zeta closed forms recovered by an integer-relation search (orders eps^-4 .. eps^5) and, at eps^6, where that search returned NULL through its whole height sequence, vs the exact coefficient of the closed form's own expansion recorded in the data file (route (A) carried out in exact arithmetic over Q[gamma_E, zeta values]; every lower order of that exact tower equals the search's closed form identically): -64793*zeta(3)*zeta(7)/21 - 63209*zeta(5)**2/50 + 46603*pi**10/22809600 + 6961*pi**4*zeta(3)**2/720 + 16907*pi**2*zeta(3)*zeta(5)/90; 5. --check: the two-precision tier -- route (A) re-run at dps + 60 and every order eps^-4 .. eps^N compared with the dps run (target >= dps - 5 digits each); verdict PASS / FAIL, folded into the exit code. Expected at the default dps 100: about 110 digits at the low orders, falling to about 78 digits at eps^6 against the numerical references (their own per-order accuracy, not the closed form's); the zeta-value control reads about 110 digits at every order including eps^6. --mutate perturbs one reduction coefficient (-3 -> -3.000000001); the run must then exit nonzero (a control that cannot fail is vacuous). DATA FILE: sudakov-data.json beside this script, sha256-pinned by DATA_SHA256 below. It vendors the three reference inputs in content -- the s = -1 run, the held-out s = -2 run (arb-style balls '[value +/- radius]', all ten integrals of the family), and the integer-relation zeta tower (with, since 2026-09-11, the exact eps^6 coefficient beside the search's NULL row at that order) -- with the sha256 of each original file recorded beside it. A mutated or truncated data file is REFUSED before any computation. Requirements: python3 + mpmath ONLY (json/hashlib/argparse stdlib). No network, no machine-path imports. mp.dps is set INSIDE main() after argparse. EXIT CODE: 0 all checks pass; 1 a check fails (this is what --mutate must produce); 2 usage error; 3 data-file pin refusal; 4 data file missing. TIMING: the default run finishes in a few seconds on a laptop-class machine; stdout is line-buffered and the first line prints immediately. CLI: python3 sudakov-evaluate.py # all checks at dps 100 python3 sudakov-evaluate.py --dps 60 # shallower (digit counts cap at dps+10) python3 sudakov-evaluate.py --point=-7/3 --eps 1/50 # + direct value at exact eps (note the '=') python3 sudakov-evaluate.py --order 8 # Laurent coefficients through eps^8 python3 sudakov-evaluate.py --check # + the two-precision tier (route A at dps+60, every order) python3 sudakov-evaluate.py --mutate # mutation control (rc != 0) """ import argparse import hashlib import json import os import re import sys import time from fractions import Fraction import mpmath as mp from mpmath import mpf sys.stdout.reconfigure(line_buffering=True) DATA_FILE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "sudakov-data.json") DATA_SHA256 = "a56ca19d2475df8cf744ce977c6e2855dfc3e1e4f290805d6d330c43951c0d76" EXIT_OK, EXIT_GATE, EXIT_USAGE, EXIT_PIN, EXIT_MISSING = 0, 1, 2, 3, 4 INDICES = [1, 0, 1, 1, 1, 1, 1] COEF = {"k1": mpf(-3), "k2": mpf(1), "k3": mpf(3)} # -3/eps, (1-2eps)^2/eps^2, 3(1-3eps)/eps^2 POLE = 4 THRESH = 30.0 class PinRefusal(RuntimeError): """Byte-level mutation of the pinned data file (exit 3).""" class MissingRefusal(RuntimeError): """Data file absent (exit 4).""" def load_data(): """sudakov-data.json, sha256-checked against DATA_SHA256 before anything is computed.""" if not os.path.isfile(DATA_FILE): raise MissingRefusal(f"{os.path.basename(DATA_FILE)} not found beside this script") raw = open(DATA_FILE, "rb").read() sha = hashlib.sha256(raw).hexdigest() if sha != DATA_SHA256: raise PinRefusal(f"sha256 of {os.path.basename(DATA_FILE)} is {sha}, pinned {DATA_SHA256}") return json.loads(raw), sha def mprat(s): fr = Fraction(s) return mpf(fr.numerator) / mpf(fr.denominator) def G(a, b, eps): d = 4 - 2 * eps return mp.gamma(a + b - d / 2) * mp.gamma(d / 2 - a) * mp.gamma(d / 2 - b) / ( mp.gamma(a) * mp.gamma(b) * mp.gamma(d - a - b)) def T(a, b, c, eps): d = 4 - 2 * eps return mp.gamma(a + b + c - d / 2) * mp.gamma(d / 2 - a - b) * mp.gamma(d / 2 - b - c) / ( mp.gamma(a) * mp.gamma(c) * mp.gamma(d - a - b - c)) def bracket(eps): k1, k2, k3 = COEF["k1"], COEF["k2"], COEF["k3"] return (k1 / eps * G(1, 2, eps) * T(1, 1 + eps, 1, eps) + k2 * (1 - 2 * eps) ** 2 / eps ** 2 * G(1, 1, eps) ** 2 + k3 * (1 - 3 * eps) / eps ** 2 * G(1, 2, eps) * G(1 + eps, 1, eps)) def ladder(s, eps): """I[1,0,1,1,1,1,1](s) at exact eps; s < 0 (Euclidean), 0 < |eps| < 1/2.""" return (-s) ** (-2 - 2 * eps) * bracket(eps) # ---- truncated power series in eps (lists of mpf, index = power) ---- def ser_exp(a, n): a = list(a) + [mpf(0)] * (n + 1 - len(a)) b = [mp.exp(a[0])] + [mpf(0)] * n for m in range(1, n + 1): b[m] = sum(k * a[k] * b[m - k] for k in range(1, m + 1)) / m return b def lngamma1(x, n): """ln Gamma(1 + x eps) through eps^n.""" return [mpf(0), -mp.euler * x] + [(-1) ** k * mp.zeta(k) * mpf(x) ** k / k for k in range(2, n + 1)] def series_analytic(nmax): """Route (A): coefficients {order: c} of the ladder at s=-1, orders -POLE..nmax.""" n = nmax + POLE L = {x: lngamma1(x, n) for x in (1, -1, 2, -2, -3)} comb = lambda w: [sum(w[x] * L[x][k] for x in w) for k in range(n + 1)] E1 = ser_exp(comb({1: 1, -1: 2, 2: 1, -2: 1, -3: -1}), n) E2 = ser_exp(comb({1: 2, -1: 4, -2: -2}), n) E3 = ser_exp(comb({-1: 3, 2: 1, -3: -1}), n) k1, k2, k3 = COEF["k1"], COEF["k2"], COEF["k3"] c = [-k1 / 4 * E1[k] + k2 * E2[k] - k3 / 2 * E3[k] for k in range(n + 1)] return {k - POLE: c[k] for k in range(n + 1)} def series_quad(nmax, radius="0.1", M=256): """Route (B): Cauchy-circle Taylor coefficients of eps^4 * bracket(eps) by the M-point trapezoidal rule (DFT) on |eps| = radius; bracket evaluated directly through mp.gamma. Aliasing error ~ (radius/0.5)^M (nearest Gamma pole at |eps| = 1/2); f(conj e) = conj f(e).""" r, n = mpf(radius), nmax + POLE w = [mp.expjpi(2 * mpf(j) / M) for j in range(M)] f = [(r * w[j]) ** POLE * bracket(r * w[j]) for j in range(M // 2 + 1)] f += [mp.conj(f[M - j]) for j in range(M // 2 + 1, M)] cs = [mp.re(sum(f[j] * w[(-k * j) % M] for j in range(M))) / M / r ** k for k in range(n + 1)] return {k - POLE: cs[k] for k in range(n + 1)} def times_exp(coef, lin, scale=1): """Laurent series coef (dict order->c) times scale * exp(lin * eps), same order range.""" lo, hi = min(coef), max(coef) ex = ser_exp([mpf(0), mpf(lin)], hi - lo) return {k: scale * sum(ex[n] * coef[k - n] for n in range(k - lo + 1)) for k in range(lo, hi + 1)} # ---- vendored references ---- def load_amflow(content): """{order: (value, radius_string)} of the ladder from a vendored auxiliary-mass-flow block.""" for r in content["result"]: if list(r["integral"]["indices"]) == INDICES: out = {} for c in r["coefficients"]: m = re.fullmatch(r"\[(\S+) \+/- (\S+)\]", c["value"]["re"]) assert m and c["value"]["im"] == "0", c out[c["order"]] = (mpf(m.group(1)), m.group(2)) return out raise KeyError(f"ladder {INDICES} not in the vendored block") def load_zeta(content): """{order: (closed-form string, mpf)} from the vendored zeta tower (pure zeta closed forms).""" ns = {"mpf": mpf, "pi2": mp.pi ** 2, "z3": mp.zeta(3), "z5": mp.zeta(5), "z7": mp.zeta(7), "z9": mp.zeta(9), "z11": mp.zeta(11)} out = {} for e in content["results"]["ladder"]: cf = e.get("closed_form_reduced") if cf is None: v = e["verdict"] if v.startswith("ZERO"): cf = "0" elif v.startswith("RATIONAL"): cf = re.search(r"->\s*\((.*)\)", v).group(1) else: continue # no closed form at this order expr = re.sub(r"(?6} | {'reference':>42} | {'closed form':>42} | {'radius':>9} | digits") for k, o, c, rad, dg in rows: print(f"{k:>6} | {mp.nstr(o, 40):>42} | {mp.nstr(c, 40):>42} | {rad:>9} | {dg:6.1f}") w = min(r[4] for r in rows) print(f" worst {w:.1f} d over orders {rows[0][0]}..{rows[-1][0]} -> {'PASS' if w >= THRESH else 'FAIL'} (>= {THRESH:.0f} d)") return w def run_checks(args, data): """Every check of the docstring at the current mp.dps. Returns (all_ok, summary line).""" nmax = max(args.order, 6) t0 = time.time() cA = series_analytic(nmax) cB = series_quad(nmax) dps0 = mp.mp.dps mp.mp.dps = dps0 + 40 cA2 = series_analytic(nmax) mp.mp.dps = dps0 print(f"\n(all digit counts below are capped by the working precision, dps+10 = {dps0})") print(f"\nLaurent coefficients at s=-1 (route A, analytic), eps^-4 .. eps^{nmax}:") for k in sorted(cA): print(f" eps^{k:<3} {mp.nstr(cA[k], min(50, args.dps)):>56} A-vs-B {digits(cB[k], cA[k]):6.1f} d" f" A(dps)-vs-A(dps+40) {digits(cA[k], cA2[k]):6.1f} d") wAB = min(digits(cB[k], cA[k]) for k in cA) wAA = min(digits(cA[k], cA2[k]) for k in cA) ok = wAB >= THRESH and wAA >= THRESH print(f" route A vs route B worst {wAB:.1f} d; two-precision worst {wAA:.1f} d -> {'PASS' if ok else 'FAIL'}") print(f" [series time {time.time() - t0:.2f} s]") src = data["sources"] r0 = load_amflow(src["amflow_s1"]["content"]) rows = [(k, r0[k][0], cA[k], r0[k][1], digits(cA[k], r0[k][0], r0[k][1])) for k in sorted(r0) if k in cA] w_r0 = table("s = -1 reference (vendored auxiliary-mass-flow run) vs closed form:", rows) ok &= w_r0 >= THRESH s2 = load_amflow(src["amflow_s2"]["content"]) pred = times_exp(cA, -2 * mp.log(2), scale=mpf(1) / 4) # (-s)^{-2-2eps} at s=-2 rows = [(k, s2[k][0], pred[k], s2[k][1], digits(pred[k], s2[k][0], s2[k][1])) for k in sorted(s2) if k in pred] w_s2 = table("HELD-OUT s = -2 reference (vendored run) vs closed form x 2^(-2-2eps) re-expanded:", rows) ok &= w_s2 >= THRESH mz = load_zeta(src["zeta_tower"]["content"]) resc = times_exp(cA, 2 * mp.euler) print("\nzeta-value control: e^(2 eps gamma_E) x coefficients vs the integer-relation closed forms:") print(f"{'order':>6} | {'closed form (pure zeta)':>52} | {'rescaled coefficient':>42} | digits") w_mz = 999.0 for k in sorted(resc): if k in mz: dg = digits(resc[k], mz[k][1]) w_mz = min(w_mz, dg) print(f"{k:>6} | {mz[k][0]:>52} | {mp.nstr(resc[k], 40):>42} | {dg:6.1f}") else: print(f"{k:>6} | {'(no closed form recovered at this order)':>52} | {mp.nstr(resc[k], 40):>42} | n/a") print(f" worst {w_mz:.1f} d -> {'PASS' if w_mz >= THRESH else 'FAIL'} (>= {THRESH:.0f} d)") ok &= w_mz >= THRESH summary = (f"SUMMARY A-vs-B {wAB:.1f} d | two-precision {wAA:.1f} d | s=-1 reference {w_r0:.1f} d | " f"held-out s=-2 reference {w_s2:.1f} d | zeta control {w_mz:.1f} d" f"{' | MUTATED k1=' + mp.nstr(COEF['k1'], 12) if args.mutate else ''}") return bool(ok), summary, cA def main(): ap = argparse.ArgumentParser(description="two-loop Sudakov vertex: planar massless ladder I[1,0,1,1,1,1,1], exact Gamma closed form") ap.add_argument("--dps", type=int, default=100, help="reported precision (working dps = this + 10); >= 40") ap.add_argument("--point", action="append", default=[], help="Euclidean s < 0 (rational string, write --point=-7/3); repeatable") ap.add_argument("--eps", default="1/100", help="eps value for --point (rational string, 0 < |eps| < 1/2)") ap.add_argument("--order", type=int, default=6, help="highest eps order of the Laurent series (>= 6 is always computed)") ap.add_argument("--check", action="store_true", help="two-precision tier: route (A) re-run at dps+60, every order compared (>= dps-5 digits)") ap.add_argument("--mutate", action="store_true", help="control: perturb the -3 reduction coefficient to -3.000000001; run must exit nonzero") args = ap.parse_args() if args.dps < 40: ap.error("--dps must be >= 40 (the checks require >= 30 agreeing digits)") print("== two-loop Sudakov vertex: planar massless ladder I[1,0,1,1,1,1,1](s) = (-s)^(-2-2eps) [ -3/eps G(1,2) T(1,1+eps,1)" " + (1-2eps)^2/eps^2 G(1,1)^2 + 3(1-3eps)/eps^2 G(1,2) G(1+eps,1) ] ==") T0 = time.time() data, sha = load_data() # pin check FIRST: exit 3/4 before any computation src = data["sources"] print(f"[data] {os.path.basename(DATA_FILE)} sha256 {sha[:16]}... matches the pin; vendored: " f"{src['amflow_s1']['original_file']} ({src['amflow_s1']['original_sha256'][:16]}...), " f"{src['amflow_s2']['original_file']} ({src['amflow_s2']['original_sha256'][:16]}...), " f"{src['zeta_tower']['original_file']} ({src['zeta_tower']['original_sha256'][:16]}...)") mp.mp.dps = args.dps + 10 if args.mutate: COEF["k1"] = mpf("-3.000000001") print("[mutate] reduction coefficient -3 -> -3.000000001 (this run MUST exit nonzero)") print(f"dps={args.dps} (+10 guard); d=4-2eps; s<0; series from eps^-4") ok, summary, cA = run_checks(args, data) if args.point: eps = mprat(args.eps) if not 0 < abs(eps) < mpf(1) / 2: ap.error("--eps must satisfy 0 < |eps| < 1/2") print(f"\ndirect values at eps={args.eps} (series check truncated at eps^{max(cA)}):") for p in args.point: s = mprat(p) if not s < 0: ap.error(f"--point={p}: Euclidean s < 0 only") v = ladder(s, eps) ser = (-s) ** (-2 - 2 * eps) * sum(cA[k] * eps ** k for k in cA) print(f" s={p:>8}: I = {mp.nstr(v, min(40, args.dps))} series {digits(ser, v):.1f} d") if args.check: hi = args.dps + 60 print(f"\n--check: rerunning route (A) at dps {hi} ...") mp.mp.dps = hi + 10 cA_hi = series_analytic(max(args.order, 6)) mp.mp.dps = args.dps + 10 tgt = args.dps - 5 ok_check = True for k in sorted(cA): d = digits(cA[k], cA_hi[k]) ok_check = ok_check and d >= tgt print(f" eps^{k:<3} stable to {d:.1f} d (target >= {tgt})") print(f"--check verdict: {'PASS' if ok_check else 'FAIL'}") if not ok_check: ok = False summary += " | --check FAIL" else: summary += " | --check PASS" print(f"\n{summary}") print(f"[wall] {time.time() - T0:.2f}s at dps {args.dps} (+10 working guard)") if not ok: print("OVERALL FAIL: a check fell below 30 agreeing digits (see the FAIL lines above)") return EXIT_GATE print("OVERALL PASS") return EXIT_OK if __name__ == "__main__": try: sys.exit(main()) except PinRefusal as e: print(f"\nREFUSED (data pin, exit {EXIT_PIN}): {e}") sys.exit(EXIT_PIN) except MissingRefusal as e: print(f"\nREFUSED (missing data, exit {EXIT_MISSING}): {e}") sys.exit(EXIT_MISSING)