#!/usr/bin/env python3 r"""eval_row33.py -- LBL3VP (QED parent: light-by-light box x dotted VP-type self-energy insertion, nu=[1,1,1,2,1,2,1,1], 3 loops), row 33: standalone arbitrary-precision evaluator, ZERO AMFlow at eval time. Imports: mpmath + sympy (exact-rational DE parsing) + lbl3disp/k33lib/rho33lib (same dir). No AMFlow / Kira / IBP / network / outside-tree reads at runtime; the only file inputs are the exact-rational class-(c) w-DEs in row33_data.json and the independent reference literals (never consumed by the computation). SCOPE (what is computed LIVE here): [A] fixed-eps dispersive value at the kinematic point --kin s,t[,m2] (default the REFERENCE point (s,t,m2)=(-1,-1/3,1); build B1, 2026-09-04): I_disp(eps) = (1/pi) int_1^inf rho_VP(w';eps) K_P(w';eps) dw' - kernel K_P(w';eps): fixed-eps Taylor-step transport (RatDE) of the exact -rational 8-master one-loop box w'-DE, seeded at w'=5 by the EXACT fixed-eps Feynman-parameter masters (k33lib.masters_at: closed tadpoles, 1-fold bubbles/triangles, 2-fold box; all Euclidean, geometric GL); the threshold Taylor block K_n(eps) is the EXACT moment 2-fold at w'=1 (k33lib.kn_moments). NO AMFlow seed, no fit. - spectral density rho_VP(w';eps): fully SYMBOLIC (rho33lib.Rho33) -- closed 2-body form on (1,9); Frobenius series rho_exact - c(eps) psi_0 on [9, wsw], with c(eps) the ANALYTIC 3-massive-cut threshold coefficient (Gamma functions; reproduces the dev-oracle-measured c at 2^-6,-7,-8 to all digits) and psi the exact-rational Frobenius solution of the {0,1} 2x2 block at w=9, exponent d-2; real-axis Im-4-vector transport of the EXACT sigvp w-DE with the symbolic seed (ImM0, ImM1, 0, ImTB) above w' = wsw (no complex deflection: the Im parts satisfy the same real rational DE on (9,inf)). - tanh-sinh panel quadrature + one Richardson step (live); compared with the independent 3-loop AMFlow reference grid (eps=2^-k, k=4..18). [B] Laurent pole layers at ARBITRARY (s,t,m2) in the documented domain: eps^-2: I_VP^(-2) = -K_2^(0)(s,t,m2), K_2^(0) extracted live by a Chebyshev-Vandermonde threshold-Taylor fit of the CLOSED dilogarithmic one-loop box kernel (the closed form in lbl3disp, Vieta-stabilized; no stored value, any Euclidean point); eps^-1: I_VP^(-1) = -K_2^(1) + B[sigma_-1]^(0). B[sigma_-1]^(0) (three tanh-sinh quadratures of the closed-form kernel against the exact 2-body eps^-1 density layer) is live at any Euclidean point; K_2^(1) (order-eps threshold Taylor coefficient) is derived at ANY point by the exact parametric threshold 2-fold (K2_threshold_eps01, the exact threshold two-fold, (s,t)-generalized; since 2026-07-05) -- the eps^-1 layer is complete everywhere in the documented domain. Depth knobs (fit degree nfit, quad/kernel dps, kn-moment count, far-tail cutoff) are wired to --dps: frozen at the reference values for dps <= 110, scaled by MEASURED rates above it, so delivered digits track the requested precision. [C] eps^0 Laurent layer via --eps0 (since 2026-07-05): pole-subtracted multi-eps Neville extraction (rational-eps grid below the oracle grid, layers at lifted precision), compared with the independent 50-d eps0 reference string with a leave-one-out error bar. The fixed-eps comparison [A] plus the two-route fit-vs-fold layer check are the certification. AGREEMENT TABLE (measured 2026-07-03, --dps 110, LMAX 6/7, N_TAY 104, SF 0.28; oracle = INDEPENDENT 3-loop AMFlow reference grid, never consumed by the computation; archived AMFlow-seeded run in parens): I_disp fixed-eps, ALL EIGHT 2^-6..2^-13 vs 3-loop AMFlow grid: 2^-6 37.84 | 2^-7 37.59 | 2^-8 37.61 | 2^-9 37.77 | 2^-10 38.01 | 2^-11 38.27 | 2^-12 38.56 | 2^-13 38.85 d (MIN 37.59 d; archived seeded 31.70-32.79 d -- the eps^20-truncation seed cap is gone, +~6 d everywhere) per-point: 353-389 s at dps 110 (seed ~90 s, K transport ~200 s, rho ~90 s); self-checks: box-seed 2-prec 122 d, K_n 2-prec 119 d, backward rho transport 68.8-73.1 d, quad Richardson 38-39 d (= the agreement limiter, work_dps=DPS-25). eps^-2 layer vs independent 50-d reference literal ~45 d eps^-1 layer vs independent 50-d reference literal ~39.6 d (K_2^(1) de-AMFlowed, 60.4 d vs the retired AMFlow constant) CROSS-CHECK: the independent evaluator lbl3vp-evaluate.py (same page) computes the SAME 3-loop object by a different route and agrees to 37.7-40.08 d with the same reference values -- consistent with this script's 37.59-38.85 d (no flag). INPUTS (row33_data.json; class-(c) exact-rational data + independent reference literals (comparison only) -- NO AMFlow seed literals enter the computation): (1) sigvp_de, box1eq_de: exact-rational 4- and 8-master w-DE matrices (Kira IBP connection, validated max-diff 0). The A-linear-in-d entries are substituted at fixed eps at runtime. (2) independent reference values (comparison only): 3-loop AMFlow fixed-eps grid + eps^-2/eps^-1 Laurent layers (50-320 d literals). The former AMFlow boundary seeds (sigvp_bnd_w5, box1eq_bnd_w5, Kn) remain in the JSON but are NO LONGER READ: the kernel seed, threshold moments and rho density are now computed symbolically (k33lib / rho33lib). STATED GAPS: * fixed-eps I_disp at a non-reference point needs that point's one-loop box w'-DE (--box-de, reconstructed from IBP samples at that point) -- the k33lib seeds are (s,u)-parameterized since build B1 but the DE is a file input; --eps0 stays reference-only (its comparison literals are the reference's). The pole layers [B] are complete (eps^-2 AND eps^-1) at any point in the domain. INTERFACE: python3 eval_row33.py [--dps N] comparison demo: I_disp(2^-6) + pole layers, all compared live python3 eval_row33.py --point 2^-9 I_disp(eps) at another eps (compared if on k=4..18 grid) python3 eval_row33.py --point=-2,-1/5,1 pole layers at another Euclidean (s,t,m2), live python3 eval_row33.py --check [...] rerun at dps+60, print two-precision agreement python3 eval_row33.py --eps0 [--dps N] eps^0 Laurent layer, live (grid override: --eps0-ks) build B1 (2026-09-04) flags: --kin s,t[,m2] kinematic point for the fixed-eps value (rationals; same domain as the 3-value --point); default reference --box-de FILE the point's one-loop box w'-DE json (schema of DATA['box1eq_de']: masters, var, dim_var, A[8][8]); required off the reference; default DATA['box1eq_de'] --oracle-json FILE AMFlow black_box output at the point (result[0] = nu [1,1,1,2,1,2,1,1] samples): compare the value live --no-oracle never compare (write the value only) --out FILE write {k, eps, I_disp (full-dps string), diag, ...} --dry-run parse, load the DE, derive its singularities, write --out WITHOUT the key I_disp (status DRY_RUN); no eval (the served front end lbl3vp-points-evaluate.py drives these flags) DOMAIN: * eps: rational, 0 < eps <= 1/16 (below-threshold grid side). * (s,t,m2): s<0, t<0, m2>0, u=-s-t<4*m2 (deep-Euclidean; closed-dilog kernel domain). m2 is handled by exact rescaling: the script evaluates Ihat(s/m2, t/m2) := I_VP(s/m2, t/m2, 1); the physical normalization is I_VP(s,t,m2) = m2^(3d/2-10) * Ihat (3 loops, sum(nu)=10, d=4-2eps), which mixes Laurent layers by -(4+3eps)*log(m2) -- printed as a reminder when m2 != 1. * --dps sets the transport working precision; Taylor order, quad levels, work-dps and the threshold patch width scale with it (VP_TAYLOR_N / VP_LMAX / VP_STEP_FRAC / VP_VS env knobs honored). The reference comparison config is --dps 110 (default); the quad Richardson (work_dps=DPS-25) is the limiter. CHANGELOG: 2026-09-04 build B1 (the (s,t)-general build of this evaluator, the fixed-eps evaluator that checks the further (s,t) points): --kin s,t[,m2] sets the module point KIN (kin_hat/is_reference); k33lib seeds receive (shat, uhat); --box-de loads the point's one-loop box w'-DE (else the reference DE in the data file); the transport singularity list is DERIVED from the loaded DE (de_singularities: roots of every denominator at fixed eps, plus 0 and 1) instead of the reference literals; oracle comparison only at the reference (DATA) or from --oracle-json; --out writes the full-dps value + diagnostics; --dry-run. Regression check: the default invocation and --point 2^-7 at dps 110 unchanged (the wall-time tokens aside). 2026-07-05 depth-knob wiring: depth knobs wired to --dps, FROZEN at the reference values for dps <= 110 (regression-checked: default output numbers byte-identical) and scaled by measured rates above -- pole-layer fit degree nfit(DPS) (1.398 d/term window rate), kn-moment count, WP_MAX far-tail cutoff, Vandermonde solve guard. K_2^(1) central difference -> exact parametric threshold 2-fold at ANY (s,t) (K2_threshold_eps01; dlt=1e-30 knob deleted; off-reference eps^-1 gap CLOSED; fit_chk re-referenced to the fold = live two-route depth check). NEW --eps0 mode: Route N pole-subtracted Neville extraction of the eps^0 Laurent layer. 2026-07-03 de-AMFlowed row 33. Box-kernel seed at w'=5 -> exact fixed-eps Feynman-parameter masters (k33lib.masters_at); threshold K_n(eps) -> exact moment 2-folds (k33lib.kn_moments); rho_VP -> fully symbolic (rho33lib: closed 2-body + analytic c(eps) Frobenius series + exact-DE real-axis Im-transport); pole-layer K_2^(1) -> central diff of the exact moment. Retired the sigvp_bnd_w5 / box1eq_bnd_w5 / Kn AMFlow seeds (kept in JSON, unread) and the PF-sweep label; only the exact-rational w-DEs remain as file inputs. Agreement 32.09 -> 37.59-38.85 d (all 8 eps, dps 110). Default --dps 70 -> 110. 2026-07-03 (interim) built from the round-1 blog machinery + lbl3disp.py; pole layers generalized to arbitrary (s,t,m2); CLI to rowscripts convention. """ import os, json, time, bisect, sys, math import mpmath as mp import sympy as sp HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, HERE) import lbl3disp as L import k33lib from rho33lib import Rho33 DATA = json.load(open(os.path.join(HERE, 'row33_data.json'))) ORC = DATA['oracles'] # ---- kinematic point of the fixed-eps evaluation (build B1, 2026-09-04): set # by --kin s,t[,m2] as exact sympy Rationals; default = the reference. # The evaluator works on Ihat(s/m2, t/m2) (see DOMAIN). mpf values are # built at CALL time, never at import time (an import-time mpf freezes that dps). ---- KIN = (sp.Integer(-1), sp.Rational(-1, 3), sp.Integer(1)) BOX_DE = None # the one-loop box w'-DE in use; None -> the reference DE in the data file BOX_DE_PATH = None NO_ORACLE = False # --no-oracle: never compare ORACLE_GRID = None # --oracle-json: {k: value string} at the point ORACLE_PATH = None def set_kin(s_, t_, m2_=1): global KIN KIN = (sp.Rational(s_), sp.Rational(t_), sp.Rational(m2_)) def kin_hat(): """(shat, that, uhat) = (s/m2, t/m2, -(s+t)/m2) as exact sympy Rationals.""" s_, t_, m2_ = KIN return s_ / m2_, t_ / m2_, -(s_ + t_) / m2_ def is_reference(): sh, th, _ = kin_hat() return sh == -1 and th == sp.Rational(-1, 3) def box_de(): return BOX_DE if BOX_DE is not None else DATA['box1eq_de'] def kin(): sh, th, _ = kin_hat() return mp.mpf(sh.p) / sh.q, mp.mpf(th.p) / th.q, mp.mpf(1) def agree_d(a, b): if a == b: return mp.inf return float(-mp.log10(abs(a - b) / abs(b))) # ---- depth-knob scaling (2026-07-05, measured rates) ------------- # At --dps <= LEGACY_DPS every depth knob freezes to the reference-run # value (the default invocation's output is unchanged -- the regression bar); # above it each knob scales with its MEASURED rate + safety margin so the # delivered digits track --dps instead of a frozen method depth. LEGACY_DPS = 110 def nfit_for(DPS): """Chebyshev-fit degree for the pole-layer kernel Taylor fit. Measured: c2 accuracy = 1.61 d/order at r=0.05, flat in dps (nfit 28/44/60 -> 45.03/70.83/96.58 d); the |u|<0.04 Horner window carries log10(1/0.04) = 1.398 d/term -- the binding (slower) rate. Target: DPS + 8 digits at the window rate.""" if DPS <= LEGACY_DPS: return 28 return max(28, math.ceil((DPS + 8) / 1.398) + 1) # ===================================================================== # Fixed-eps rational DE transport (Taylor-step ODE integrator), ported from # the fixed-eps dispersive prototype of this work. # ===================================================================== class RatDE: def __init__(self, de_json, eps_rat, sings): self.masters = [tuple(m) for m in de_json['masters']] self.n = len(self.masters) W, D = sp.symbols('w d') d_sym = sp.Rational(4) - 2 * eps_rat self.funcs = [[None] * self.n for _ in range(self.n)] for i in range(self.n): for j in range(self.n): e = sp.sympify(de_json['A'][i][j], locals={'w': W, 'd': D}) if e == 0: continue e = sp.cancel(e.subs(D, d_sym)) num, den = sp.fraction(e) nc = [mp.mpf(sp.Rational(c).p) / mp.mpf(sp.Rational(c).q) for c in sp.Poly(num, W).all_coeffs()] dc = [mp.mpf(sp.Rational(c).p) / mp.mpf(sp.Rational(c).q) for c in sp.Poly(den, W).all_coeffs()] self.funcs[i][j] = (nc, dc) self.sings = [mp.mpc(s) for s in sings] def nearest(self, w): w = mp.mpc(w) return min(abs(w - s) for s in self.sings) def _shift(self, coefs, w0): p = [mp.mpc(coefs[0])] for c in coefs[1:]: pn = [mp.mpc(0)] * (len(p) + 1) for m, a in enumerate(p): pn[m] += a * w0 pn[m + 1] += a pn[0] += c p = pn return p def step(self, M, w0, h, N): n = self.n ser = [[] for _ in range(N + 1)] for i in range(n): for j in range(n): cd = self.funcs[i][j] if cd is None: continue nc, dc = cd ph = self._shift(nc, w0) qh = self._shift(dc, w0) q0 = qh[0] r = [mp.mpc(1) / q0] for m in range(1, N + 1): s = mp.mpc(0) for l in range(1, min(m, len(qh) - 1) + 1): s += qh[l] * r[m - l] r.append(-s / q0) for m in range(N + 1): s = mp.mpc(0) for l in range(min(m, len(ph) - 1) + 1): s += ph[l] * r[m - l] if s != 0: ser[m].append((i, j, s)) C = [list(M)] for m in range(N): s = [mp.mpc(0)] * n for l in range(m + 1): cv = C[m - l] for (ri, ci, a) in ser[l]: s[ri] += a * cv[ci] inv = mp.mpf(1) / (m + 1) C.append([x * inv for x in s]) v = [mp.mpc(0)] * n hp = mp.mpc(1) for m in range(N + 1): cm = C[m] for i in range(n): v[i] += cm[i] * hp hp *= h return v, C def path(self, M0, w_from, targets, N, sf, deflect=None): n = self.n M = list(M0) w = mp.mpc(w_from) out = {} ti = 0 nT = len(targets) end = mp.mpc(targets[-1]) sgn = 1 if end.real > w.real else -1 wps = [] for s in sorted(deflect or [], key=lambda x: sgn * x): if min(w.real, end.real) < s < max(w.real, end.real): r = mp.mpf('0.35') wps += [mp.mpc(s - sgn * r, 0), mp.mpc(s, r), mp.mpc(s + sgn * r, 0)] wps.append(end) for wp in wps: while abs(w - wp) > mp.mpf('1e-200'): d = self.nearest(w) hmax = sf * d dirn = wp - w stp = dirn if abs(dirn) <= hmax else dirn / abs(dirn) * hmax Mnew, C = self.step(M, w, stp, N) w_next = w + stp if abs(w.imag) < mp.mpf('1e-50') and abs(w_next.imag) < mp.mpf('1e-50') and stp != 0: while ti < nT: t = mp.mpc(targets[ti]) frac = ((t - w) / stp).real if frac < -mp.mpf('1e-50') or frac > 1 + mp.mpf('1e-50'): break hh = t - w v = [mp.mpc(0)] * n hp = mp.mpc(1) for m in range(len(C)): cm = C[m] for i in range(n): v[i] += cm[i] * hp hp *= hh out[targets[ti]] = v ti += 1 M = Mnew w = w_next w = wp while ti < nT: out[targets[ti]] = list(M) ti += 1 return out, M def to(self, M0, w_from, w_to, N, sf, deflect=None): _, M = self.path(M0, w_from, [w_to], N, sf, deflect) return M def de_singularities(de_json, eps_rat, dps, extra=(0, 1)): """Singular points of the fixed-eps rational DE for the Taylor-step size rule (RatDE.nearest): the distinct roots of every denominator of A after the substitution d = 4-2eps and cancellation (exactly RatDE's parsing), factored over Q, each irreducible factor's roots by mp.polyroots at dps+20 (linear factors exactly), parts below 10^-(dps+5) zeroed, duplicates merged, plus `extra` (0 and 1: the tadpole/threshold points, kept as in the record's list). Returned rounded to the CALLER's precision. At the reference DE (denominators w, w-1, w^2+1, w^2+4, 3w^2-2w+3) this reproduces the record's literal list {0, 1, +-i, +-2i, 1/3 +- i sqrt(8)/3}.""" W, D = sp.symbols('w d') d_sym = sp.Rational(4) - 2 * sp.Rational(eps_rat) factors = {} for row in de_json['A']: for ent in row: e = sp.sympify(ent, locals={'w': W, 'd': D}) if e == 0: continue _, den = sp.fraction(sp.cancel(e.subs(D, d_sym))) for f, _m in sp.factor_list(sp.Poly(den, W))[1]: p = sp.Poly(f, W) if p.degree() < 1: continue factors[tuple(sp.Rational(c) for c in p.all_coeffs())] = p old = mp.mp.dps mp.mp.dps = dps + 20 try: tol = mp.mpf(10) ** (-(dps + 5)) roots = [mp.mpc(x) for x in extra] for key in sorted(factors, key=lambda k: (len(k), [str(c) for c in k])): if len(key) == 2: r = -key[1] / key[0] rts = [mp.mpc(mp.mpf(r.p) / r.q)] else: rts = mp.polyroots([mp.mpf(c.p) / c.q for c in key], maxsteps=500, extraprec=mp.mp.prec) for r in rts: r = mp.mpc(r) r = mp.mpc(r.real if abs(r.real) > tol else 0, r.imag if abs(r.imag) > tol else 0) if all(abs(r - q) > tol for q in roots): roots.append(r) finally: mp.mp.dps = old return [+r for r in roots] def vs_for(DPS): """Half-width of the threshold panel patch served by the reference 20-term K_n(eps) Taylor series (VP_VS overrides). Truncation of that series, ~|c_20| vs^18 * G/(9 |I|), is an analytic-input digit cap; it is lifted by shrinking vs (DE transport covers the rest of the panel).""" env = os.environ.get('VP_VS') if env: return mp.mpf(env) return mp.mpf('0.02') if DPS <= 90 else mp.mpf('0.004') if DPS <= 200 else mp.mpf('0.0015') def load_laur(block): """Stored boundary Laurent data {idx: {order: [re,im]}} -> mpc dict.""" out = {} for idx, orders in block.items(): key = tuple(int(x) for x in idx.split(',')) out[key] = {int(o): mp.mpc(mp.mpf(v[0]), mp.mpf(v[1])) for o, v in orders.items()} return out def at_eps(l, eps, kmin, kmax): return sum(l.get(k, mp.mpc(0)) * eps ** k for k in range(kmin, kmax + 1)) def ts_nodes(a, b, level, dps): old = mp.mp.dps mp.mp.dps = dps + 20 ts = mp.calculus.quadrature.TanhSinh(mp.mp) prec = mp.mp.prec pts = set() for deg in range(1, level + 1): for x, w in ts.get_nodes(mp.mpf(a), mp.mpf(b), deg, prec): pts.add(+mp.mpf(x)) mp.mp.dps = old return sorted(pts) def quad_full(f, a, b, level, dps): old = mp.mp.dps mp.mp.dps = dps + 20 ts = mp.calculus.quadrature.TanhSinh(mp.mp) prec = mp.mp.prec results = [] for deg in range(1, level + 1): nodes = ts.get_nodes(mp.mpf(a), mp.mpf(b), deg, prec) results.append(ts.sum_next(f, nodes, deg, prec, results, False)) mp.mp.dps = old return results[-1] # ===================================================================== # Closed-form (dilog-level) one-loop box kernel K(w') at eps^0, at ARBITRARY # Euclidean (shat, that) = (s/m2, t/m2), m2 rescaled to 1: delegated to the # shared lib lbl3disp (Vieta-stabilized 4-root partial fractions, # accurate to arbitrary M5). # ===================================================================== def box_closed_dilog(shat, that, M5, dps=40): return L.box1(shat, that, mp.mpf(1), M5, dps) def kernel_taylor_live(shat, that, nfit=28, r='0.05', dps=100): """Taylor coefficients c_n of K(1+u) at u=0 from the CLOSED-FORM kernel: degree-nfit polynomial fit on Chebyshev nodes |u|<=r (Vandermonde solve). Nearest kernel singularity is at w'=0 (|u|=1), so the degree-(nfit) truncation limits c_2 to about (nfit-1)*log10(1/r) digits (measured 1.61 d/order). The solve guard grows with nfit (conditioning headroom) for nfit > 28; at the legacy depth it is the reference +40.""" sg = 40 + (0 if nfit <= 28 else nfit // 2) old = mp.mp.dps mp.mp.dps = dps + sg r = mp.mpf(r) us = [r * mp.cos(mp.pi * j / nfit) for j in range(nfit + 1)] vals = mp.matrix([[box_closed_dilog(shat, that, 1 + u, dps)] for u in us]) mp.mp.dps = dps + sg V = mp.matrix(nfit + 1, nfit + 1) for i, u in enumerate(us): p = mp.mpf(1) for j in range(nfit + 1): V[i, j] = p p *= u c = mp.lu_solve(V, vals) mp.mp.dps = old return [c[j] for j in range(nfit + 1)] # ===================================================================== # Per-eps dispersive assembly (ported from disp_vp_fixedeps.py, this work) # ===================================================================== def I_disp_at_eps(eps_rat, sig_de, box_de, DPS, N_TAY, SF, LMAX, work_dps=None, verbose=True, s_rat=None, u_rat=None): """I_disp(eps) at (shat, that) = (s/m2, t/m2) -- the module point KIN unless (s_rat, u_rat) are given as exact rationals (u = -s-t) -- for any rational 0 < eps <= 1/16; box_de = the one-loop box w'-DE AT THAT POINT. ZERO AMFlow at eval time: the box-kernel DE seed at w'=5 and the threshold Taylor block K_n(eps) are exact fixed-eps Feynman-parameter integrals (k33lib, (s,u)-parameterized); rho_VP is symbolic (rho33lib, (s,t)-free). Only the exact-rational sigvp/box w-DEs (class-c IBP data) and the independent reference values are file inputs. Everything below RUNS now. work_dps (quadrature assembly precision) scales with DPS unless pinned.""" t_e = time.time() if s_rat is None or u_rat is None: sh, th, uh = kin_hat() s_rat = sh if s_rat is None else sp.Rational(s_rat) u_rat = uh if u_rat is None else sp.Rational(u_rat) if work_dps is None: work_dps = max(45, DPS - 25) eps_rat = sp.Rational(eps_rat) mp.mp.dps = DPS eps = mp.mpf(eps_rat.p) / eps_rat.q G = mp.gamma(1 + eps) * mp.gamma(1 - eps) / (eps * mp.gamma(2 - 2 * eps)) # --- de-AMFlowed kernel seeds: exact fixed-eps Feynman-parameter masters at # w'=5 (k33lib.masters_at) + exact threshold Taylor moments K_n(eps) # (k33lib.kn_moments); no fit, no AMFlow. --- seed_work = DPS + 30 n_box = max(120, int(1.65 * (DPS + 25))) # 0.63 d/node measured n_kn = max(100, int(1.10 * (DPS + 25))) # 1.0 d/node measured # threshold Taylor moment COUNT: frozen legacy 20 at dps <= LEGACY_DPS # (byte-identical default); above, scaled with the |u|= vs] + ts_nodes(mp.mpf(1) + vs, bb, LN, work_dps)), reverse=True) K_tg_mid = sorted([v for v in n_p1b if v > 5]) K_tg_hi = sorted([v for v in (n_p2 + n_tail) if 9 < v < WP_MAX]) s_tg_hi = list(K_tg_hi) mp.mp.dps = DPS # === box-kernel transport K(w';eps): 8-master exact DE, seed = k33lib # Feynman-parameter masters at w'=5 (Mb0 above, in the DE master order); # the step-size singularity list is derived from the loaded DE (build # B1; the reference literal list is its value at the reference DE) === b_sings = de_singularities(box_de, eps_rat, DPS) de_b = RatDE(box_de, eps_rat, b_sings) assert [tuple(m) for m in de_b.masters] == [ (1, 1, 1, 1), (1, 0, 0, 0), (0, 1, 0, 0), (1, 0, 1, 0), (0, 1, 0, 1), (1, 1, 1, 0), (1, 1, 0, 1), (0, 1, 1, 1)], \ "box1eq DE master order != k33lib.masters_at order" K_at = {} t0 = time.time() for tg in (K_tg_lo, K_tg_mid): if tg: res, _ = de_b.path(Mb0, mp.mpf(5), tg, N_TAY, SF) for v in tg: K_at[v] = mp.re(res[v][0]) Mb95 = de_b.to(Mb0, mp.mpf(5), mp.mpf('9.5'), N_TAY, SF) for tg in [sorted([v for v in K_tg_hi if v <= mp.mpf('9.5')], reverse=True), sorted([v for v in K_tg_hi if v > mp.mpf('9.5')])]: if tg: res, _ = de_b.path(Mb95, mp.mpf('9.5'), tg, N_TAY, SF) for v in tg: K_at[v] = mp.re(res[v][0]) tK = time.time() - t0 K_keys = sorted(K_at.keys()) miss = [] def _near(keys, wp): i = bisect.bisect_left(keys, wp) for j in (i - 1, i, i + 1): if 0 <= j < len(keys) and abs(keys[j] - wp) < abs(wp) * mp.mpf('1e-40'): return keys[j] return None def K_P(wp): wp = mp.mpf(wp) u = wp - 1 if abs(u) < vs: acc = mp.mpf(0) for c in reversed(Kne[2:]): acc = acc * u + c return acc if wp >= WP_MAX: return mp.mpf(0) Kv = K_at.get(wp) if Kv is None: kk = _near(K_keys, wp) if kk is not None: Kv = K_at[kk] if Kv is None: miss.append(('K', wp)) if wp > mp.mpf('1e20'): return mp.mpf(0) raise KeyError(mp.nstr(wp, 20)) return (Kv - K0e - K1e * u) / (u * u) # === rho_VP(w';eps): fully SYMBOLIC (rho33lib), NO AMFlow, no complex # deflection. Closed 2-body form on (1,9); Frobenius series # rho_exact - c(eps) psi_0 on [9, wsw]; real-axis Im-4-vector transport # of the EXACT sigvp DE with the symbolic seed for w' > wsw. === R33 = Rho33(sys.modules[__name__], eps_rat, DPS) rho_at = {} t0 = time.time() ser_nodes = [v for v in s_tg_hi if v <= R33.wsw] tail_nodes = [v for v in s_tg_hi if v > R33.wsw] for v in ser_nodes: rho_at[v] = R33.rho_series(v) if tail_nodes: rho_at.update(R33.transport_tail(tail_nodes, N_TAY, SF)) tR = time.time() - t0 # live internal check: backward real-axis transport wsw -> 9.6 vs series if R33._de is None: R33._de = RatDE(sig_de, eps_rat, [0, 1, 9, -3]) back = R33._de.path(R33.seed_at(R33.wsw), R33.wsw, [mp.mpf('9.6')], N_TAY, SF)[0] xchk = agree_d(-mp.re(back[mp.mpf('9.6')][0]), R33.rho_series(mp.mpf('9.6'))) def rho_exact(wp, e): return R33.rho_exact(mp.mpf(wp)) diag = {'nK': len(K_at), 'tK': tK, 'nR': len(rho_at), 'tR': tR, 'xchk': xchk, 'seed_selfchk': seed_selfchk, 'kn_selfchk': kn_selfchk, 't_seed': t_seed, 's_rat': str(s_rat), 'u_rat': str(u_rat), 'n_box': n_box, 'n_kn': n_kn, 'nkn': nkn, 'seed_work': seed_work, 'work_dps': work_dps, 'vs': mp.nstr(vs, 6), 'WP_MAX': mp.nstr(WP_MAX, 6), 'b_sings_30d': [mp.nstr(x, 30) for x in b_sings]} if verbose: print(f" K transport {len(K_at)} pts {tK:.0f}s | rho {len(rho_at)} pts {tR:.0f}s" f" | backward-transport rho self-check {xchk:.1f} d" f" | seed 2-prec {seed_selfchk:.0f} d, Kn 2-prec {kn_selfchk:.0f} d ({t_seed:.0f}s)") rho_keys = sorted(rho_at.keys()) def rho(wp): wp = mp.mpf(wp) if wp <= 1: return mp.mpf(0) if wp < mp.mpf(9): return rho_exact(wp, eps) if wp >= WP_MAX: return mp.mpf(0) if wp in rho_at: return rho_at[wp] kk = _near(rho_keys, wp) if kk is not None: return rho_at[kk] miss.append(('rho', wp)) if wp > mp.mpf('1e20'): return mp.mpf(0) raise KeyError(mp.nstr(wp, 25)) # threshold panel [1,1.4]: peel the u^{-2eps} mode analytically (holds the poles) def f1(u): u = mp.mpf(u) return (2 + u) / (1 + u) ** (1 - eps) * K_P(1 + u) def f2(u): u = mp.mpf(u) return K_P(1 + u) / (1 + u) ** (1 - eps) f10 = 2 * K2e def assem(level): old = mp.mp.dps mp.mp.dps = work_dps + 20 A = f10 * Lb ** (-2 * eps) / (-2 * eps) A += quad_full(lambda u: (mp.mpf(u)) ** (-2 * eps) * (f1(u) - f10) / mp.mpf(u), mp.mpf(0), Lb, level, work_dps) A *= G * (1 - 2 * eps) B = -G * eps / 6 * quad_full(lambda u: mp.mpf(u) ** (1 - 2 * eps) * f2(u), mp.mpf(0), Lb, level, work_dps) P1 = A + B P2 = quad_full(lambda wp: rho(wp) * K_P(wp) / mp.pi, bb, mp.mpf(9), level, work_dps) P3 = quad_full(lambda wp: rho(wp) * K_P(wp) / mp.pi, mp.mpf(9), mp.mpf(200), level, work_dps) def ft(t): wp = mp.mpf(200) + mp.mpf(200) * (1 + t) / (1 - t) return rho(wp) * K_P(wp) / mp.pi * mp.mpf(200) * 2 / (1 - t) ** 2 P4 = quad_full(ft, mp.mpf(-1), mp.mpf(1), level, work_dps) mp.mp.dps = old return P1 + P2 + P3 + P4 I5 = assem(LMAX) I6 = assem(LMAX + 1) R56 = 2 * I6 - I5 diag['quad_agree'] = agree_d(I5, I6) diag['t_total'] = time.time() - t_e if verbose: print(f" quad levels L{LMAX}/L{LMAX+1} agree {diag['quad_agree']:.1f} d" f" -> Richardson ({diag['t_total']:.0f}s total)") return R56, diag # ===================================================================== # Exact parametric threshold 2-fold: K_2^{(0)}, K_2^{(1)} at ANY Euclidean # (shat, that) -- (s,t)-generalized (2026-07-05). # Closes the former "K_2^(1)(s,t) unknown off the reference point" gap and # retires the kn_moments central difference (and its dlt=1e-30 knob). # ===================================================================== def K2_threshold_eps01(shat, that, dps, guard=25): r"""eps^0 and eps^1 of the 2nd threshold Taylor coefficient K_2(eps) of the box kernel at w'=1, from the exact parametric 2-fold (no seeds, no fit). K analytic at w'=1 (A >= 1-u/4 > 0 on the simplex for u=-s-t < 4); differentiating K = Gamma(2+eps) II [A^{-1-eps}-C^{-1-eps}]/((1+eps)B) twice in w' (dA/dw' = dC/dw' = x1, C = A + BX) gives, at w'=1, K_2(eps) = [Gamma(2+eps)(2+eps)/2] II x1^2 [A^{-3-eps}-C^{-3-eps}]/B K_2^{(0)} = Q0, K_2^{(1)} = (3-2*euler)/2 * Q0 + Q1 with the cancellation-free divided-difference integrands (C-A = BX): Q0: x1^2 X (A^2+AC+C^2)/(A^3 C^3) Q1: x1^2 X [A^{-3} log1p(BX/A)/(BX) - (A^2+AC+C^2)/(A^3 C^3) ln C] Tensor Gauss-Legendre on (x1,x2)=(xi,eta(1-xi)); integrand analytic on the closed square -> geometric convergence (n = 2.2*dps + 40 nodes; rate calibrated at u=4/3, degrades smoothly as u -> 4 -- the fit_chk cross-check in pole_layers() catches any degradation live). Verified: two-precision 161.4 d at dps 160; vs stored AMFlow strings 121.0/120.6 d = full stored length.""" work = dps + guard n = int(mp.mpf('2.2') * dps) + 40 old = mp.mp.dps mp.mp.dps = work try: s = mp.mpf(shat) t = mp.mpf(that) u = -s - t nodes, weights = k33lib.gl_nodes(n, work) one = mp.mpf(1) half = one / 2 xs = [half * (one + tk) for tk in nodes] ws = [half * wk for wk in weights] Q0 = mp.mpf(0) Q1 = mp.mpf(0) for xi, wxi in zip(xs, ws): # xi = x1 om = one - xi x1sq_j = xi * xi * om * wxi # x1^2 * Jacobian * weight acc0 = mp.mpf(0) acc1 = mp.mpf(0) for et, wet in zip(xs, ws): x2 = et * om X = om - x2 A = one - u * x2 * X B = -s * xi + u * x2 BX = B * X C = A + BX A3 = A ** -3 R = (A * A + A * C + C * C) * A3 * C ** -3 ddlog = mp.log1p(BX / A) / BX if BX != 0 else 1 / A acc0 += wet * X * R acc1 += wet * X * (A3 * ddlog - R * mp.log(C)) Q0 += x1sq_j * acc0 Q1 += x1sq_j * acc1 K2_0 = Q0 K2_1 = (3 - 2 * mp.euler) / 2 * Q0 + Q1 mp.mp.dps = dps K2_0 = +K2_0 # round at dps BEFORE restoring caller precision K2_1 = +K2_1 # (module-level dps=15 would truncate otherwise) finally: mp.mp.dps = old return K2_0, K2_1 # ===================================================================== # Laurent pole layers, computed at runtime (route of this work's check_em1.py) # ===================================================================== def pole_layers(shat, that, at_ref=False, qdps=55, kdps=90, nfit=28): """eps^-2 and eps^-1 Laurent layers of I_VP at ARBITRARY Euclidean (shat, that) (m2 rescaled to 1), live from the closed-form eps^0 kernel. K_2^(0)/K_2^(1) are derived at ANY point by the exact parametric threshold 2-fold (K2_threshold_eps01) -- the eps^-1 layer is complete everywhere in the documented domain (the former reference-point-only central difference and its dlt knob are retired). fit_chk checks the Chebyshev-fit c2 against the fold: two INDEPENDENT routes, so it sees the nfit truncation that same-depth two-precision reruns cannot. Depth knobs (qdps, kdps, nfit) are wired to --dps by the caller and freeze to the reference values at the default; at_ref is retained for interface compatibility (it no longer changes the computation).""" mp.mp.dps = kdps + 40 c = kernel_taylor_live(shat, that, nfit=nfit, r='0.05', dps=kdps) K2_live = c[2] # de-AMFlowed K_2 layers at ANY Euclidean point: exact parametric 2-fold. K2_0_fold, K2_1 = K2_threshold_eps01(shat, that, kdps) # cross-check the live closed-dilog Taylor fit vs the independent fold fit_chk = agree_d(K2_live, K2_0_fold) K0, K1 = c[0], c[1] def K_P0(wp): wp = mp.mpf(wp) u = wp - 1 if abs(u) < mp.mpf('0.04'): acc = mp.mpf(0) for cc in reversed(c[2:]): acc = acc * u + cc return acc old = mp.mp.dps v = (box_closed_dilog(shat, that, wp, kdps - 20) - K0 - K1 * u) / (u * u) mp.mp.dps = old return +v def f1(wp): # [K_P(w')-K_2]/(w'-1), regular at threshold wp = mp.mpf(wp) u = wp - 1 if abs(u) < mp.mpf('0.04'): acc = mp.mpf(0) for cc in reversed(c[3:]): acc = acc * u + cc return acc return (K_P0(wp) - K2_live) / u def f2(wp): return K_P0(wp) / (mp.mpf(wp) - 1) def f3(wp): return K_P0(wp) / mp.mpf(wp) b = mp.mpf('1.4') mp.mp.dps = qdps I1 = mp.quad(f1, [1, b], method='tanh-sinh') I2 = mp.quad(f2, [b, 9], method='tanh-sinh') + mp.quad(f2, [9, 200], method='tanh-sinh') I2 += mp.quad(lambda t: f2(200 + 200 * (1 + t) / (1 - t)) * 200 * 2 / (1 - t) ** 2, [-1, 1], method='tanh-sinh') I3 = mp.quad(f3, [1, 9], method='tanh-sinh') + mp.quad(f3, [9, 200], method='tanh-sinh') I3 += mp.quad(lambda t: f3(200 + 200 * (1 + t) / (1 - t)) * 200 * 2 / (1 - t) ** 2, [-1, 1], method='tanh-sinh') Bsm1 = 2 * K2_live * (mp.euler + mp.log(b - 1)) + 2 * I1 + 2 * I2 - I3 # K_2^(1) de-AMFlowed above (exact parametric threshold 2-fold; any (s,t)) return {'I_em2': -K2_live, 'I_em1': -K2_1 + Bsm1, 'Bsm1': Bsm1, 'fit_chk': fit_chk, 'K2_0_fold': K2_0_fold, 'K2_1': K2_1} def parse_eps(txt): """eps from '2^-k', 'p/q' or an exact decimal string -> sympy Rational.""" txt = txt.strip() if txt.startswith('2^-'): return sp.Rational(1, 2 ** int(txt[3:])) return sp.Rational(txt) def eps_label(e): if e.p == 1 and (e.q & (e.q - 1)) == 0: return f"2^-{e.q.bit_length() - 1}" return f"{e.p}/{e.q}" def parse_ball(txt): """Arb ball string '[mid +/- rad]' -> mid as a decimal string, or None when the ball carries no digit ('[+/- r]'). A plain decimal passes through.""" t = txt.strip() if t.startswith('['): t = t.strip('[]') mid = t.split('+/-')[0].strip() return mid if mid else None return t def load_amflow_grid(path): """AMFlow black_box output at the point -> {k: value string} for the nu=[1,1,1,2,1,2,1,1] integral (result[0]; the row-33 object), samples at eps = 2^-k only; uncertified balls are dropped.""" G = json.load(open(path)) res = G['result'][0] idx = res['integral']['indices'] if idx[:8] != [1, 1, 1, 2, 1, 2, 1, 1]: raise SystemExit(f"--oracle-json: result[0] is not nu=[1,1,1,2,1,2,1,1]: {idx}") out = {} for smp in res['samples']: old = mp.mp.dps mp.mp.dps = 50 try: e = mp.mpf(smp['eps']['re'] if isinstance(smp['eps'], dict) else smp['eps']) k = int(mp.nint(-mp.log(e, 2))) ok = (mp.mpf(2) ** (-k) == e) finally: mp.mp.dps = old if not ok: continue v = smp['value']['re'] if isinstance(smp['value'], dict) else smp['value'] mid = parse_ball(v) if mid is not None: out[k] = mid return out def oracle_for(e): """Independent reference string for eps, if any: the --oracle-json grid at the point when given; else the DATA grid (k=4..18) AT THE REFERENCE POINT ONLY; nothing under --no-oracle or at a non-reference point.""" if e.p == 1 and (e.q & (e.q - 1)) == 0: k = e.q.bit_length() - 1 if ORACLE_GRID is not None: return ORACLE_GRID.get(k) if NO_ORACLE or not is_reference(): return None return ORC['I_AMF_grid'].get(str(k)) return None def run_gate(pts, DPS, N_TAY, SF, LMAX, verbose=True): """Evaluate I_disp at each rational eps in pts (sequential); where a independent reference value exists, recompute the agreement live as -log10 rel diff. Returns (min agreement over compared points -- mp.inf if none, {label: value}, {label: {'diag', 'oracle', 'agreement_d'}}).""" mind = mp.inf vals = {} info = {} for e in pts: val, diag = I_disp_at_eps(e, DATA['sigvp_de'], box_de(), DPS, N_TAY, SF, LMAX, verbose=False) mp.mp.dps = max(80, DPS + 20) vals[eps_label(e)] = val info[eps_label(e)] = {'diag': diag, 'oracle': None, 'agreement_d': None} if verbose: print(f" eps = {eps_label(e)}:") print(f" K transport {diag['nK']} pts {diag['tK']:.0f}s | rho {diag['nR']}" f" pts {diag['tR']:.0f}s | backward-transport rho self-check" f" {diag['xchk']:.1f} d | seed 2-prec {diag['seed_selfchk']:.0f} d," f" Kn 2-prec {diag['kn_selfchk']:.0f} d") print(f" quad levels L{LMAX}/L{LMAX+1} agree {diag['quad_agree']:.1f} d ->" f" Richardson ({diag['t_total']:.0f}s)") print(f" I_disp = {mp.nstr(val, 40)} (computed now)") ostr = oracle_for(e) if ostr is None: if verbose: if ORACLE_GRID is None and (NO_ORACLE or not is_reference()): print(" no oracle consumed at this point (--no-oracle / off the" " reference without --oracle-json); value computed, no comparison\n") else: print(" no independent reference value at this eps (grid: eps=2^-k, k=4..18);" " value computed, no comparison\n") continue orac = mp.mpf(ostr) d = agree_d(val, orac) mind = min(mind, d) info[eps_label(e)]['oracle'] = ostr info[eps_label(e)]['agreement_d'] = float(d) if verbose: if ORACLE_GRID is not None: print(f" I_AMF = {mp.nstr(orac, 40)} (independent 3-loop AMFlow reference grid at" f" the point, never consumed by the computation: {os.path.basename(ORACLE_PATH)})") print(f" agreement = {d:.2f} d (reference-point comparison at dps 110:" f" 37.59-38.85 d, for scale)\n") else: print(f" I_AMF = {mp.nstr(orac, 40)} (independent 3-loop reference value)") print(f" agreement = {d:.2f} d (symbolic reference target 37.59-38.85 d;" f" archived AMFlow-seeded 31.70-32.79 d)\n") return mind, vals, info def settings_for(DPS): """Precision-scaled defaults (env knobs override): Taylor order, tanh-sinh level, quadrature work-dps all grow with the requested transport dps so that raising --dps genuinely raises the reachable digits.""" N_TAY = int(os.environ.get('VP_TAYLOR_N', 0)) or max(80, int(0.95 * DPS)) LMAX = int(os.environ.get('VP_LMAX', 0)) or (5 if DPS <= 90 else 6 if DPS <= 200 else 7) return N_TAY, LMAX def one_pass(DPS, pts, kinpt, verbose=True): """One full evaluation pass at working precision DPS. pts: list of rational eps for fixed-eps I_disp at the reference point; kinpt: (s,t,m2) sympy Rationals for the live pole layers (or None). Returns {name: mp value} for the --check two-precision diff.""" mp.mp.dps = DPS out = {} if pts: SF = mp.mpf(os.environ.get('VP_STEP_FRAC', '0.28')) N_TAY, LMAX = settings_for(DPS) if verbose: if is_reference() and BOX_DE is None: print("fixed-eps dispersive value at the REFERENCE point s=-1, t=-1/3, m^2=1:") else: s_, t_, m2_ = KIN sh, th, uh = kin_hat() print(f"fixed-eps dispersive value at the point s={s_}, t={t_}, m^2={m2_}" f" [(shat,that)=({sh},{th}), uhat={uh}" f"{' = the reference point' if is_reference() else ''}]:") print(f" one-loop box w'-DE: {BOX_DE_PATH if BOX_DE_PATH else 'the reference DE in the data file'}") print(" I_disp(eps) = (1/pi) int_1^inf rho_VP(w';eps) K_P(w';eps) dw' (computed now)") print(" ZERO AMFlow: box-kernel DE seed at w'=5 = exact Feynman-parameter masters") print(" (k33lib); K_n(eps) = exact threshold moments; rho_VP symbolic (rho33lib:") print(" closed 2-body on (1,9), Frobenius rho_exact-c(eps)psi_0 on [9,wsw], exact-DE") print(" real-axis Im-transport above). Only the exact-rational w-DEs are file inputs.") print(f" settings: transport dps={DPS}, Taylor order={N_TAY}, quad levels" f" {LMAX}/{LMAX+1} (reference comparison run: dps 110, levels 6/7)") mind, vals, info = run_gate(pts, DPS, N_TAY, SF, LMAX, verbose) out['gate_min_d'] = mind out['_info'] = info for k, v in vals.items(): out[f'I_disp({k})'] = v if verbose and mp.isfinite(mind): print(f" min agreement over {{{', '.join(eps_label(e) for e in pts)}}}:" f" {float(mind):.2f} d [archived seeded run: >= 31.70 d]\n") if kinpt is not None: s_, t_, m2_ = kinpt sh_r, th_r = s_ / m2_, t_ / m2_ is_ref = (sh_r == -1 and th_r == sp.Rational(-1, 3)) qdps = max(55, DPS - 15) kdps = max(90, DPS + 20) mp.mp.dps = kdps + 60 shat = mp.mpf(sh_r.p) / sh_r.q that = mp.mpf(th_r.p) / th_r.q if verbose: tag = " [= reference point after m2-rescaling]" if is_ref else "" print(f"Laurent pole layers at (s,t,m2) = ({s_}, {t_}, {m2_}), i.e." f" (shat,that) = ({sh_r}, {th_r}){tag} (computed now from the" f" closed-form kernel; quad dps {qdps}, kernel dps {kdps}):") res = pole_layers(shat, that, is_ref, qdps, kdps, nfit=nfit_for(DPS)) out['I_em2'] = res['I_em2'] out['Bsm1'] = res['Bsm1'] out['I_em1'] = res['I_em1'] out['_pole'] = {'fit_chk_d': float(res['fit_chk']), 'qdps': qdps, 'kdps': kdps, 'nfit': nfit_for(DPS), 'shat': str(sh_r), 'that': str(th_r)} if verbose: mp.mp.dps = 60 print(f" [check] live closed-dilog K_2^(0) fit vs exact parametric 2-fold" f" K_2^(0): {res['fit_chk']:.1f} d") print(f" eps^-2: -K_2^(0) = {mp.nstr(res['I_em2'], 36)} (computed)") if is_ref: o2 = mp.mpf(ORC['laurent']['eps-2']) print(f" oracle = {mp.nstr(o2, 36)} (independent reference, 50 d)") print(f" agreement = {agree_d(res['I_em2'], o2):.2f} d" f" (archived: 45.76 d)") print(f" eps^-1: -K_2^(1)+B[s_-1]^(0) = {mp.nstr(res['I_em1'], 36)} (computed;" f" K_2^(1) = exact parametric threshold 2-fold)") if is_ref: o1 = mp.mpf(ORC['laurent']['eps-1']) print(f" oracle = {mp.nstr(o1, 36)} (independent reference, 50 d)") print(f" agreement = {agree_d(res['I_em1'], o1):.2f} d" f" (archived: 39.57 d)") if m2_ != 1: print(f" [m2 != 1] printed layers are of Ihat(s/m2,t/m2,1);" f" I_VP(s,t,m2) = m2^(3d/2-10) * Ihat, d=4-2eps (layer mixing by" f" -(4+3eps)*log(m2)).") print(" eps^0 Laurent layer: available via --eps0 (pole-subtracted multi-eps" " Neville extraction, computed live; archived quad-floor value 20.90 d).") return out def neville0(xs, ys): """Neville polynomial extrapolation of (xs, ys) to x=0.""" T = list(ys) n = len(xs) for m in range(1, n): for i in range(n - m): T[i] = (xs[i + m] * T[i] - xs[i] * T[i + 1]) / (xs[i + m] - xs[i]) return T[0] def eps0_extract(DPS, ks, verbose=True): """eps^0 Laurent layer at the reference point: I0 = Neville(eps->0) of [ I_disp(eps) - I^(-2)/eps^2 - I^(-1)/eps ], with the pole layers computed at LIFTED precision (LDPS-class knobs) so the subtraction noise sits far below the per-point I_disp accuracy. Grid eps = 2^-k; any rational eps in (0, 1/16] is legal input, so extraction points need no oracle -- choose them below (or at the deep end of) the independent reference grid for truncation control. Error bar = leave-one-out spread; final comparison vs the independent 50-d eps0 reference string. Depth-varied certification legs (vary VP_LMAX and --dps between invocations) are harness-level: run twice and compare I0 (the per-point Richardson pair and the fit-vs-fold layer check are the in-run depth checks).""" t0 = time.time() SF = mp.mpf(os.environ.get('VP_STEP_FRAC', '0.28')) N_TAY, LMAX = settings_for(DPS) n_jobs = int(os.environ.get('VP_JOBS', 0)) or min(len(ks), max(1, (os.cpu_count() or 4) // 2)) # (1) pole layers at lifted precision (derived fold + scaled-depth Bsm1) LDPS = max(160, DPS + 20) qdps_l, kdps_l, nfit_l = max(55, LDPS - 15), max(90, LDPS + 20), nfit_for(LDPS) if verbose: print(f"eps^0 extraction (Route N): grid eps = 2^-k, k in {ks}; transport dps={DPS}," f" Taylor {N_TAY}, quad levels {LMAX}/{LMAX+1}, {n_jobs} worker(s)") print(f" [1/3] pole layers at lifted precision (dps {LDPS}: qdps {qdps_l}," f" kdps {kdps_l}, nfit {nfit_l}) ...", flush=True) mp.mp.dps = kdps_l + 60 shat, that, _ = kin() res = pole_layers(shat, that, True, qdps_l, kdps_l, nfit=nfit_l) em2, em1 = res['I_em2'], res['I_em1'] o2 = mp.mpf(ORC['laurent']['eps-2']) o1 = mp.mpf(ORC['laurent']['eps-1']) if verbose: print(f" layer comparisons vs independent 50-d reference strings: eps^-2 {agree_d(em2, o2):.2f} d," f" eps^-1 {agree_d(em1, o1):.2f} d (string-capped); fit-vs-fold" f" {res['fit_chk']:.1f} d", flush=True) # (2) I_disp at each grid eps (forked workers, one per point) if verbose: print(f" [2/3] I_disp at {len(ks)} grid points ...", flush=True) pts = [sp.Rational(1, 2 ** k) for k in ks] def run_pt(ee): val, diag = I_disp_at_eps(ee, DATA['sigvp_de'], box_de(), DPS, N_TAY, SF, LMAX, verbose=False) return mp.nstr(val, DPS + 10), diag results = {} if n_jobs > 1 and len(pts) > 1: import multiprocessing ctx = multiprocessing.get_context('fork') def worker(i, conn): try: conn.send((i,) + run_pt(pts[i])) except Exception as ex: conn.send((i, None, {'error': repr(ex)})) conn.close() pending = list(range(len(pts))) pipes, procs, active = {}, {}, [] while pending or active: while pending and len(active) < n_jobs: i = pending.pop(0) pr, pw = ctx.Pipe(False) p = ctx.Process(target=worker, args=(i, pw)) p.start() pipes[i], procs[i] = pr, p active.append(i) i = active.pop(0) ii, vstr, diag = pipes[i].recv() procs[i].join() results[ii] = (vstr, diag) else: for i, ee in enumerate(pts): results[i] = run_pt(ee) mp.mp.dps = max(240, DPS + 80) xs, ys = [], [] for i, k in enumerate(ks): vstr, diag = results[i] if vstr is None: raise SystemExit(f"eps=2^-{k} FAILED: {diag.get('error')}") val = mp.mpf(vstr) line = (f" 2^-{k}: quad L{LMAX}/L{LMAX+1} agree {diag['quad_agree']:.1f} d," f" {diag['t_total']:.0f}s") ostr = ORC['I_AMF_grid'].get(str(k)) if ostr is not None: line += f" | vs independent reference value {agree_d(val, mp.mpf(ostr)):.2f} d" if verbose: print(line, flush=True) xs.append(mp.mpf(1) / 2 ** k) ys.append(val - em2 * mp.mpf(4) ** k - em1 * mp.mpf(2) ** k) # (3) Neville to eps=0 + leave-one-out spread + comparison with the independent reference value I0 = neville0(xs, ys) loos = [neville0(xs[:j] + xs[j + 1:], ys[:j] + ys[j + 1:]) for j in range(len(xs))] spread = max(abs(l - I0) for l in loos) o0 = mp.mpf(ORC['laurent']['eps0']) d0 = agree_d(I0, o0) if verbose: print(f" [3/3] Neville(eps->0) over {len(ks)} points:") print(f" eps^0 layer = {mp.nstr(I0, 40)} (computed now)") print(f" oracle = {mp.nstr(o0, 40)} (independent reference, 50 d)") print(f" agreement = {d0:.2f} d (archived quad-floor run: 20.90 d)") print(f" leave-one-out spread = {mp.nstr(spread, 3)}" f" (~{float(-mp.log10(spread / abs(I0))):.1f} d bar)") print(f" ({time.time()-t0:.0f}s total)") return {'I0': I0, 'gate_d': d0, 'loo_spread': spread} def _sha256(path): import hashlib with open(path, 'rb') as fh: return hashlib.sha256(fh.read()).hexdigest() def write_out(path, DPS, pts, out1, T0, dry=False): """--out json: the fixed-eps value(s) as full-dps strings (DPS+10 digits, the record's own serialization in eps0_extract) + diagnostics + provenance. The key I_disp is ABSENT on a dry run (a resume test on that key must fail on a dry file).""" import subprocess s_, t_, m2_ = KIN sh, th, uh = kin_hat() rec = {'record': 'I_DISP_FIXED_EPS (row 33 (s,t) point; build B1 evaluator)', 'status': 'DRY_RUN' if dry else 'DONE', 'point': {'s': str(s_), 't': str(t_), 'msq': str(m2_), 'u': str(-(s_ + t_)), 'shat': str(sh), 'that': str(th), 'uhat': str(uh), 'is_reference': bool(is_reference())}, 'dps': DPS, 'box_de': {'path': BOX_DE_PATH or f"{os.path.join(HERE, 'row33_data.json')}#box1eq_de", 'sha256': _sha256(BOX_DE_PATH) if BOX_DE_PATH else _sha256(os.path.join(HERE, 'row33_data.json')), 'masters': box_de()['masters'], 'n_samples': box_de().get('n_samples'), 'mismatches': box_de().get('mismatches')}, 'oracle': ({'source': ORACLE_PATH, 'sha256': _sha256(ORACLE_PATH), 'ks': sorted(ORACLE_GRID)} if ORACLE_GRID is not None else (None if (NO_ORACLE or not is_reference()) else {'source': 'DATA[oracles][I_AMF_grid] (reference literals)'}))} if pts: e = pts[0] rec['eps'] = str(e) rec['eps_label'] = eps_label(e) rec['k'] = (e.q.bit_length() - 1) if (e.p == 1 and (e.q & (e.q - 1)) == 0) else None if len(pts) > 1: rec['eps_all'] = [str(x) for x in pts] if dry: mp.mp.dps = DPS rec['b_sings_30d'] = [mp.nstr(x, 30) for x in de_singularities(box_de(), pts[0] if pts else sp.Rational(1, 64), DPS)] else: mp.mp.dps = DPS + 20 if pts: lab = eps_label(pts[0]) val = out1[f'I_disp({lab})'] inf = out1['_info'][lab] rec['I_disp'] = mp.nstr(val, DPS + 10) rec['I_disp_40d'] = mp.nstr(val, 40) rec['I_disp_is'] = 'Ihat(shat, that) at m2=1' + ('' if m2_ == 1 else '; physical I_VP(s,t,m2) = m2^(3d/2-10) * Ihat, d=4-2eps (module DOMAIN note)') if m2_ != 1: epsv = mp.mpf(pts[0].p) / pts[0].q m2v = mp.mpf(m2_.p) / m2_.q rec['I_disp_phys'] = mp.nstr(val * m2v ** (-4 - 3 * epsv), DPS + 10) rec['gate'] = {'oracle_value': inf['oracle'], 'agreement_d': inf['agreement_d']} d = dict(inf['diag']) for kk in ('tK', 'tR', 't_seed', 't_total', 'xchk', 'seed_selfchk', 'kn_selfchk', 'quad_agree'): if kk in d: d[kk] = float(d[kk]) rec['diag'] = d if len(pts) > 1: rec['I_disp_all'] = {eps_label(x): mp.nstr(out1[f'I_disp({eps_label(x)})'], DPS + 10) for x in pts} if 'I_em2' in out1: rec['pole_layers'] = {'I_em2': mp.nstr(out1['I_em2'], DPS), 'I_em1': mp.nstr(out1['I_em1'], DPS), 'Bsm1': mp.nstr(out1['Bsm1'], DPS), **out1.get('_pole', {})} rec['digits_note'] = ('the strings carry the working precision, not the accuracy; the accuracy is the ' 'comparison (agreement_d here when an oracle is consumed, else the recorded run)') rec['module_sha256'] = {f: _sha256(os.path.join(HERE, f)) for f in ('eval_row33.py', 'k33lib.py', 'lbl3disp.py', 'rho33lib.py', 'row33_data.json')} rec['wall_s'] = round(time.time() - T0, 2) rec['PRODUCER'] = {'script': os.path.abspath(__file__), 'sha256': _sha256(os.path.abspath(__file__)), 'argv': sys.argv[1:], 'host': os.uname().nodename, 'python': sys.version.split()[0], 'mpmath': mp.__version__, 'sympy': sp.__version__, 'stamp': subprocess.run(['date', '-u', '+%Y-%m-%dT%H:%M:%SZ'], capture_output=True, text=True).stdout.strip()} if path: tmp = path + '.tmp' with open(tmp, 'w') as fh: json.dump(rec, fh, indent=1) os.replace(tmp, path) return rec def main(): import argparse ap = argparse.ArgumentParser( description="row 33 (LBL3VP) standalone evaluator (ZERO AMFlow): fixed-eps " "dispersive value at the reference point (s,t,m2)=(-1,-1/3,1) compared with " "the independent 3-loop AMFlow reference grid (never consumed by the computation), + live eps^-2/eps^-1 pole layers at " "arbitrary Euclidean (s,t,m2). See module docstring for scope/gaps.") ap.add_argument('--dps', type=int, default=110, help="transport working precision in decimal digits (default 110 = the reference comparison run)") ap.add_argument('--point', default=None, metavar='P', help="one token, comma/space separated inside quotes. 1 value: eps " "('2^-k', 'p/q' or exact decimal) -> I_disp(eps) at the reference " "point (compared if on the k=4..18 grid). 3 values: s,t,m2 " "(rationals/decimals; s<0, t<0, m2>0, -s-t<4*m2) -> pole layers, " "e.g. --point=-2,-1/5,1") ap.add_argument('--check', action='store_true', help="rerun everything at dps+60 and print two-precision agreement") ap.add_argument('--eps0', action='store_true', help="compute the eps^0 Laurent layer at the reference point live. " "DEFAULT construction (2026-07-05 specification v2, " "ROW33_SUBTRACTION_SPEC): fully analytic plus-distribution " "subtraction at s0=7/5 -- c0 = exact Gamma-series terms + " "finite integrals under refine-until-bound quadrature " "(row33_eps0.py). No eps grid.") ap.add_argument('--legacy-grid', action='store_true', help="with --eps0: use the RETIRED multi-eps Neville grid " "extraction instead (provenance only; grid-class method, " "cert ceiling ~39 d)") ap.add_argument('--sabotage-k2', action='store_true', help="with --eps0: mutation test -- perturb K2^(0) by 1e-30; " "the comparison must FAIL with exit code 1") ap.add_argument('--eps0-ks', default=None, metavar='KS', help="comma list of grid exponents k (eps=2^-k) for " "--eps0 --legacy-grid (default: 21..26 at dps>=140, else 14..19)") ap.add_argument('--kin', default=None, metavar='S,T[,M2]', help="(build B1) kinematic point of the fixed-eps value, rationals " "(s<0, t<0, m2>0, -s-t<4*m2); default the reference -1,-1/3,1") ap.add_argument('--box-de', default=None, metavar='FILE', help="(build B1) the point's one-loop box w'-DE json (" "schema of DATA['box1eq_de']); required off the reference") ap.add_argument('--oracle-json', default=None, metavar='FILE', help="(build B1) AMFlow black_box output at the point (independent reference grid): " "compare the fixed-eps value live against result[0]") ap.add_argument('--no-oracle', action='store_true', help="(build B1) never consume an oracle: compute and write the value only") ap.add_argument('--out', default=None, metavar='FILE', help="(build B1) write {k, eps, I_disp , diag, ...} json") ap.add_argument('--dry-run', action='store_true', help="(build B1) parse, load the DE, derive its singularities, write --out " "without the key I_disp (status DRY_RUN); evaluate nothing") args = ap.parse_args() T0 = time.time() DPS = args.dps global BOX_DE, BOX_DE_PATH, NO_ORACLE, ORACLE_GRID, ORACLE_PATH if args.kin: toks = args.kin.replace(',', ' ').split() if len(toks) not in (2, 3): raise SystemExit("--kin takes s,t or s,t,m2") set_kin(*[sp.Rational(x) for x in toks]) ks_, kt_, km_ = KIN if not (ks_ < 0 and kt_ < 0 and km_ > 0 and -(ks_ + kt_) < 4 * km_): raise SystemExit(f"--kin (s,t,m2)=({ks_},{kt_},{km_}) outside the documented " "deep-Euclidean domain: need s<0, t<0, m2>0, u=-s-t<4*m2") if args.box_de: BOX_DE_PATH = os.path.abspath(args.box_de) BOX_DE = json.load(open(BOX_DE_PATH)) for key in ('masters', 'A'): if key not in BOX_DE: raise SystemExit(f"--box-de {BOX_DE_PATH}: missing key {key}") if len(BOX_DE['A']) != 8 or any(len(r) != 8 for r in BOX_DE['A']) or len(BOX_DE['masters']) != 8: raise SystemExit(f"--box-de {BOX_DE_PATH}: A must be 8x8 over 8 masters") if BOX_DE.get('var', 'wp') != 'wp' or BOX_DE.get('dim_var', 'd') != 'd': raise SystemExit(f"--box-de {BOX_DE_PATH}: var/dim_var must be wp/d") NO_ORACLE = bool(args.no_oracle) if args.oracle_json: ORACLE_PATH = os.path.abspath(args.oracle_json) ORACLE_GRID = load_amflow_grid(ORACLE_PATH) if args.eps0 and not is_reference(): raise SystemExit("--eps0 is reference-only (its comparison literals are the reference's)") if args.eps0: if args.legacy_grid: ks = ([int(x) for x in args.eps0_ks.split(',')] if args.eps0_ks else (list(range(21, 27)) if DPS >= 140 else list(range(14, 20)))) print(f"row 33 (LBL3VP) eps^0 Laurent layer (LEGACY grid mode, " f"provenance only); dps={DPS}") eps0_extract(DPS, ks, verbose=True) else: import row33_eps0 row33_eps0.run(DPS, sabotage=('k2' if args.sabotage_k2 else None)) print(f"\ntotal wall time: {time.time()-T0:.1f}s") return pts = [] kinpt = None if args.point is None: pts = [sp.Rational(1, 64)] # comparison point 2^-6 kinpt = KIN # layers at the point else: toks = args.point.replace(',', ' ').split() if len(toks) == 1: pts = [parse_eps(toks[0])] elif len(toks) == 3: kinpt = tuple(sp.Rational(x) for x in toks) else: raise SystemExit("--point takes 1 value (eps) or 3 values (s,t,m2); see --help") for e in pts: if not (0 < e <= sp.Rational(1, 16)): raise SystemExit(f"eps={e} outside the supported domain 0 < eps <= 1/16 " "(oracle grid side; see module docstring)") if kinpt is not None: s_, t_, m2_ = kinpt if not (s_ < 0 and t_ < 0 and m2_ > 0 and -(s_ + t_) < 4 * m2_): raise SystemExit(f"(s,t,m2)=({s_},{t_},{m2_}) outside the documented " "deep-Euclidean domain: need s<0, t<0, m2>0, u=-s-t<4*m2") if pts and not is_reference() and BOX_DE is None: raise SystemExit("fixed-eps value off the reference point needs --box-de (the point's " "one-loop box w'-DE)") if args.dry_run: write_out(args.out, DPS, pts, None, T0, dry=True) print(f"dry run: point {KIN}, DE {BOX_DE_PATH or 'DATA[box1eq_de]'}, eps {[eps_label(e) for e in pts]}," f" out {args.out} (no I_disp key); nothing evaluated") return print(f"row 33 (LBL3VP) standalone evaluator (ZERO AMFlow at eval time); dps={DPS}") print("(kernel seed + K_n = exact Feynman-parameter integrals; rho_VP symbolic;") print(" only inputs = exact-rational w-DEs + independent reference literals (comparison only), see header)\n") out1 = one_pass(DPS, pts, kinpt, verbose=True) if args.out: write_out(args.out, DPS, pts, out1, T0) if args.check: print(f"\n--check: re-running everything at dps={DPS + 60} ...") t1 = time.time() out2 = one_pass(DPS + 60, pts, kinpt, verbose=False) mp.mp.dps = DPS + 90 for k in sorted(out1): if k == 'gate_min_d' or k.startswith('_') or k not in out2: continue print(f" {k}: two-precision agreement {agree_d(out1[k], out2[k]):.2f} d") print(f" (--check leg wall: {time.time()-t1:.0f}s)") print(f"\ntotal wall time: {time.time()-T0:.1f}s") if __name__ == '__main__': main()