#!/usr/bin/env python3 """massless-dbox-evaluate.py — the massless planar double box of the mpl-suite page (Smirnov hep-ph/9905323): standalone evaluator of the uniform-weight HPL/GPL tower and its checks. Object (the results paper, Sec. 2.3): I_row(s,t;eps) = LS * g(x,eps), LS = 1/(-2 s^2 t), x = t/s, g(x,eps) = sum_{j=-4}^{0} eps^j g_{4+j}(x), g_w pure weight w, g_w = sum c_a * zeta_k * G(a_1..a_m; x), letters a_i in {0,-1} (alphabet {x, 1+x}), exact rational c_a vendored in row02_data.json. Normalization (documented, derived once in the build notes): g = -K/4 with Smirnov's K(x,eps) ==> g = 1/eps^4 + O(eps^-3), c_empty=1; DBox_raw = -8 (i pi^{d/2} e^{-gammaE eps})^2 (-s)^{-2eps} * LS * g. Coefficient provenance: derived by exact Fraction shuffle algebra from the closed form transcribed from the arXiv TeX source of hep-ph/9905323 (eqs. (2box), (K0t), the final K_{1t} block; S_{1,2}=Li_{2,1}, S_{2,2}=Li_{3,1}), by a one-time derivation script. NOTE: the TeX source has "-(1/2)ln^4 y" COMMENTED OUT in K_{0s}; the probe (probe_sform_mb.py) confirms the printed (excluded) form is the correct one (including it breaks the x>1 match at 2e-3; excluding it matches the tower to 1e-61). What runs at runtime (mpmath only, no network, no machine-path imports): 1. TOWER: Goncharov transport evaluator (Taylor-stepped iterated-integral ODE d/dx G(a,w;x) = G(w;x)/(x-a) from x=0, SAFETY=0.25 step ratio, Taylor order 1.7*dps+30 as the starting guess) + the vendored exact rational coefficients. Every step carries a CERTIFIED geometric tail bound (ratio 1/4 by the step rule); N is grown (up to 8x) until the bound is < 10^-(dps+8), else RuntimeError. The accumulated bound is reported with the values (a bound, not an estimate). 2. ORACLE (independent polylog stack): Smirnov's closed form evaluated with mpmath.polylog + direct nested sums for Li_{2,1}, Li_{3,1}; x<1 uses the K_{0t}+K_{1t} form, x>1 the K_{0s}+K_{1s} (y=1/x) form. AMFlow ran directly on this family: sixteen staggered goal-40/70 pairs; two Euclidean points held out of the fit gate the served tower (item 5 below). Tower and oracle share the same analytic source but are evaluated on independent code paths/function stacks; the independent numerics are item 5. 3. CONTROL (independent numerical integration): Smirnov's one-fold Mellin-Barnes representation of K_{1t} (eq. (K1tMB)), integrated numerically with mp.quad on Re w = -1/2 — checks the nontrivial eps^-1 and eps^0 content against the tower with no polylog evaluation at all. 4. Positive control: transported G(0,-1;x) vs -mpmath.polylog(2,-x). 5. AMFLOW HELD-OUT GATE (independent numerics): the tower at x = 2/5 and x = 3/4 -- the two Euclidean points held out of the fit that produced this tower (ten fit points, x in {1/7, ..., 8/9}; both lists and the source lines in the data file's references.amflow_heldout block) -- mapped g -> J by J_j = (1/x) [e^(-2 gamma_E eps) K]_j with K = -4 g at s = -1, t = -x (the AMFlow measure d^d l/(i pi^(d/2)); the map's constants are read from that block, not typed here), compared per order eps^-4..eps^0 with BOTH the goal-40 and the goal-70 raw output strings of the AMFlow runs of record (vendored verbatim in the data file, arb ball form '[midpoint +/- radius]', about 110 printed digits). PASS iff at every order and both points the relative agreement with both strings is >= bar = min(pair-certified digits of that order, working dps - 5), the pair-certified digits being the whole digits on which the goal-40 and goal-70 strings agree (a digit counts only when both runs agree on it; capped by the string length) as counted in the record and vendored per order in the data file; the same count is recomputed from the shipped strings at run time and must agree with the vendored figure (so a changed string cannot lower its own bar). At the default dps the bar is the tower's own precision (55); at --dps >= 100 the bar is the pair's, at least 89 digits at every order. Only the tower is evaluated at the two points (the same code path as --point); the oracle and the control are not rerun there. Domain: Euclidean planar region s<0, t<0 ==> x = t/s > 0, arbitrary (transport cost is logarithmic in x). Oracle nested sums converge geometrically in min(x,1/x); for x in (0.9, 1.1) the Nielsen pieces fall back to the mp.quad integral representation (soft precision cap ~ 100d there). EXIT CODE: 0 only when every eps order agrees with the Smirnov oracle to within 5 digits of the working dps at every point, (unless --no-control) the Mellin-Barnes control reaches >= 40 digits at eps^-1 and eps^0, and the AMFLOW HELD-OUT GATE passes at both points and every order (and, with --check, the dps+60 rerun is stable and the gate passes again there, its bar min(pair-certified digits, dps+55)); 1 otherwise. 3 when the data file does not match its sha256 pin (refused by name before any computation); 4 when the data file is missing; 2 on a usage error. --mutate perturbs one vendored tower coefficient and must exit nonzero (the oracle comparison and the gate both fail by name). Usage: python3 massless-dbox-evaluate.py # checks at x=3/11 and x=9/7, dps 50 python3 massless-dbox-evaluate.py --dps 150 # deep run python3 massless-dbox-evaluate.py --point 5/9 # extra point(s) python3 massless-dbox-evaluate.py --check # rerun at dps+60 and diff python3 massless-dbox-evaluate.py --no-control # skip the MB control (faster) python3 massless-dbox-evaluate.py --mutate # mutation control (rc != 0) python3 massless-dbox-evaluate.py --dps 110 # the gate's bar becomes the pair's (the record's own count) Measured walls (one process, nice 10, a shared 96-core host under load; the delivered bytes): default 30.45 s; --check 35.21 s; --point 2/5 31.71 s; --no-control 3.14 s; --dps 110 31.42 s; --dps 150 28.24 s; --mutate 10.93 s. The two gate points' tower evaluations: 0.37 s together at working dps 60 (the default) and 0.26 s at working dps 120 (the --check rerun). """ import argparse import hashlib import json import os import sys import time from fractions import Fraction from mpmath import (mp, mpf, mpc, log, polylog, zeta, pi, gamma, psi, quad, euler, factorial) sys.stdout.reconfigure(line_buffering=True) DATA_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "massless-dbox-data.json") GATE_POINTS = ("3/11", "9/7") # held out of the derivation solve/verify grid # sha256 of massless-dbox-data.json (script-emitted from the shipped bytes); a # mismatch is refused by name before any computation (rc 3), a missing file rc 4 DATA_SHA256 = "93a42d0da28c6f0da8ca769d30ce677050dadbee42b3874d08d2c929a284a489" FAILS = [] # --------------------------------------------------------------------- tower def load_data(): if not os.path.exists(DATA_PATH): print(f"REFUSED (rc 4): the data file {os.path.basename(DATA_PATH)} is " f"missing beside the script; nothing computed") sys.exit(4) with open(DATA_PATH, "rb") as f: raw = f.read() got = hashlib.sha256(raw).hexdigest() if got != DATA_SHA256: print(f"REFUSED (rc 3): {os.path.basename(DATA_PATH)} sha256 {got[:16]}... " f"does not match the pin {DATA_SHA256[:16]}... (the vendored tower " f"coefficients or the AMFlow reference strings changed); nothing " f"computed") sys.exit(3) data = json.loads(raw.decode()) orders = {} for j, blk in data["orders"].items(): terms = [] for c, zk, p, ws in blk["reduced"]: word = tuple(0 if a == "0" else 1 for a in ws.split(",")) if ws else () terms.append((Fraction(c), zk, p, word)) orders[int(j)] = terms return data, orders def transport(words, x): """G(u;x) for words u over letters (0,-1) with nonzero trailing letter, Taylor-stepped from x=0 where they all vanish. x > 0. CERTIFIED TAIL BOUND: the step rule (h = dmin/4) guarantees series ratio r <= 1/4 for every kernel, so the truncation tail of each step obeys tail <= max(|c_n h^n| over the last TAILWIN computed terms) * r/(1-r). The bound is checked against 10^-(dps+GUARD); N = 1.7*dps+30 stays the starting guess (fast path unchanged when it suffices) and is DOUBLED — extra terms by exact continuation of the same recurrences — until the bound passes, up to a hard cap of 8x the initial N; hitting the cap unconverged raises RuntimeError naming the step, x-position and achieved bound. Per-step bounds accumulate. Returns (vals, err_bound): err_bound is a certified BOUND (not an estimate) on the absolute truncation error of every returned G value.""" LET = (mpf(0), mpf(-1)) need = set() for w in words: for i in range(len(w)): need.add(w[i:]) order = sorted(need, key=len) N0 = int(1.7 * mp.dps) + 30 # starting guess; SAFETY=0.25 -> r=1/4 NCAP = 8 * N0 # hard cap on Taylor order GUARD = 8 TAILWIN = 8 tol = mpf(10) ** (-(mp.dps + GUARD)) rfac = mpf(1) / 3 # r/(1-r) at certified ratio r = 1/4 vals = {w: mpf(0) for w in order} err_total = mpf(0) x0 = mpf(0) xt = mpf(x) step = 0 while x0 < xt: step += 1 h = (mpf(1) if x0 == 0 else min(x0, x0 + 1)) / 4 if xt <= x0 + h: h = xt - x0 N = N0 C = {(): [mpf(1)] + [mpf(0)] * N} for w in order: a = LET[w[0]] rest = C[w[1:]] c = [mpf(0)] * (N + 1) c[0] = vals[w] if a == 0 and x0 == 0: for n in range(1, N + 1): c[n] = rest[n] / n else: d = x0 - a 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[w] = c while True: # certified geometric tail: ratio <= 1/4 by the step rule bound = mpf(0) hn0 = h ** (N - TAILWIN + 1) for w in order: c = C[w] hn = hn0 t = mpf(0) for n in range(N - TAILWIN + 1, N + 1): tn = abs(c[n]) * hn if tn > t: t = tn hn = hn * h if t > bound: bound = t bound = bound * rfac if bound < tol: break if N >= NCAP: raise RuntimeError( "row02 transport: certified tail bound NOT met at " f"step {step}, x0 = {mp.nstr(x0, 20)}, h = {mp.nstr(h, 10)}" f": achieved bound {mp.nstr(bound, 5)} >= " f"tol {mp.nstr(tol, 5)} at Taylor order N = {N} " f"(cap {NCAP} = 8 x initial {N0})") N2 = min(2 * N, NCAP) # extend all series by exact continuation of the recurrences C[()].extend([mpf(0)] * (N2 - N)) for w in order: # suffix-closure order: rest first a = LET[w[0]] rest = C[w[1:]] c = C[w] c.extend([mpf(0)] * (N2 - N)) if a == 0 and x0 == 0: for n in range(N + 1, N2 + 1): c[n] = rest[n] / n else: d = x0 - a p0, r = 1 / d, -1 / d s = c[N] * N # recurrence state: c[n] = s/n for n in range(N + 1, N2 + 1): s = r * s + p0 * rest[n - 1] c[n] = s / n N = N2 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 err_total += bound x0 += h return vals, err_total def tower_orders(orders, x): """g_j(x), j=-4..0, from the vendored exact rational coefficients. Returns (out, tv, certs): certs[j] is a certified bound on the absolute error of out[j] from transport truncation (per-G bound from the transport gate times the l1-norm of the exact prefactors).""" x = mpf(x) words = [w for terms in orders.values() for (_, _, _, w) in terms if w] tv, terr = transport(words, x) L = log(x) zval = {0: mpf(1), 2: zeta(2), 3: zeta(3), 4: zeta(4)} out = {} certs = {} for j, terms in orders.items(): v = mpf(0) b = mpf(0) for c, zk, p, u in terms: gv = tv[u] if u else mpf(1) coef = mpf(c.numerator) / c.denominator * zval[zk] * L ** p v += coef * gv if u: b += abs(coef) * terr out[j] = v certs[j] = b return out, tv, certs # -------------------------------------------------------------------- oracle def nested_li(ms, z): """Li_{m1..mk}(z) = sum_{n1>..>nk>=1} z^{n1}/prod n_i^{m_i}, |z|<1.""" tol = mpf(10) ** (-(mp.dps - 5)) N = 64 while True: B = [mpf(0)] * (N + 1) run = mpf(0) for n in range(1, N + 1): B[n] = run run += mpf(1) / mpf(n) ** ms[-1] for j in range(len(ms) - 1, 1, -1): C = [mpf(0)] * (N + 1) run = mpf(0) for n in range(1, N + 1): C[n] = run run += (mpf(1) / mpf(n) ** ms[j - 1]) * B[n] B = C s = mpf(0) xp = mpf(1) for n in range(1, N + 1): xp *= z s += xp / mpf(n) ** ms[0] * (B[n] if len(ms) > 1 else 1) tail = abs(xp) / (1 - abs(z)) if tail < tol or N > 400000: return s N = int(N * 1.8) + 32 def nielsen_quad(n, p, z): """Nielsen S_{n,p}(z) by the integral rep (Smirnov eq. (Sab)); fallback for |z| ~ 1 where the nested sums converge too slowly.""" z = mpf(z) pref = mpf(-1) ** (n - 1 + p) / (gamma(n) * gamma(p + 1)) f = lambda t: log(t) ** (n - 1) * log(1 - z * t) ** p / t return pref * quad(f, [0, 1]) def _s12_s22(z, x_ratio): """S_{1,2}(z)=Li_{2,1}(z), S_{2,2}(z)=Li_{3,1}(z); z=-x_ratio, 01: K_{0s}+K_{1s} (y=1/x). g = -K/4. Independent stack: mpmath.polylog + nested sums / quad.""" x = mpf(x) z3 = zeta(3) P = pi ** 2 if x <= 1: L = log(x) L1 = log(1 + x) li2, li3, li4 = polylog(2, -x), polylog(3, -x), polylog(4, -x) s12, s22 = _s12_s22(-x, x) K = {-4: mpf(-4), -3: 5 * L, -2: -(2 * L ** 2 - mpf(5) / 2 * P), -1: -(mpf(2) / 3 * L ** 3 + mpf(11) / 2 * P * L - mpf(65) / 3 * z3), 0: (mpf(4) / 3 * L ** 4 + 6 * P * L ** 2 - mpf(88) / 3 * z3 * L + mpf(29) / 30 * pi ** 4)} K[-1] += -2 * (2 * li3 - 2 * L * li2 - (L ** 2 + P) * L1) K[0] += (-4 * (s22 - L * s12) + 44 * li4 - 4 * (L1 + 6 * L) * li3 + 2 * (L ** 2 + 2 * L * L1 + mpf(10) / 3 * P) * li2 + (L ** 2 + P) * L1 ** 2 - mpf(2) / 3 * (4 * L ** 3 + 5 * P * L - 6 * z3) * L1) else: y = 1 / x L = log(y) L1 = log(1 + y) li2, li3, li4 = polylog(2, -y), polylog(3, -y), polylog(4, -y) s12, s22 = _s12_s22(-y, y) K = {-4: mpf(-4), -3: -5 * L, -2: -(2 * L ** 2 - mpf(5) / 2 * P), -1: mpf(7) / 2 * P * L + mpf(65) / 3 * z3, 0: (mpf(1) / 3 * P * L ** 2 + mpf(76) / 3 * z3 * L - mpf(83) / 90 * pi ** 4)} K[-1] += -2 * (2 * li3 - 2 * L * li2 - (L ** 2 + P) * L1) K[0] += (-4 * (s22 - L * s12) - 36 * li4 - 4 * (L1 - 5 * L) * li3 - 2 * (L ** 2 - 2 * L * L1 + mpf(10) / 3 * P) * li2 + (L ** 2 + P) * L1 ** 2 + 2 * (L ** 3 + mpf(2) / 3 * P * L + 2 * z3) * L1) return {j: -K[j] / 4 for j in range(-4, 1)} # ------------------------------------------------------------------- control def mb_k1t(x): """Smirnov eq. (K1tMB): one-fold MB for K_{1t}, numerically on Re w = -1/2. Returns (eps^-1 coeff, eps^0 coeff). Valid for any x>0.""" x = mpf(x) def f(t, mode): w = mpc(mpf(-1) / 2, t) B = x ** (w + 1) / (1 + w) * gamma(1 + w) ** 3 * gamma(-w) ** 3 if mode == 0: return B.real return (B * (-5 / (1 + w) + 3 * psi(0, 1 + w) - 4 * psi(0, -w) - euler)).real m1 = (4 / pi) * quad(lambda t: f(t, 0), [0, 8, 30]) z0 = (4 / pi) * quad(lambda t: f(t, 1), [0, 8, 30]) return m1, z0 def k0t_logs(x): """Elementary K_{0t} pieces at eps^-1, eps^0 (logs and constants only).""" L = log(mpf(x)) P = pi ** 2 z3 = zeta(3) return {-1: -(mpf(2) / 3 * L ** 3 + mpf(11) / 2 * P * L - mpf(65) / 3 * z3), 0: (mpf(4) / 3 * L ** 4 + 6 * P * L ** 2 - mpf(88) / 3 * z3 * L + mpf(29) / 30 * pi ** 4)} # ------------------------------------------------- 5. AMFLOW HELD-OUT GATE def ball_mid(s): """midpoint string of an arb-style ball '[mid +/- rad]' (a plain number string is returned as is).""" s = s.strip() if s.startswith("["): s = s[1:-1].split("+/-")[0].strip() return s def sig_digits(mid): """significant decimal digits carried by a printed midpoint string.""" return len(mid.lstrip("-").replace(".", "").lstrip("0")) def rel_digits(value, ref_mid): """relative agreement -log10(|value - ref| / |ref|) of a value with a printed reference string, capped by the string's own significant digits (the digit law of the record's count); evaluated at a precision that resolves every printed digit of the string.""" cap = sig_digits(ref_mid) with mp.workdps(max(mp.dps, cap + 10)): ref = mpf(ref_mid) d = abs(value - ref) if d == 0: return float(cap) return min(float(-log(d / abs(ref)) / log(10)), float(cap)) def j_from_g(g, x, norm): """g_j (j=-4..0) -> J_j in the AMFlow normalisation at (s,t) = (-1,-x): K = kfac * g, J_j = (1/x) sum_{m>=0} (cgam*gamma_E)^m/m! K_{j-m}; kfac and cgam are read from the data file's references block (K_from_g_factor -4, exponent_gamma_coefficient -2), never typed here.""" kfac = Fraction(norm["K_from_g_factor"]) cgam = Fraction(norm["exponent_gamma_coefficient"]) kf = mpf(kfac.numerator) / kfac.denominator c = mpf(cgam.numerator) / cgam.denominator * euler K = {j: kf * g[j] for j in g} J = {} for j in range(-4, 1): s = mpf(0) for m in range(0, j + 5): s += c ** m / factorial(m) * K[j - m] J[j] = s / x return J def amflow_heldout_gate(data, orders, label=""): """Stage 5: the served tower at the two points held out of the fit vs the raw goal-40 / goal-70 AMFlow strings of record, per order, in the J normalisation; PASS iff every (point, order) agrees with BOTH strings to >= min(pair-certified digits, working dps - 5). Appends to FAILS by name; returns (n_cells_ok, n_cells, min_agreement, tower_walls).""" ref = data["references"]["amflow_heldout"] norm = ref["normalisation"] bar_dps = mp.dps - 5 print(f"\n== 5. AMFLOW HELD-OUT GATE{label}: the served tower at the two " f"points held out of the fit vs the raw goal-40 / goal-70 AMFlow " f"strings (working dps {mp.dps}; bar = min(pair-certified digits, " f"{bar_dps})) ==") print(f" data file sha256 {DATA_SHA256[:16]}... = pin (verified before " f"any computation); g -> J: K = {norm['K_from_g_factor']} g, " f"J_j = (1/x) [e^({norm['exponent_gamma_coefficient']} gamma_E eps) K]_j " f"at s = -1, t = -x") print(f" {'point':>5} {'order':>6} {'tower J':<34} {'goal-40 d':>9} " f"{'goal-70 d':>9} {'pair d':>7} {'(record)':>8} {'bar':>4} verdict" f" [pair d = the goal-40 vs goal-70 agreement recomputed from the " f"shipped strings; (record) = the vendored count, the bar's source]") n_ok = n_all = 0 dmin = float("inf") walls = [] for xs, pt in ref["points"].items(): fx = Fraction(xs); x = mpf(fx.numerator) / fx.denominator t0 = time.time() tw, tv, certs = tower_orders(orders, x) # the same code path as --point walls.append(time.time() - t0) J = j_from_g(tw, x, norm) for j in range(-4, 1): m40 = ball_mid(pt["g40"]["orders"][str(j)]["re"]) m70 = ball_mid(pt["g70"]["orders"][str(j)]["re"]) d40 = rel_digits(J[j], m40) d70 = rel_digits(J[j], m70) with mp.workdps(max(sig_digits(m40), sig_digits(m70)) + 10): pair = min(rel_digits(mpf(m40), m70), float(sig_digits(m40))) rec = pt["pair_certified_digits_per_order"][str(j)] rec_ok = abs(round(pair, 2) - float(rec)) <= 0.011 bar = min(int(rec), bar_dps) # the record's count, not the live one ok = d40 >= bar and d70 >= bar and rec_ok n_all += 1 n_ok += ok dmin = min(dmin, d40, d70) print(f" {xs:>5} eps^{j:+d} {mp.nstr(J[j], 30):<34} {d40:9.2f} " f"{d70:9.2f} {pair:7.2f} {('(' + str(rec) + ')'):>8} {bar:4d} " f"{'PASS' if ok else 'FAIL'}") short = [f"goal-{g} {d:.2f}d" for g, d in (("40", d40), ("70", d70)) if d < bar] if short: FAILS.append(f"AMFLOW HELD-OUT GATE{label} x={xs} eps^{j:+d}: " f"{' and '.join(short)} < bar {bar}d") if not rec_ok: FAILS.append(f"AMFLOW HELD-OUT GATE{label} x={xs} eps^{j:+d}: pair " f"digits from the strings {pair:.2f} != the data " f"file's figure {rec}") print(f" x = {xs}: certified transport error |Delta g_j| <= " f"{mp.nstr(max(certs.values()), 3)} for every eps order; tower " f"{walls[-1]:.2f}s") verdict = "PASS" if n_ok == n_all else "FAIL" print(f" AMFLOW HELD-OUT GATE{label}: {verdict} ({n_ok}/{n_all} (point, " f"order) cells at or above the bar against both strings; min " f"agreement {dmin:.2f} digits; the two points' tower " f"{sum(walls):.2f}s)") return n_ok, n_all, dmin, walls # ---------------------------------------------------------------------- main def digits(a, b): d = abs(a - b) if d == 0: return mp.dps r = d / max(abs(b), mpf(1)) return float(-log(r) / log(10)) def run(points, orders, do_control=True, label=""): results = {} for xs in points: fx = Fraction(xs); x = mpf(fx.numerator) / fx.denominator t0 = time.time() tw, tv, certs = tower_orders(orders, x) t_tower = time.time() - t0 t0 = time.time() orc = oracle_orders(x) t_orc = time.time() - t0 LS = 1 / (2 * x) # LS = 1/(-2 s^2 t) at (s,t)=(-1,-x) print(f"\n== point x = {xs}{label} " f"(tower {t_tower:.2f}s, oracle {t_orc:.2f}s) ==") print(f" LS = 1/(-2 s^2 t) = {mp.nstr(LS, 20)} at (s,t)=(-1,{mp.nstr(-x,8)})") dmin = mp.dps for j in range(-4, 1): dj = digits(tw[j], orc[j]) dmin = min(dmin, dj) print(f" eps^{j:+d} g = {mp.nstr(tw[j], 30)} " f"agreement vs oracle: {dj:.1f} digits") print(f" certified transport error: |Delta g_j| <= " f"{mp.nstr(max(certs.values()), 3)} for every eps order " f"(accumulated per-step geometric tail bound, r=1/4, guard 8; " f"BOUND, not estimate)") print(f" min agreement at this point: {dmin:.1f} digits " f"(ceiling = working dps {mp.dps}; need >= {mp.dps - 5})") if dmin < mp.dps - 5: FAILS.append(f"x={xs}: min agreement {dmin:.1f}d < {mp.dps-5}d") results[xs] = (tw, orc, dmin) if do_control: t0 = time.time() with mp.workdps(min(mp.dps, 60)): m1, z0 = mb_k1t(x) k0 = k0t_logs(x) c_m1 = -4 * tw[-1] - k0[-1] c_0 = -4 * tw[0] - k0[0] d1, d0 = digits(m1, c_m1), digits(z0, c_0) print(f" CONTROL Mellin-Barnes (independent numerical " f"integration, {time.time()-t0:.1f}s, dps<=60):") print(f" K1t eps^-1: {d1:.1f} digits K1t eps^0: {d0:.1f} digits" f" (need >= 40)") if d1 < 40.0 or d0 < 40.0: FAILS.append(f"x={xs}: MB control {d1:.1f}/{d0:.1f}d < 40d") # positive control pc = digits(tv[(0, 1)], -polylog(2, -x)) print(f" positive control G(0,-1;x) vs -Li2(-x): {pc:.1f} digits") return results def main(): ap = argparse.ArgumentParser(description=__doc__.splitlines()[0], epilog=__doc__[__doc__.index("Usage:"):], formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--dps", type=int, default=50) ap.add_argument("--point", action="append", default=[], help="extra x=t/s>0 (fraction or decimal), repeatable") ap.add_argument("--check", action="store_true", help="rerun tower+oracle at dps+60 and report drift") ap.add_argument("--no-control", action="store_true", help="skip the Mellin-Barnes control integration") ap.add_argument("--mutate", action="store_true", help="perturb one vendored tower coefficient; run must exit nonzero") args = ap.parse_args() T0 = time.time() mp.dps = args.dps + 10 data, orders = load_data() if args.mutate: c, zk, pw, word = orders[0][0] orders[0][0] = (c + 1, zk, pw, word) print("[mutate] one vendored eps^0 tower coefficient +1 " "(this run MUST exit nonzero)") print(f"massless double box; {data['name']}") print(f"data: {DATA_PATH}") print(f"normalization: {data['normalization']['g']}") print(f"domain: Euclidean planar region, x = t/s > 0 " f"(oracle soft-capped ~100d for x in (0.9,1.1))") print(f"working dps = {mp.dps} (requested {args.dps} + 10 guard)") points = list(GATE_POINTS) + args.point res = run(points, orders, do_control=not args.no_control) amflow_heldout_gate(data, orders) if args.check: print(f"\n== --check: rerun at dps+60 = {args.dps + 70} ==") with mp.workdps(args.dps + 70): drift = mpf(0) for xs in points: fx = Fraction(xs); x = mpf(fx.numerator) / fx.denominator tw2, _, _ = tower_orders(orders, x) orc2 = oracle_orders(x) tw1 = res[xs][0] for j in range(-4, 1): drift = max(drift, abs(tw1[j] - tw2[j]) / max(abs(tw2[j]), mpf(1))) dj = digits(tw2[j], orc2[j]) print(f" x={xs}: deep-check min digits " f"{min(digits(tw2[j], orc2[j]) for j in range(-4,1)):.1f}") print(f" max tower drift vs dps+60 rerun: {mp.nstr(drift, 3)} " f"(need <= 1e-{args.dps - 5})") if drift > mpf(10) ** (-(args.dps - 5)): FAILS.append(f"--check drift {mp.nstr(drift, 3)} > 1e-{args.dps-5}") amflow_heldout_gate(data, orders, label=" (--check rerun at dps+60)") print(f"\ntotal wall time: {time.time() - T0:.2f}s") if FAILS: print(f"OVERALL FAIL: {'; '.join(FAILS)}") return 1 print("OVERALL PASS") return 0 if __name__ == "__main__": sys.exit(main())