#!/usr/bin/env python3 r"""LBL3KP (QED kite parent: light-by-light box x 2-loop crossed self-energy, nu3=2 doubled stub) -- COMPLIANT final form, computed at runtime. Deps: python3 + mpmath (pip) + the sibling kite_de_exact_parse.py (the exact parser of the shipped connection; no sympy). THE FAMILY. 9-propagator 3-loop gamma-gamma->gamma-gamma family (V=8, E=10, planar, one closed fermion loop): the 1-loop equal-mass LbL box with one rung dressed by the 2-loop crossed-photon electron self-energy Sigma_kite, flanked by nu3 = 2, 1, 0 powers of the stub propagator (l^2-m^2). Three targets: I_KP = J[1,1,1,2,1,1,1,1,1] (the parent, nu3=2) I_K1 = J[1,1,1,1,1,1,1,1,1] (single stub, nu3=1) I_SE = J[1,1,1,0,1,1,1,1,1] (= LBL3SE, nu3=0) FINAL FORM (IBP-reduction of the parent onto kite-family masters + one-fold integrals with closed-form kernels -- no numeric node cache anywhere): I_nu3(s,t,m^2) = (1/pi) \int_{m^2}^{inf} dw K_nu3(s,t;w) * rho(w), [PF] The stub block is partial-fractioned EXACTLY ((l^2-m^2) and (l^2-w) share the momentum l): K_P(w) = [K(w) - K(m^2)]/(w-m^2)^2 - K'(m^2)/(w-m^2) (nu3=2) K_1(w) = [K(w) - K(m^2)]/(w-m^2) (nu3=1) K_0(w) = K(w) (nu3=0) all regular at w=m^2 (Taylor route there, raw PF outside, live seam checks). K(m^2), K'(m^2) are COMPUTED here from (s,t) by two independent routes each (direct / Chebyshev-fit / Richardson FD). [K] K(w) = Box1(s,t; m^2,m^2,m^2, M^2=w), the one-loop massive box in its CLOSED dilogarithmic 1-dim form (4-root partial fractions + one fixed Gauss-Legendre rule; every value computed here from (s,t,m2,w)). [rho] rho(w) = -Im J_top(w+i0), the discontinuity of the eps^0 kite-family TOP master J[1,1,1,1,1] -- the IBP-onto-kite-masters content. NOT a stored table: evaluated at runtime, to the requested precision, as a Picard-Fuchs / IBP-connection SERIES SOLUTION -- the exact rational 8x8 system d/dw J = A(w,d) J (lbl3kp-kite-de.json) is eps-graded and solved by adaptive local Taylor series along (1,inf), seeded once at w=5. The elliptic content is the equal-mass sunrise block of A (Gamma_1(6), cusps {0,1,9,inf}); the cusp expansions at w=1 (turn-on -2pi*(w-1)log(w-1), coefficient CLOSED FORM) and w=9 (finite jump) are that curve's Eichler-word data. The non-elliptic kite masters have Gamma/2F1 closed forms, VERIFIED live against the transported Laurent blocks (section [2]). [rho closed form] (2026-09-11) the same density written out in closed form -- the write-up's eqs. (rho_elem), (rho_full), (F23); m^2 = 1: rho(w) = -(pi/w) [2 ln(w-1) ln w + 3 Li2(1-w)] 1 < w <= 9 rho(w) = the same - (1/w) int_9^w v F23(v) dv w > 9 F23(v) = [2 Im S(v) - (v+3) Im Sdot(v)] / (v (v-1)^2), Im S / Im Sdot = the d=4 equal-mass sunrise three-body cut and the dotted-sunrise cut, one-folds over s in [4, (sqrt v - 1)^2] (the formulas in full in the block comment of section C2 below). THIS is the density the script evaluates by default at the density gate points w = 5, 12, 100, in every tier (the fast-start tier included: seconds), gated against the AMFlow solve_integrals references at the served floors; the series solution of [rho] (which the dispersion quadrature consumes -- the value path of the integrals is unchanged) is printed beside as the cross-check and the agreement of the two densities gated at the rho self-check bar. --density series prints the series density alone (the behaviour before 2026-09-11); --sum-rule checks the large-l^2 normalisation (1/pi) int_1^inf rho dw = 6 zeta_3 on the closed form at a finite cutoff W (default 1e50), the no-tail residual gated against the cutoff's own truncation (the omitted tail (4 ln W + 10)/W of the sum, 4.7e-48 at W = 1e50; the bar capped at the inner working precision) and the tail-restored residual against the inner working precision (a cutoff too low for a raised --dps, i.e. whose next omitted tail (4 ln W + 8)/W^2 does not clear that precision, is refused by name). [quad] singularity-subtracted tanh-sinh quadrature (Dispersify): the w=1 turn-on model is subtracted and added back as exact power-log moments x kernel-Taylor coefficients; panel break at the w=9 Gamma_1(6) cut; mapped real tail. rho values are memoized ACROSS the three ladders (runtime memo, nothing stored on disk). ARBITRARY PRECISION: no numeric node cache anywhere (the round-1 1669-node rho_kite_cache.json is DELETED). Precision is set by --dps (= QUAD_DPS): LEVEL/DPS/NORD and the truncation guards (U_MIN, V_MIN, W_MAX -- printed with derived bounds) rescale from it, so more quadrature levels / more series terms deliver more digits of the SAME analytic objects. SERVED PRECISION RANGE AND ITS CAP (2026-09-06, Q20f; the same closed-box form as the sibling lbl3se-evaluate.py, on whose kernel the figures were measured): the closed dilog-box kernel's far tail (BoxKernel._tail) relaxes the box quadrature to a floor of 18 digits (work 36) past w ~ 1e18; the four-root partial fraction then resolves 1/M5 only to (36 - log10 w) digits, so the kernel's tail accuracy is capped near 1e-62..1e-66 from dps ~62 on (2.2e-61 at dps 60, 3.8e-64 at dps 66, 1.6e-66 at dps 70, summed over the tail nodes), whatever --dps says; past w ~ 1e37.5 that floor returned NaN (log 0), which every --dps >= 72 reaches (W_MAX = 10^(dps//2+2)); cured 2026-09-06 by a dps floor of ceil(log10 w) + 12 at w >= 1e36 with the value asserted finite. The --point P1 tier at dps 60 (W_MAX = 1e32) stays below the floor's switch; the other documented tiers have W_MAX <= 1e32; digits past ~62 are not claimed at any --dps until a later cure lifts the relaxation below the onset. INPUT DATA (this dir; nothing else read at runtime): lbl3kp-kite-de.json exact rational 9x9 kite IBP connection A(w,d) (Kira; reconstructed from 43 mode:"diffeq" samples, max validation diff 0). Byte-copy (md5 cf204e14dc8fb1e804a7571b3a4427fd) of /phys_lbl3se/numeric/kite_de_symbolic.json lbl3kp-w5-derived.json the single transport SEED: eps-Laurent of the kite masters at w=5, DERIVED AMFlow-free (p^2=0 vacuum closed forms + Frobenius at w=0 + exact-DE Taylor march; emitted once at dps 80 -- rerun lbl3kp-w5-seed.py --dps N --json to regenerate, or pass --boundary-recompute [DPS] to THIS script to do that inline: the stored strings then demote to a byte-agreement fast-start cache). Byte-copy (md5 c9bfec1cdeb086222a5d771923fc21a7) of ../lbl3se/lbl3se-w5-derived.json. The pipeline is now free of AMFlow inputs; AMFlow artifacts appear only as held-out oracles. lbl3kp-kite-boundary-w5.json the RETIRED AMFlow w=5 seed (solve_integrals, goal_digits=140; md5 9e7bbca1e6ceb548a2322717de6f4797, byte-copy of /phys_lbl3se/numeric/amflow_kite_bnd_w5_out.json). Held-out cross-check only: the run prints the recomputed derived-vs-AMFlow seed agreement; it is never consumed upstream. HELD-OUT ORACLES (literals below, byte-traced by grep to banked artifacts, NEVER used in the computation; agreement RECOMPUTED live as -log10|f-o|/|o|): ORACLE_KP / ORACLE_K1 / ORACLE_SE Neville eps->0 extraction (kmin=0, the recorded gate.py form) of the deeper nine-propagator AMFlow black-box grid of THIS family at the reference point (2026-09-06): 20 eps samples 2^-4..2^-23, x_order 500 with an x_order-400 companion, working_pre 340; order-convergence bars 68.73/68.56/68.17 relative (69.68/69.23/68.45 absolute), leave-one-out bars 74.45/74.28/73.89 relative (75.40/74.95/74.17 absolute), x_order pair 93.52/93.99/98.23 d relative -- relative = -log10(|delta|/|value|), the paper's convention; absolute = -log10|delta|, the comparator of record's (2026-09-09); 60 digits stored. <- the grid outs by sha256: LBL3KP-REF_x500_n20.json 0514c51a24b70c5f990edf78564e7b036e357e13a44a6f6eb728131d0d69a5c7 LBL3KP-REF_x400_n20.json cc0c893f05a18e2a3c1bc0df8d431ee0abbdb3841ffcff00bdd565f303915b68 the compare receipt COMPARE_item33_REF_20260906T064534Z.json 27d7adc0cb4bc01b9754609d7d6c585b30bccb33ec08fb480eee9f5ef6b50262 The earlier grid (x_order 200, fifteen samples; conv 40.1/40.5/39.9 d, LOO 44.3/44.7/44.1 d) sits 39.09/39.55/43.42 d from these peels: the 2026-07 count of 39/39 was that grid's own floor. GT_AMF same LBL3SE number from the prior INDEPENDENT 8-prop-family AMFlow eps-grid run (42-digit artifact) <- /phys_lbl3se/amflow_v2/RESULT_epsgrid_eps0.json RHO_W12 rho(12) from an independent AMFlow solve_integrals run at w=12 (never used here; the transport seed is the w=5 run; 163-digit string -- the dps-growth exhibit) <- /phys_lbl3se/numeric/amflow_kite_bnd_w12_out.json RHO5 / RHO100 rho(5) and rho(100) from the AMFlow solve_integrals runs at w=5 (140 d; the density at the seed point) and w=100 -- the sibling lbl3se-evaluate.py's density spot references, carried verbatim with its floors (GATE_FLOORS); references for the closed-form density gate only (2026-09-11) <- /phys_lbl3se/numeric/amflow_kite_bnd_w5_out.json, amflow_kite_bnd_w100_out.json AMF_BOX_M / AMF_BOX_DOT archived AMFlow values of the equal-mass box nu=[1,1,1,1] and [2,1,1,1] eps^0 at the demo point (138-digit strings); validate the LIVE K(1), K'(1) only -- the quadrature uses the computed values <- /phys_lbl3_qed_parents/kp/amflow_box1_eq_out.json USAGE (evaluation interface -- kinematic point and precision caller-chosen): python3 lbl3kp-evaluate.py gate demo (default) python3 lbl3kp-evaluate.py --dps 60 same point, more digits: the three Neville rows read 46.08/ 46.01/46.59 d vs the 60-digit literals (this script's own digits at dps 60; measured wall 2125 s on the build host); rho(12) and the kernel constants keep growing python3 lbl3kp-evaluate.py --point -2 -0.2 I_*(s,t,m2=1) at a different kinematic point [+ --dps D] python3 lbl3kp-evaluate.py --no-fastcache force the full live run at the gate point (~10-40 min) python3 lbl3kp-evaluate.py --density series any tier with the density gate points reported by the series density alone (the closed form not evaluated: the prints of the version before 2026-09-11) python3 lbl3kp-evaluate.py --sum-rule [LOG10_W] the spot gate (1/pi) int_1^inf rho dw = 6 zeta_3 on the closed- form density at the cutoff W = 10^LOG10_W (default 50; accepted 34..240 where the next omitted tail (4 ln W + 8)/W^2 clears the inner working precision at the running --dps -- the whole range at the default dps -- else refused by name with the smallest LOG10_W that does; minutes, the wall printed): the no-tail residual gated against the cutoff's own truncation capped at the inner working precision, the tail-restored residual against the inner working precision, and rho(5), rho(12), rho(100) from the same grid against the AMFlow references; a tier of its own python3 lbl3kp-evaluate.py --point P1 the tagged point P1 = (-1/2,-1,1): I_KP, I_K1, I_SE computed live and GATED vs the deeper grid's peels in points/P1.json (at dps 60: agreements 45.17/46.09/45.65 d vs the tier bars 44/45/44 d; measured wall 2232 s on the build host at loadavg 94.77). The counts of record 50/50/59 are the recorded evaluations' (dps 45/45/55 vs the full- precision peel), printed beside; the same value path read 45.2/44.7/50.7 d (I_KP) at dps 60/70/80: the digits past ~45 are realisation-dependent here (the add-back kernel-Taylor fit, FAST-START CACHE below); at dps 80 the LEVEL-6 refinement's last difference ~1e-56 fails this tier's own test 10^-(dps-15) (rc 1), so dps 60 is the measured dps with the largest worst-target agreement among those passing it FAST-START CACHE (2026-07-06): lbl3kp-fastcache.json in this directory holds the gate-demo values (I_KP, I_K1, I_SE, rho(12), kernel constants) BANKED from one full live run of THIS script (provenance + sha pins inside the file). CERTIFIED DEPTH (2026-09-06): the stored I strings are this script's own dps-45 output (QUAD_DPS + 8 = 53 significant digits); the depth certified is the two-precision pair recorded in the cache header (certified_depth_d, the floor of the smallest recorded agreement in banked_d: 38 d for the 2026-09-06 bank at dps 45), the fast-path prints say so and the printed gate agreements are capped there, and no digit beyond it is claimed. Cause (one sentence): the even Gauss-Legendre kernel rule (blog commit 3115ef81, 2026-09-05) moved this script's ill-conditioned 92-term kernel Taylor fit and its 87-of-92 add-back, so digits past ~41 are realisation-dependent -- the 2026-07-07 bank and today's bytes differ there, and neither is accurate past the pair. Default behavior at the banked gate point for --dps <= the banked depth: the banked values print instantly, the Neville-oracle held-out gates are re-asserted NOW on the banked strings (any mismatch, incl. a 1e-30 mutation of a cached string, exits nonzero), and live cheap cross-checks run (w=5 seed held-out gate; the closed dilog-box kernel rebuilt live vs its banked K(1)). The cache is a fast start ONLY -- the live machinery in this file remains the definition: any --point, --dps above the banked depth, --boundary-recompute, --no-fastcache, or QUAD_DPS/LEVEL/DPS/ NORD/SF env override runs it unchanged. DOMAIN: deep-Euclidean s<0, t<0, u=-s-t<4 in units m2=1 (the w=1,9 cusps of rho fix the mass scale; other m2>0 by dimensional rescaling of s,t: each I_nu3 is homogeneous, I_KP ~ (m2)^{-4-3eps} etc.). Off the demo point no banked oracle exists -- the script prints ladder self-consistency plus the still-held-out rho(12) gate and the kernel-constant cross-route checks. Default settings (QUAD_DPS=40, LEVEL=5, DPS=55, NORD=80): single core, wall time measured and printed at the end. The archived campaign record for this page (/phys_lbl3_qed_parents/RESULT.md): gate 40.04d (nu3=2), 40.22d (nu3=1), 43.32d (nu3=0 cross-check) -- quoted as history, not measured here. """ import argparse import bisect import hashlib import importlib.util import json import math import os import sys import time import mpmath as mp HERE = os.path.dirname(os.path.abspath(__file__)) T0 = time.time() _ap = argparse.ArgumentParser( description="LBL3KP compliant final-form evaluator (see module docstring; " "DOMAIN: s<0, t<0, -s-t<4, units m2=1)", epilog="tiers: (no flags) gate demo served from the sha-pinned fast-start " "cache, seconds; --no-fastcache the live gate run at dps 40, " "~15 min; --point S T a bare point (ungated), the live wall; " "--point P1 the tagged point gated vs the deeper grid's peels in points/P1.json (at dps 60; measured wall 2232 s on the build host); " "--boundary-recompute [DPS] the live run seeded by the shipped " "seed shim; --density series any tier with the series density alone at the density gate points (the closed form, evaluated and gated there by default since 2026-09-11, switched off); --sum-rule [LOG10_W] the 6 zeta_3 sum-rule spot gate on the closed-form density, a tier of its own (minutes). Walls measured on a shared 96-core host.") _ap.add_argument('--point', nargs='+', metavar='S_T_or_TAG', default=None, help="evaluate I_*(s,t,m2=1) at this deep-Euclidean point " "(two numbers S T; ungated) instead of only the gate " "demo -- or ONE tag naming a shipped reference record: " "--point P1 = (-1/2,-1,1), gated vs points/P1.json at " "its recorded dps (60 unless --dps is given)") _ap.add_argument('--dps', type=int, default=None, help="target quadrature digits (default 40 = gate demo); " "LEVEL/DPS/NORD/guards are re-derived unless set by env") _ap.add_argument('--boundary-recompute', nargs='?', const=80, type=int, metavar='DPS', default=None, help="re-derive the w=5 seed by running the shipped " "lbl3kp-w5-seed.py shim (-> ../lbl3se/lbl3se-w5-seed.py) " "at DPS (default 80) instead of reading the cached " "strings; the stored JSON demotes to a byte-agreement " "cache check (2026-07-05 wiring, the LBL3KP build spec)") _ap.add_argument('--no-fastcache', action='store_true', help="skip the sha-pinned fast-start cache " "(lbl3kp-fastcache.json) and run the live machinery " "even at the banked gate point; the cache is a " "fast-start convenience ONLY -- the live machinery in " "this file remains the definition") _ap.add_argument('--fastcache-bank', action='store_true', help="maintainer mode: run fully live and, on a PASSING " "gate run, write/refresh lbl3kp-fastcache.json " "(banked values + measured agreements + sha pins)") _ap.add_argument('--density', choices=('closed', 'series'), default='closed', help="which density reports the density gate points w = 5, 12, 100 " "(2026-09-11): closed (default) = the closed form (weight-two " "dilogarithms below w = 9, the one-fold over the equal-mass " "sunrise cut above) evaluated live in every tier and gated " "against the AMFlow references, the series density printed " "beside as the cross-check; series = the Picard-Fuchs series " "density alone (the closed form not evaluated). The dispersion " "integrand is the series density under both settings") _ap.add_argument('--sum-rule', nargs='?', const=50, type=int, metavar='LOG10_W', default=None, help="the spot gate (1/pi) int_1^inf rho dw = 6 zeta_3 on the closed-form " "density at the cutoff W = 10^LOG10_W (default 50; accepted 34..240 where " "the next omitted tail (4 ln W + 8)/W^2 clears the tail-restored bar at the " "running --dps -- the whole range at the default dps; a cutoff too low for a " "raised --dps is refused by name with the smallest LOG10_W that passes): " "a tier of its own (minutes), rc 0 PASS / 1 FAIL / 2 usage or refused cutoff") _ARGS, _ = _ap.parse_known_args() # sha256 of the pinned sibling module (computed by the producer that cut this file, never typed); a byte change # or a missing file is REFUSED by name (exit 3) before any value is printed SCRIPT = "lbl3kp-evaluate" PINS = { "kite_de_exact_parse.py": "0cc967d6362b63963b2e7ddd3e41f9edc1a4a0cab9cd2d7ebd7aef45aa83d7c1", } def _pinned_path(name): """Path of a sibling file; a PINNED file (the exact parser module) is refused on any byte change (exit 3, the recorded and recomputed sha256 named).""" path = os.path.join(HERE, name) want = PINS.get(name) if want is not None: if not os.path.exists(path): sys.stderr.write(f"{SCRIPT} REFUSED: pinned file {name} is missing (recorded sha256 {want})\n") raise SystemExit(3) have = hashlib.sha256(open(path, "rb").read()).hexdigest() if have != want: pos = next((k + 1 for k, (x, y) in enumerate(zip(want, have)) if x != y), 0) sys.stderr.write(f"{SCRIPT} REFUSED: {name} integrity pin mismatch (recorded {want}, recomputed " f"{have}; first differing hex position {pos} of 64, 1-based) -- the shipped file " f"was altered\n") raise SystemExit(3) return path def _import_pinned_module(name): """Import a sibling module through the pin check (a byte change is REFUSED, exit 3).""" spec = importlib.util.spec_from_file_location(name[:-3], _pinned_path(name)) mod = importlib.util.module_from_spec(spec) spec.loader.exec_module(mod) return mod kx = _import_pinned_module("kite_de_exact_parse.py") # the exact parser of the connection strings (no sympy) # ---- tagged point tier (2026-09-05): --point P1 ------------------------------ # A tag names a shipped reference record points/.json in this directory: # the point (s,t,m2=1), the tier's default dps and the reference value strings # (recorded by an independent implementation of the same one-fold and gated # there against an independent AMFlow eps-grid; provenance and shas inside the # record). Unlike a bare `--point S T` (ungated), the tagged tier compares the # value(s) computed HERE against the record and exits 1 by name below the # record's bar. The record carries an integrity pin over its own body; a # record that fails its pin is REFUSED (rc 3). POINT_TAGS = {'P1': 'points/P1.json'} def _point_refuse(msg): import sys print(f"[point-tag] REFUSED: {msg} (rc 3)", file=sys.stderr) raise SystemExit(3) def _load_point_record(tag): import hashlib rel = POINT_TAGS[tag] path = os.path.join(HERE, rel) if not os.path.exists(path): _point_refuse(f"reference record {rel} for --point {tag} is missing from " "this directory") with open(path) as fh: rec = json.load(fh) body = {k: v for k, v in rec.items() if k != '_sha256_body'} want = rec.get('_sha256_body', '') got = hashlib.sha256(json.dumps(body, sort_keys=True, separators=(',', ':'), ensure_ascii=False).encode()).hexdigest() if got != want: _point_refuse(f"{rel} integrity pin mismatch (recorded {want}, " f"recomputed {got}; first differing hex position " f"{1 + next((i for i, (a, b) in enumerate(zip(want, got)) if a != b), min(len(want), len(got)))} " f"of 64, 1-based) -- the reference record was altered; " f"not serving --point {tag}") if rec.get('tag') != tag: _point_refuse(f"{rel} carries tag {rec.get('tag')!r}, not {tag!r}") return rec def _resolve_point_tag(args): """--point P1 -> (tag, record) with args.point rewritten to the record's decimal (s, t); --point S T -> (None, None), the ungated point mode as before; anything else -> usage error.""" if args.point is None: return None, None if len(args.point) == 1: tag = args.point[0] if tag not in POINT_TAGS: _ap.error(f"--point {tag}: unknown point tag (known: " f"{', '.join(sorted(POINT_TAGS))}); a bare point is given as " "--point S T (m2 = 1)") rec = _load_point_record(tag) args.point = [rec['point']['s_decimal'], rec['point']['t_decimal']] return tag, rec if len(args.point) != 2: _ap.error("--point takes S T (a bare point, m2 = 1) or one tag " f"({', '.join(sorted(POINT_TAGS))})") return None, None _POINT_TAG, _POINT_REC = _resolve_point_tag(_ARGS) QUAD_DPS = _ARGS.dps if _ARGS.dps else ( int(_POINT_REC['tier_dps']) if _POINT_REC is not None # the tag's recorded dps else int(os.environ.get("QUAD_DPS", "40"))) # tanh-sinh digits roughly double per level: 5 covers the 40d gate, +1/doubling LEVEL = int(os.environ.get( "LEVEL", str(5 + max(0, math.ceil(math.log2(QUAD_DPS / 40.0)))))) DPS = int(os.environ.get("DPS", str(QUAD_DPS + 15))) # rho-transport digits # Taylor order per transport step: per-step truncation ~ SF^NORD NORD = int(os.environ.get("NORD", str(max(80, int(1.45 * DPS) + 1)))) SF = mp.mpf(os.environ.get("SF", "0.25")) # step = SF * dist-to-singularity # truncation guards, DERIVED from QUAD_DPS (bounds printed in main): # U_MIN : below w=1+U_MIN the subtracted remainder (rho - A1 u log u) ~ B*u, # |B|<~10 -> dropped < 10*U_MIN^2*max|K_nu3(1)|/pi # V_MIN : within V_MIN of the w=9 finite-jump cusp rho is frozen at the edge # value; |rho'| ~ |log V_MIN| -> dropped < |K(9)|*V_MIN^2*log/pi # W_MAX : tail truncation; rho ~ 4pi log w/w^2 and K_P ~ K_1/w ~ K/w^2 decay # faster than the nu3=0 case bounded in the lbl3se sibling analysis U_MIN = mp.mpf(10) ** -(QUAD_DPS // 2 + 2) V_MIN = mp.mpf(10) ** -(QUAD_DPS // 2 + 1) W_MAX = mp.mpf(10) ** (QUAD_DPS // 2 + 2) # ---------------------------------------------------------------- oracles --- # Held-out literals (provenance in the docstring). Parsed at high dps inside # main() -- module-level mp.mpf parse at dps=15 is the trap that once # capped the sibling LBL3SE gate at 17.6 digits. ORACLE_KP_STR = "0.112062481433958330568395694749093143352780028985830100546698" ORACLE_K1_STR = "-0.216286487952560629274051222247763446168074545743523831915117" ORACLE_SE_STR = "0.524957776781144632332996415282646042436170972240670568090275" GT_AMF_STR = "0.524957776781144632332996415282646042436171" RHO_W12_STR = ("0.39118160908247272140297662324912768960644354162754176827840" "05818116281784996731939692480451865506078571628503177960305228" "830215624118821245684483066354860486343504") AMF_BOX_M_STR = ("0.17805022679323064388457145940884606798964278254995890951970672224460" "8363756215261852028001046479803625777795740996738835167421365134266104") AMF_BOX_DOT_STR = ("-0.08182371131729398217126910224859083105346020966938410581877343661644" "37153710261103325092253533222644247850151209999369312290789556836833") # 2026-09-11: the sibling lbl3se-evaluate.py's two other density spot # references (the AMFlow solve_integrals runs at w = 5 and w = 100), carried # verbatim with its floors keyed on the --dps target (QUAD_DPS); used by the # closed-form density gate of section C2 only -- the series density's gates # in this file are unchanged. RHO5_STR = ("1.663479584361080173948716060213024634373432582368693673876490533944980" "63819387518116473235729360240114759165722492372800189805496698030269") RHO100_STR = ("0.0077916795514592605197494159377621432968434210826071616597299343160" "8722568461946565680546674989984783444333061578558621081701943628764877") GATE_FLOORS = { 'rho(5)': lambda dps: min(float(dps), 79.0), 'rho(100)': lambda dps: min(0.45 * dps + 14.0, 40.0), } def agree_digits(a, b): d = mp.fabs(a - b) if d == 0: return mp.inf return -mp.log10(d / mp.fabs(b)) def point_tag_gate(tag, rec, values, ok_run): """The tagged tier's gate (2026-09-05; the deeper grid 2026-09-06): every target in the shipped record points/.json vs the value computed here, agreement = -log10|a-b|/|b| with the record strings parsed at the ambient precision. Two references per target: the deeper independent AMFlow grid's Neville peel (grid_peel_eps0_60d, 59-60 stored digits: the GATE, PASS iff every target reaches the record's tier_bar_d AND the run's own certificates held) and the recorded evaluator value (value: REPORTED, its own two-precision agreement named as that string's ceiling). Prints the computed and recorded strings side by side at min(62, dps+2) digits plus the recorded provenance; returns ok.""" mp.mp.dps = 90 n = min(62, QUAD_DPS + 2) pt = rec['point'] print(f"[P] [{time.time()-T0:6.1f}s] tagged point {tag}: (s,t,m2) = ({pt['s']},{pt['t']},{pt['m2']}) -- " f"computed here vs the shipped reference record {POINT_TAGS[tag]}") src = rec.get('sources', {}) ev = src.get('evaluator', {}) g5, g4 = src.get('independent_grid_x500', {}), src.get('independent_grid_x400', {}) print(f" reference: {rec.get('record', '')}") print(f" the deeper grid: {g5.get('name', '?')} {g5.get('sha256', '?')[:16]} (x_order {g5.get('x_order', '?')}) " f"with {g4.get('name', '?')} {g4.get('sha256', '?')[:16]} (x_order {g4.get('x_order', '?')}); " f"{g5.get('n_eps', '?')} eps samples {g5.get('eps_samples', '?')}; compare receipt " f"{src.get('compare_receipt', {}).get('sha256', '?')[:16]}") print(f" recorded values by {ev.get('name', '?')} {ev.get('sha256', '?')[:16]} + {ev.get('library', '?')} " f"{ev.get('library_sha256', '?')[:16]} (I_SE: {src.get('evaluator_row28', {}).get('name', '?')} " f"{src.get('evaluator_row28', {}).get('sha256', '?')[:16]})") if QUAD_DPS < int(rec['tier_dps']): print(f" NOTE: --dps {QUAD_DPS} is below the tier's recorded dps {rec['tier_dps']}; " "the bar is the recorded claim and may not be reachable here") ok = ok_run for name, tgt in rec['targets'].items(): if name not in values: continue ours = values[name] peel = mp.mpf(tgt['grid_peel_eps0_60d']) ref = mp.mpf(tgt['value']) d_peel = agree_digits(ours, peel) d_ref = agree_digits(ours, ref) bar = float(tgt['tier_bar_d']) bar_rec = float(tgt['bar_d']) two = float(tgt['check']['two_precision_agreement_d']) stored = _sig_digits(tgt['grid_peel_eps0_60d']) print(f" {name:5s} computed = {mp.nstr(ours, n)}") print(f" {name:5s} deeper-grid peel = {mp.nstr(peel, n)} ({stored} stored digits; " f"order-convergence bar {float(tgt['grid_peel_conv_d_relative']):.2f} relative ({float(tgt['grid_peel_conv_d_absolute']):.2f} absolute), " f"leave-one-out bar {float(tgt['grid_peel_loo_d_relative']):.2f} relative ({float(tgt['grid_peel_loo_d_absolute']):.2f} absolute), " f"x_order pair {float(tgt['x_order_pair_400_vs_500_d']):.2f} d relative; the paper prints the relative convention)") print(f" {name:5s} agreement (peel) = {mp.nstr(d_peel, 6)} d >= tier bar {bar:.0f} d " f"(the count of record {bar_rec:.0f} d: {tgt['bar_d_member']}) -- {'OK' if d_peel >= bar else 'FAIL'}") print(f" {name:5s} recorded value = {mp.nstr(ref, n)} (recorded at dps {tgt['value_dps']}; " f"two-precision {two:.2f} d = that string's own depth, the ceiling of the next line)") print(f" {name:5s} agreement (recorded value) = {mp.nstr(d_ref, 6)} d (reported)") ok = ok and (d_peel >= bar) return ok def _sig_digits(lit): return len(lit.lstrip('-+').replace('.', '').lstrip('0')) # ============== fast-start cache (2026-07-06) ================================ # Doctrine: cache = sha-pinned + compare-gated fast start; the live machinery # in this file remains the DEFINITION. See module docstring, FAST-START CACHE. FASTCACHE_PATH = os.path.join(HERE, 'lbl3kp-fastcache.json') _FC_GATE_KEY = 's=-1,t=-1/3' _FC_ENV_OVERRIDES = ('QUAD_DPS', 'LEVEL', 'DPS', 'NORD', 'SF') _FC_PIN_FILES = ('lbl3kp-evaluate.py', 'lbl3kp-kite-de.json', 'lbl3kp-w5-derived.json', 'lbl3kp-kite-boundary-w5.json') def _fc_sha(fname): import hashlib with open(os.path.join(HERE, fname), 'rb') as fh: return hashlib.sha256(fh.read()).hexdigest() def _fc_str_sha(s): import hashlib return hashlib.sha256(s.encode()).hexdigest() def _fc_eligible(): live_note = ("live machinery run (the definition) -- expect the full " "wall (~10-20 min at default --dps 40; more above)") if _ARGS.no_fastcache or _ARGS.fastcache_bank: return None, None if _ARGS.point is not None: return None, f"requested --point is not the banked gate point; {live_note}" if _ARGS.boundary_recompute is not None: return None, f"--boundary-recompute forces the live path; {live_note}" envs = [k for k in _FC_ENV_OVERRIDES if k in os.environ] if envs: return None, f"env override {envs} set; {live_note}" if not os.path.exists(FASTCACHE_PATH): return None, f"no lbl3kp-fastcache.json; {live_note}" try: with open(FASTCACHE_PATH) as fh: fc = json.load(fh) pt = fc['points'][_FC_GATE_KEY] except Exception as ex: return None, f"fast cache unreadable ({ex!r}); {live_note}" if QUAD_DPS > pt['cached_dps']: return None, (f"requested --dps {QUAD_DPS} exceeds the banked depth " f"{pt['cached_dps']}; {live_note}") return fc, None def _fc_fail(msg): raise RuntimeError(f"[fast-cache] COMPARE GATE FAIL: {msg} -- refusing " "to serve the cache; rerun with --no-fastcache for " "the live machinery (the definition)") def _fc_gate(label, d_now, d_banked, floor, cap=None): ok = (float(d_now) >= floor) and (abs(float(d_now) - float(d_banked)) <= 0.05) # 2026-09-06: the PRINTED agreement is min(measured, certified_depth_d), the # cap named; the gate rule above compares the measured figure, unchanged shown = float(d_now) if cap is None else min(float(d_now), float(cap)) capnote = ("" if cap is None else f" [capped at the certified depth {cap} d; measured {float(d_now):.4f} d]") print(f" [fast-cache gate] {label}: {shown:.4f} d{capnote} " f"(banked {d_banked:.4f} d, floor {floor} d) -- " f"{'OK' if ok else 'FAIL'}") if not ok: _fc_fail(f"{label}: recomputed {float(d_now):.4f} d vs banked " f"{d_banked:.4f} d (floor {floor} d)") def _fc_fast_path(fc): pt = fc['points'][_FC_GATE_KEY] prov = fc['_provenance'] cd = fc.get('certified_depth_d') # 2026-09-06: the header's certified depth cdtxt = (f"certified depth {cd} d by the two-precision pair (cache header " "certified_depth_d); the printed gate agreements are capped there" if cd is not None else "no certified_depth_d in the cache header " "(an earlier bank block): refused below") mp.mp.dps = 90 print("LBL3KP FAST-CACHE mode -- banked gate-demo values (sha-pinned, " f"compare-gated); {cdtxt}.") print(f" cache: {os.path.basename(FASTCACHE_PATH)} banked " f"{prov['generated_utc']} by: {prov['generator_cmd']}") print(f" banked live wall {prov['wall_s']:.0f}s; this fast start serves " f"--dps <= {pt['cached_dps']} (requested {QUAD_DPS})") print(" cache = fast start ONLY; the live machinery in this file is the " "definition (--no-fastcache)\n") # --- [p1] integrity pins -------------------------------------------- for fname, want in fc['_pins'].items(): if _fc_sha(fname) != want: _fc_fail(f"sha256 pin mismatch on {fname} (definition changed " "since banking; cache is STALE)") for key, want in pt['_string_pins'].items(): if _fc_str_sha(pt[key]) != want: _fc_fail(f"cached string '{key}' fails its sha pin (mutated)") print(f"[p1] integrity: {len(fc['_pins'])} file pins + " f"{len(pt['_string_pins'])} cached-string pins OK") if cd is None: # 2026-09-06: no certified depth recorded -> fail closed _fc_fail("cache header lacks certified_depth_d (written by an earlier " "bank block); cache is STALE") # 2026-09-06 (Q20f): the header's certified depth is never TRUSTED -- it is # re-derived from the banked_d rows the compare gates below re-verify, and a # header that disagrees with its own rows is refused by name (the same # refusal path as the header-less case) cd_rows = int(math.floor(min(float(v) for v in pt['banked_d'].values()))) if cd != cd_rows: _fc_fail(f"cache header certified_depth_d {cd} != floor(min banked_d) " f"{cd_rows} (the header disagrees with the cache's own rows); " "cache is STALE") # --- [p2] Neville-oracle held-out gates recomputed NOW ---------------- print("[p2] held-out oracle gates recomputed NOW from the banked strings " "vs the in-script literals:") rows = [ ("I_KP vs independent 9-prop AMFlow (Neville)", 'I_KP', ORACLE_KP_STR), ("I_K1 vs independent 9-prop AMFlow (Neville)", 'I_K1', ORACLE_K1_STR), ("I_SE vs independent 9-prop AMFlow (Neville)", 'I_SE', ORACLE_SE_STR), ("I_SE vs prior 8-prop-family AMFlow eps-grid ", 'I_SE', GT_AMF_STR), ] for i, (label, key, lit) in enumerate(rows): v = mp.mpf(pt[key]) if i < 3: print(f" {key} = {mp.nstr(v, 44)} (banked; certified depth " f"{cd} d by the two-precision pair; {len(pt[key])} stored " "characters)") _fc_gate(label, agree_digits(v, mp.mpf(lit)), pt['banked_d'][f'row{i}'], 30, cap=cd) r12_c = mp.mpf(pt['rho12']) print(f" rho(12) = {mp.nstr(r12_c, 50)} (banked)") _fc_gate("rho(12) vs held-out AMFlow w=12", agree_digits(r12_c, mp.mpf(RHO_W12_STR)), pt['banked_d']['d12'], min(QUAD_DPS - 10, 30), cap=cd) # --- [p3] live cheap cross-checks ------------------------------------- print("[p3] live cheap cross-checks (machinery exercised NOW):") t0 = time.time() funcs, masters = load_kite_de() M0, d_seed, seed_dps = load_kite_boundary(masters) bar_seed = min(seed_dps, 128) - 12 print(f" live w=5 seed held-out gate: {mp.nstr(d_seed, 5)} d >= " f"bar {bar_seed} d ({time.time()-t0:.1f}s) -- RAISES on fail") if not d_seed >= bar_seed: _fc_fail(f"live seed held-out gate {mp.nstr(d_seed, 5)} d < {bar_seed} d") t0 = time.time() kdps_live = 28 Kl = BoxKernel(mp.mpf(-1), mp.mpf(-1) / 3, mp.mpf(1), kdps_live, verbose=False) # verify() bars against the GLOBAL # QUAD_DPS; the fast path gates the rebuild itself, right below dK = agree_digits(Kl.tay1[0], mp.mpf(pt['Box_m'])) dKa = agree_digits(Kl.tay1[0], mp.mpf(AMF_BOX_M_STR)) bar_K = kdps_live - 4 print(f" live dilog-box kernel rebuild (dps {kdps_live}): K(1) vs " f"banked constant {mp.nstr(dK, 5)} d, vs archived AMFlow " f"nu=[1,1,1,1] {mp.nstr(dKa, 5)} d >= bar {bar_K} d " f"({time.time()-t0:.1f}s)") if not (dK >= bar_K and dKa >= bar_K): _fc_fail(f"live kernel rebuild vs archived K(1): " f"{mp.nstr(dK, 5)}/{mp.nstr(dKa, 5)} d < bar {bar_K} d") # --- [p4] the closed-form density at the density gate points (2026-09-11, section C2) --- if _ARGS.density == 'closed': print("[p4] the closed-form density evaluated NOW at the density gate points " "(the cached rho(12) series string of [p2] is the cross-check):") closed_density_gates({'rho(12)': pt['rho12']}, f"the cached rho(12) string of the --dps {pt['cached_dps']} run above", bar12=min(QUAD_DPS - 10, 30), bar_series=min(QUAD_DPS - 10, 30), indent=' ') print(f"\ntotal wall time: {time.time()-T0:.1f}s") print("PASS (fast-cache) -- banked values re-gated vs all held-out " "Neville/eps-grid oracles + live seed/kernel cross-checks; for " "the full live run use --no-fastcache") raise SystemExit(0) def _fc_bank(payload, wall_s): import datetime pt = dict(payload) pt['_string_pins'] = {k: _fc_str_sha(pt[k]) for k in ('I_KP', 'I_K1', 'I_SE', 'rho12', 'Box_m', 'Box_dot')} fc = { '_doc': ("fast-start cache for lbl3kp-evaluate.py: banked from ONE " "full live run (command below). The live machinery is the " "definition; this file only fast-starts the banked gate " "demo and is refused on any sha-pin or compare-gate " "mismatch. Regenerate: --fastcache-bank."), '_provenance': { 'generated_utc': datetime.datetime.now( datetime.timezone.utc).isoformat(), 'generator_cmd': 'python3 lbl3kp-evaluate.py --dps ' f'{QUAD_DPS} --fastcache-bank', 'settings': {'QUAD_DPS': QUAD_DPS, 'LEVEL': LEVEL, 'DPS': DPS, 'NORD': NORD}, 'wall_s': wall_s, 'source': '/axis3_wave/disp-fastcache/', }, # 2026-09-06: the certified depth = floor of the smallest recorded agreement # in banked_d (the two-precision pair the cache itself records); read, never typed 'certified_depth_d': int(math.floor(min(float(v) for v in pt['banked_d'].values()))), 'depth_note': (f"the stored strings are the evaluator's own dps-{QUAD_DPS} output " f"(mp.nstr at QUAD_DPS + 8 = {QUAD_DPS + 8} significant digits); " "certified_depth_d is the floor of the smallest agreement in " "banked_d (the two-precision pair the cache records, 2026-09-06); " "no digit beyond it is claimed"), '_pins': {f: _fc_sha(f) for f in _FC_PIN_FILES}, 'points': {_FC_GATE_KEY: pt}, } with open(FASTCACHE_PATH, 'w') as fh: json.dump(fc, fh, indent=1) print(f"[fast-cache] BANKED {FASTCACHE_PATH} (dps {pt['cached_dps']}; certified " f"depth {fc['certified_depth_d']} d by the two-precision pair; " f"{len(fc['_pins'])} file pins)") # ======================= A. Box1 kernel: closed dilogarithmic form =========== # One-loop massive box, four on-shell massless legs, masses (M5,m2,m2,m2), # pySecDec normalization. Cheng-Wu + degenerate (x1,x2) quadratic reduce the # Feynman parametrisation to a single dv integral of a 4-root partial-fraction # log sum g(v); complex singularities of g sit at v=-1 and |v|=1, independent # of M5, so one fixed Gauss-Legendre rule converges uniformly in M5=w. # (Port of /phys_lbl3se/numeric/box1_dilog.py, verified >=50d there.) _GL_CACHE = {} def _gl(n, work): key = (n, work) xw = _GL_CACHE.get(key) if xw is None: old = mp.mp.dps mp.mp.dps = work xw = mp.gauss_quadrature(n, 'legendre') mp.mp.dps = old _GL_CACHE[key] = xw return xw def _g_of_v(v, s, u, m2, M5, a): vp1 = v + 1 p1 = a * v + (M5 + m2) r1 = m2 * vp1 * vp1 p2 = (M5 + m2) * vp1 r2 = r1 - u * v sd1 = mp.sqrt(p1 * p1 - 4 * M5 * r1) sd2 = mp.sqrt(p2 * p2 - 4 * M5 * r2) twoM5 = 2 * M5 mal = (p1 - sd1) / twoM5; malp = (p1 + sd1) / twoM5 mbe = (p2 - sd2) / twoM5; mbep = (p2 + sd2) / twoM5 Ral = -1 / ((u + s * mal) * sd1); Ralp = 1 / ((u + s * malp) * sd1) Rbe = 1 / ((u + s * mbe) * sd2); Rbep = -1 / ((u + s * mbep) * sd2) return -(Ral * mp.log(mal) + Ralp * mp.log(malp) + Rbe * mp.log(mbe) + Rbep * mp.log(mbep)) def box1(s, t, m2, M5, dps=40): """eps^0 of the one-loop massive box, >= dps digits, computed at runtime. Deep-Euclidean region: s<0, t<0, m2>0, M5>0, u=-s-t < 4 m2.""" work = dps + 18 old_dps = mp.mp.dps mp.mp.dps = work s = mp.mpf(s); t = mp.mpf(t); m2 = mp.mpf(m2); M5 = mp.mpf(M5) u = -s - t a = M5 + m2 - s if not (s < 0 and t < 0 and m2 > 0 and M5 > 0) or u >= 4 * m2: mp.mp.dps = old_dps raise ValueError("box1: need s<0, t<0, m2>0, M5>0, u=-s-t<4m2") n = work + 12 n += n % 2 # even rule: an odd rule places a node at the midpoint v = 2 exactly, where # (u + s*malp) and (u + s*mbep) vanish together as M5 -> inf at kinematic # points with -u/s = 3 (t = 2s), a removable singularity that the node would # divide by; an even rule has no node there. nodes, weights = _gl(n, work) L = mp.mpf(2) one = mp.mpf(1) tot = mp.mpf(0) for tk, wk in zip(nodes, weights): v = L * (one + tk) / (one - tk) dv = L * 2 / (one - tk) ** 2 tot += wk * _g_of_v(v, s, u, m2, M5, a) * dv mp.mp.dps = old_dps return +tot # ---- fast kernel evaluation: runtime-built interpolants of the CLOSED form. # K is analytic on (0,inf) (Euclidean s,t); on each finite region a Chebyshev # interpolant (built and VERIFIED at runtime, degree set by QUAD_DPS via the # Bernstein-ellipse rate to the nearest singularity w=0) evaluates it fast. # This is an evaluation device for the closed form, not stored data. class BoxKernel: def __init__(self, s, t, m2, qd, taylor_nmax=None, verbose=True): self.s, self.t, self.m2 = s, t, m2 self.kd = qd + 8 # kernel target digits self.nmax = taylor_nmax or int(1.9 * qd) + 6 self._direct_cache = {} old = mp.mp.dps mp.mp.dps = self.kd + 20 # region A: Taylor about w=1 on [0.62,1.38] (also the addback coeffs) self.tay1 = kernel_taylor(self._direct, mp.mpf(1), self.nmax, radius=mp.mpf('0.4'), dps=self.kd + 12) # regions B, C1, C2: barycentric Chebyshev (numerically stable) ln10 = mp.log(10) self.bary = [] for (a, b) in ((mp.mpf('1.38'), mp.mpf(9)), (mp.mpf(9), mp.mpf(40)), (mp.mpf(40), mp.mpf(200))): c, L = (a + b) / 2, (b - a) / 2 rho_e = c / L + mp.sqrt((c / L) ** 2 - 1) # sing at w=0 N = int(mp.ceil(qd * ln10 / mp.log(rho_e))) + 15 xs = [c + L * mp.cos(mp.pi * j / N) for j in range(N + 1)] fs = [self._direct(x) for x in xs] ws = [(mp.mpf(1) if j % 2 == 0 else mp.mpf(-1)) * (mp.mpf('0.5') if j in (0, N) else mp.mpf(1)) for j in range(N + 1)] self.bary.append((a, b, xs, fs, ws)) mp.mp.dps = old if verbose: self.verify() def _direct(self, w, dps=None): return box1(self.s, self.t, self.m2, w, dps=dps or self.kd) # 2026-09-06 (Q20f): the far-tail floor. The closed dilog box's four-root # partial fraction resolves mal = (p1 - sd1)/(2 M5) ~ 1/M5 only while its working digits # (dps + 18) exceed log10 w: past that the discriminant rounds to p1^2, # sd1 == p1, mal = mbe = 0, log 0 = -inf and Ral*log(mal) + Rbe*log(mbe) # is inf - inf = NaN (measured onset w ~ 1e37.5 at the served floor of 18, # dps-independent; reached by every --dps >= 72, whose W_MAX = 10^(dps//2+2) # exceeds it). Below TAIL_W_SWITCH the served relaxation is kept EXACTLY # (every documented tier <= dps 55 has W_MAX <= 1e29 there, so no such tier # can change a bit); at or above it the box dps is floored at # ceil(log10 w) + TAIL_MARGIN, which resolves mal to TAIL_MARGIN + 18 = 30 # relative digits (the artifact's value-claim bar) at every tail node, and # every tail value is asserted finite (rc 1 naming the node and the # quantity: fail closed). Below the onset the served relaxation still caps # the kernel's tail accuracy (see the module docstring); this floor does # not lift that cap. TAIL_W_SWITCH = 10 ** 36 # a half-decade under the measured onset 1e37.5 TAIL_MARGIN = 12 def _tail(self, w): """direct eval with dps relaxed by the tail suppression (w/200)^-2; at or above TAIL_W_SWITCH the dps is floored at ceil(log10 w) + TAIL_MARGIN and the value asserted finite (2026-09-06, Q20f).""" key = mp.nstr(w, 30) v = self._direct_cache.get(key) if v is None: drop = int(2 * mp.log10(w / 200)) d = max(18, self.kd - drop) if w >= self.TAIL_W_SWITCH: d = max(d, int(mp.ceil(mp.log10(w))) + self.TAIL_MARGIN) v = self._direct(w, dps=d) if not mp.isfinite(v): raise RuntimeError( f"[kernel] K(w) (the closed dilog box on the far tail) is " f"not finite (nan/inf) at the tail node w = {mp.nstr(w, 8)} " f"(box dps {d}) -- fail closed") self._direct_cache[key] = v return v def __call__(self, w): w = mp.mpf(w) if mp.mpf('0.62') <= w <= mp.mpf('1.38'): u = w - 1 acc = mp.mpf(0) for cn in reversed(self.tay1): acc = acc * u + cn return acc for (a, b, xs, fs, ws) in self.bary: if a <= w <= b: num = den = mp.mpf(0) for xj, fj, wj in zip(xs, fs, ws): d = w - xj if d == 0: return fj r = wj / d num += r * fj den += r return num / den return self._tail(w) def verify(self): """live self-check: interpolants vs the closed form at interior points.""" worst = mp.inf for w in ('1.03', '1.21', '2.71', '6.9', '13.7', '33.3', '77.7', '181'): worst = min(worst, agree_digits(self(mp.mpf(w)), self._direct(mp.mpf(w)))) print(f" kernel interpolant self-check (8 pts): " f">= {mp.nstr(worst, 5)} d vs direct closed form") if worst < QUAD_DPS + 3: raise SystemExit("kernel interpolant below target accuracy -- STOP") # -------- the exact partial-fraction stub kernels K_P, K_1 (regular at w=1) -- class KernelPF: """K_P (npow=2) or K_1 (npow=1): [K(w) - K0 - (npow-1)*K1*(w-1)]/(w-1)^npow. Taylor route (coefficients = K.tay1 shifted by npow) inside |w-1|= " f"{mp.nstr(worst, 5)} d") if worst < QUAD_DPS - 2: raise SystemExit(f"{label} seam check below target -- STOP") return worst def box_dot_fd(K, dps): """K'(1) by Richardson-extrapolated central finite differences on the closed form -- the route INDEPENDENT of the Chebyshev fit.""" old = mp.mp.dps mp.mp.dps = dps one = mp.mpf(1) hs = [mp.mpf(10) ** (-e) for e in (8, 10, 12, 14, 16, 18, 20)] raw = [(K._direct(one + h, dps=dps) - K._direct(one - h, dps=dps)) / (2 * h) for h in hs] r4 = mp.mpf('1e-4') # (h_{k+1}/h_k)^2 R1 = [(raw[k + 1] - r4 * raw[k]) / (1 - r4) for k in range(len(raw) - 1)] mp.mp.dps = old return R1[-2] # ============ B. subtracted tanh-sinh dispersion quadrature ================== # Port of tools/dispersion_subtracted/disp_sub.py (Dispersify): subtract the # closed-form threshold model S(w) on a sub-panel, add back \int S*K in closed # form (kernel Taylor x exact power-log moments), tanh-sinh the analytic # remainder. Exponentially convergent in the level L. def moment_powerlog(p, m, L): r"""\int_0^L u^p (log u)^m du, closed form (p>-1, m in {0,1,2}).""" p = mp.mpf(p); L = mp.mpf(L) pp = p + 1 Lp = L ** pp lnL = mp.log(L) if m == 0: return Lp / pp if m == 1: return Lp / pp * (lnL - 1 / pp) if m == 2: return Lp / pp * (lnL * lnL - 2 * lnL / pp + 2 / pp ** 2) raise ValueError("m in {0,1,2}") def kernel_taylor(K, wstar, nmax, radius=None, dps=40): """Taylor coeffs of the analytic kernel about w* via Chebyshev collocation (real evaluations only), computed at runtime.""" wstar = mp.mpf(wstar) radius = mp.mpf(radius if radius is not None else '0.4') old = mp.mp.dps mp.mp.dps = dps + 25 N = nmax + 1 us = [radius * mp.cos(mp.pi * (j + mp.mpf('0.5')) / N) for j in range(N)] fs = [mp.mpf(K(wstar + u)) for u in us] V = mp.matrix(N, N) for j in range(N): p = mp.mpf(1) for n in range(N): V[j, n] = p p *= us[j] a = mp.lu_solve(V, mp.matrix(fs)) mp.mp.dps = old return [a[n] for n in range(N)] def addback_endpoint(coeffs, Kc, L, side='left'): r"""Closed-form \int_panel S(w) K(w) dw from exact power-log moments.""" L = mp.mpf(L) tot = mp.mpf(0) sgn = mp.mpf(1) if side == 'left' else mp.mpf(-1) for (c, alpha, m) in coeffs: c = mp.mpf(c); alpha = mp.mpf(alpha) for n, Kn in enumerate(Kc): tot += c * Kn * (sgn ** n) * moment_powerlog(alpha + n, m, L) return tot def _Smodel(coeffs, wstar, side, wp): s = (wp - wstar) if side == 'left' else (wstar - wp) if s <= 0: return mp.mpf(0) val = mp.mpf(0) for (c, alpha, m) in coeffs: t = mp.mpf(c) * s ** mp.mpf(alpha) if m: t *= mp.log(s) ** int(m) val += t return val def disp_subtracted(rho, K, panels, dps=40, kernel_nmax=None, maxdegree=None, Kc_pre=None): r"""(1/pi) \int rho K dw -- the quadrature RUNS here at every call. Kc_pre: precomputed kernel-Taylor coeffs about the threshold (mandatory here: each stub kernel passes its own exact tay1 slice).""" work = dps + 20 old = mp.mp.dps mp.mp.dps = work if kernel_nmax is None: kernel_nmax = int(1.7 * dps) + 10 def _ts(f, a, b): if maxdegree is None: return mp.quad(f, [a, b], method='tanh-sinh') return mp.quad(f, [a, b], method='tanh-sinh', maxdegree=maxdegree) tot = mp.mpf(0) for P in panels: a = mp.mpf(P['a']) thr = P.get('thresh', None) if P.get('map') == 'tail': scale = mp.mpf(P.get('scale', a if a > 0 else mp.mpf(1))) def ft(t, a=a, scale=scale): wp = a + scale * (1 + t) / (1 - t) r = rho(wp) if r == 0: return mp.mpf(0) # beyond W_MAX guard return r * K(wp) * scale * 2 / (1 - t) ** 2 tot += _ts(ft, mp.mpf(-1), mp.mpf(1)) continue b = mp.mpf(P['b']) if thr is None: tot += _ts(lambda wp: rho(wp) * K(wp), a, b) continue side, wstar, coeffs = thr[0], mp.mpf(thr[1]), thr[2] delta = mp.mpf(P.get('sub_width', min(mp.mpf('0.5'), (b - a) / 2))) if side == 'left': sa, sb = a, a + delta ra, rb = a + delta, b else: sa, sb = b - delta, b ra, rb = a, b - delta gsm = lambda wp: (rho(wp) - _Smodel(coeffs, wstar, side, wp)) * K(wp) tot += _ts(gsm, sa, sb) Kc = Kc_pre if Kc_pre is not None else \ kernel_taylor(K, wstar, kernel_nmax, radius=P.get('radius'), dps=work) tot += addback_endpoint(coeffs, Kc[:kernel_nmax + 1], delta, side=side) if rb > ra: tot += _ts(lambda wp: rho(wp) * K(wp), ra, rb) res = tot / mp.pi mp.mp.dps = old return +res def lbl3kp_panels(A1, Wcut=200): """w=1 turn-on subtraction (A1=-2pi closed form); break at the w=9 Gamma_1(6) cut (finite jump); real [Wcut,inf) tail map. Same panels for all three kernels -> the memoized rho is shared across the ladders.""" return [{'a': 1, 'b': 9, 'thresh': ('left', 1, [(A1, mp.mpf(1), 1)]), 'sub_width': mp.mpf(0.4), 'radius': mp.mpf('0.4')}, {'a': 9, 'b': Wcut, 'thresh': None}, {'a': Wcut, 'map': 'tail'}] # ====== C. rho(w) as a runtime Picard-Fuchs series solution ================== # The exact rational 9x9 kite system (master 4 = [1,1,3,0,0] is a decoupled # sink, dropped -> 8x8) is eps-expanded at d=4-2eps to eps^{0,1,2} and the # 24-component Laurent-block system is solved by adaptive local Taylor series # from the single w=5 seed. ONE monotone sweep per zone covers (1,inf); every # step's local series is kept, so rho at ANY w is a Horner evaluation of an # analytic series -- no stored numeric nodes, precision set by (DPS, NORD, SF). _KEEP = [0, 1, 2, 3, 5, 6, 7, 8] _TOP = 7 # local index of J[1,1,1,1,1] in the kept list def load_kite_de(): """eps-expand the exact rational DE: the sibling kite_de_exact_parse.py parses the shipped symbolic system and grades it at d = 4 - 2 eps in exact rational arithmetic (fractions.Fraction; the same integer coefficient lists the earlier sympy parse produced, in lowest terms; no numerics happen there); a string outside the grammar is REFUSED by name (exit 3).""" J = json.load(open(os.path.join(HERE, 'lbl3kp-kite-de.json'))) n = len(_KEEP) funcs = [[[None] * n for _ in range(n)] for _ in range(3)] try: lists = kx.graded_lists(J['A'], _KEEP, 2) except kx.ParseError as ex: sys.stderr.write(f"{SCRIPT} REFUSED: lbl3kp-kite-de.json: {ex}\n") raise SystemExit(3) for (i, j, k), (num, den) in lists.items(): ii, jj = _KEEP.index(i), _KEEP.index(j) # descending coefficient lists, each integer converted as mpf(p) / mpf(q) nc = [mp.mpf(c.numerator) / mp.mpf(c.denominator) for c in reversed(num)] dc = [mp.mpf(c.numerator) / mp.mpf(c.denominator) for c in reversed(den)] funcs[k][ii][jj] = (nc, dc) return funcs, [tuple(J['masters'][i]) for i in _KEEP] _ARB_RE = None def _pb(sv): global _ARB_RE if _ARB_RE is None: import re _ARB_RE = re.compile(r'^\s*\[\s*([^\s\]]*)\s*\+/-') m = _ARB_RE.match(str(sv).strip()) return mp.mpf(m.group(1)) if m else mp.mpf(str(sv).strip()) def _seed_cache_check(st, rg): """Stored-vs-regenerated seed string check (comparator logic: /small_feeds/compare_seed_strings.py, w5 mode): PASS per component iff byte-prefix either way, OR final-digit round-to-nearest, OR numerical zero (each value below its own run's noise floor 10^-(dps+5) -- the mathematically-vanishing components: sub-threshold Im parts and the finite kite's eps^-2/eps^-1).""" assert st['masters'] == rg['masters'] and st['orders'] == rg['orders'] ok = ncomp = 0 bad = [] with mp.workdps(max(int(st['dps']), int(rg['dps'])) + 20): zt_st = mp.mpf(10) ** (-(int(st['dps']) + 5)) zt_rg = mp.mpf(10) ** (-(int(rg['dps']) + 5)) for mi, mname in enumerate(st['masters']): for oi, o in enumerate(st['orders']): for c, part in enumerate(('re', 'im')): s_st = st['values'][mi][oi][c] s_new = rg['values'][mi][oi][c] ncomp += 1 short, long_ = sorted((s_st, s_new), key=len) if long_.startswith(short): ok += 1 continue p = 0 for x, y in zip(s_st, s_new): if x != y: break p += 1 if p >= min(len(s_st), len(s_new)) - 1: # last-digit rounding ok += 1 continue if abs(mp.mpf(s_st)) < zt_st and abs(mp.mpf(s_new)) < zt_rg: ok += 1 # numerical zeros continue bad.append((str(mname), o, part, f"prefix {p}/{min(len(s_st), len(s_new))}")) return ok, ncomp, bad def load_kite_boundary(masters): """Laurent coeffs eps^{-2..0} of the 8 masters at w=5 -- the transport SEED, now the DERIVED AMFlow-free vector lbl3kp-w5-derived.json (emitted once at dps 80 by the sibling lbl3se-w5-seed.py, rerunnable from this directory via the lbl3kp-w5-seed.py shim: p^2=0 vacuum closed forms + Frobenius at w=0 + exact-DE Taylor march to w=5). With --boundary-recompute [DPS] the seed is REGENERATED live through the shim and the stored strings demote to a byte-agreement cache check. The retired AMFlow w=5 vector (lbl3kp-kite-boundary-w5.json) is NOT consumed upstream: it is loaded only as a held-out cross-check of the derived seed, and the recomputed agreement is returned for printing.""" D = json.load(open(os.path.join(HERE, 'lbl3kp-w5-derived.json'))) if getattr(_ARGS, 'boundary_recompute', None): import subprocess import sys as _sys import tempfile _bdps = _ARGS.boundary_recompute print(f" [boundary-recompute] deriving the w=5 seed live: " f"lbl3kp-w5-seed.py --dps {_bdps} --json ...") _tf = tempfile.NamedTemporaryFile(suffix='.json', delete=False) _tf.close() subprocess.run([_sys.executable, os.path.join(HERE, 'lbl3kp-w5-seed.py'), '--dps', str(_bdps), '--json', _tf.name, '-q'], check=True) D_LIVE = json.load(open(_tf.name)) os.unlink(_tf.name) _ok, _ncomp, _bad = _seed_cache_check(D, D_LIVE) print(f" [boundary-recompute] stored-string cache check: " f"{_ok}/{_ncomp} -> {'PASS' if _ok == _ncomp else 'FAIL'}") if _ok != _ncomp: raise SystemExit(f"[boundary-recompute] cache check FAILED: {_bad[:5]}") D = D_LIVE dmast = [tuple(m) for m in D['masters']] n = len(masters) out = {k: [mp.mpc(0)] * n for k in (-2, -1, 0)} for loc, idx in enumerate(masters): i = dmast.index(idx) for t, o in enumerate(D['orders']): if o in out: re_s, im_s = D['values'][i][t] out[o][loc] = mp.mpc(mp.mpf(re_s), mp.mpf(im_s)) # held-out cross-check: derived seed vs the retired AMFlow vector J = json.load(open(os.path.join(HERE, 'lbl3kp-kite-boundary-w5.json'))) worst = mp.inf for r in J['result']: idx = tuple(r['integral']['indices']) if idx not in masters: continue loc = masters.index(idx) for c in r['coefficients']: o = c['order'] if o not in out: continue ref = mp.mpc(_pb(c['value']['re']), _pb(c['value']['im'])) dv = mp.fabs(out[o][loc] - ref) if dv != 0: worst = min(worst, -mp.log10(dv / max(mp.fabs(ref), mp.mpf(1)))) return out, worst, int(D['dps']) def _shiftpoly(coefs, x0, zero): """ascending Taylor coeffs about x0 of the poly with DESCENDING coefs.""" p = [coefs[0] + zero] for c in coefs[1:]: pn = [zero] * (len(p) + 1) for m_, a in enumerate(p): pn[m_] += a * x0 pn[m_ + 1] += a pn[0] += c p = pn return p def _sparse_A_taylor(funcs, w0, N, real): """Per-eps-order Taylor series of A about w0; real arithmetic on the real axis (2x faster complex-times-real products in the state recursion).""" zero = mp.mpf(0) if real else mp.mpc(0) n = len(_KEEP) out = [] for jj in range(3): ser_m = [[] for _ in range(N + 1)] for i in range(n): for jc in range(n): cd = funcs[jj][i][jc] if cd is None: continue nc, dc = cd ph = _shiftpoly(nc, w0, zero) qh = _shiftpoly(dc, w0, zero) q0 = qh[0] r = [1 / q0] for m_ in range(1, N + 1): sacc = zero for l in range(1, min(m_, len(qh) - 1) + 1): sacc += qh[l] * r[m_ - l] r.append(-sacc / q0) for m_ in range(N + 1): sacc = zero for l in range(min(m_, len(ph) - 1) + 1): sacc += ph[l] * r[m_ - l] if sacc != 0: ser_m[m_].append((i, jc, sacc)) out.append(ser_m) return out def taylor_step(funcs, M, w0, h, N, real): """One local-series step of the Laurent-block system. Returns the state at w0+h AND the eps^0 TOP-component series (for later rho evaluation).""" n = len(_KEEP) Aser = _sparse_A_taylor(funcs, w0, N, real) C = {k: [list(M[k])] for k in (-2, -1, 0)} for m_ in range(N): for k in (-2, -1, 0): sacc = [mp.mpc(0)] * n for j in range(3): if k - j < -2: continue Clist = C[k - j] Alist = Aser[j] for l in range(m_ + 1): cv = Clist[m_ - l] for (ri, ci, a) in Alist[l]: sacc[ri] += a * cv[ci] inv = mp.mpf(1) / (m_ + 1) C[k].append([x * inv for x in sacc]) out = {} for k in (-2, -1, 0): v = [mp.mpc(0)] * n hp = mp.mpc(1) if not real else mp.mpf(1) for m_ in range(N + 1): cm = C[k][m_] for i in range(n): v[i] += cm[i] * hp hp *= h out[k] = v top_series = [C[0][m_][_TOP] for m_ in range(N + 1)] return out, top_series class RhoPF: """rho(w) = -Im J_top(w+i0) as a runtime PF-series solution on (1,inf).""" def __init__(self, funcs, M0, dps=DPS, N=NORD, sf=SF, verbose=True): self.funcs, self.N, self.sf = funcs, N, sf self.dps = dps self.segs = [] # (a, b, w0, top_eps0_series) self.checkpoints = {} # w -> full state (for closed-form master checks) self._starts = None self._memo = {} # runtime memo shared by the three ladders old = mp.mp.dps mp.mp.dps = dps t0 = time.time() nsteps = 0 # zone M: 5 -> 9-V_MIN (up); snapshot the state near 8.4 for the arc hop _, snap = self._leg(M0, mp.mpf(5), 9 - V_MIN, snap_at=mp.mpf('8.4')) nsteps += self._ns self.rho_9m = self._scan_eval(9 - V_MIN) # zone L: 5 -> 1+U_MIN (down) self._leg(M0, mp.mpf(5), 1 + U_MIN) nsteps += self._ns # arc hop: snapshot -> 9.5 via the Im(w)>0 semicircle (Feynman s+i0 # sheet; the Im<0 arc gives the conjugate branch -- campaign-verified) wsn, Marc = snap for wp in (mp.mpc('8.5', '0.5'), mp.mpc(9, '0.5'), mp.mpc('9.5', '0.5'), mp.mpc('9.5')): Marc, _ = self._leg(Marc, wsn, wp, store=False) wsn = wp nsteps += self._ns # zone R2: 9.5 -> 9+V_MIN (down) self._leg(Marc, mp.mpf('9.5'), 9 + V_MIN) nsteps += self._ns self.rho_9p = self._scan_eval(9 + V_MIN) # zone R3: 9.5 -> W_MAX (up); checkpoint the full state at w=12 self._leg(Marc, mp.mpf('9.5'), W_MAX, want=(mp.mpf(12),)) nsteps += self._ns self.segs.sort(key=lambda s: s[0]) self._starts = [s[0] for s in self.segs] self._wcap = max(s[1] for s in self.segs) self.nsteps = nsteps mp.mp.dps = old if verbose: print(f" PF sweep: {nsteps} local-series steps " f"(dps={dps}, order={N}, step={float(sf)}*dist), " f"{len(self.segs)} stored segments, {time.time()-t0:.1f}s") # -- internals ------------------------------------------------------------ def _leg(self, M, w_from, w_to, store=True, snap_at=None, want=()): """straight leg w_from -> w_to (real or complex), adaptive local series. Returns (final state, snapshot (w, state) if snap_at was crossed).""" self._ns = 0 M = {k: list(M[k]) for k in M} w = w_from snap = None want = list(want) real = (mp.im(mp.mpc(w)) == 0 and mp.im(mp.mpc(w_to)) == 0) if real: w, w_to = mp.re(mp.mpc(w)), mp.re(mp.mpc(w_to)) tol = mp.mpf(10) ** (-self.dps + 8) while abs(w_to - w) > tol * max(1, abs(w_to)): d = min(abs(w), abs(w - 1), abs(w - 9)) hmax = self.sf * d dirn = w_to - w h = dirn if abs(dirn) <= hmax else dirn / abs(dirn) * hmax if real: h = mp.re(mp.mpc(h)) M2, tser = taylor_step(self.funcs, M, w, h, self.N, real) if real: lo, hi = (w, w + h) if h > 0 else (w + h, w) if store: self.segs.append((lo, hi, w, tser)) if snap_at is not None and snap is None and lo <= snap_at <= hi: snap = (w + h, {k: list(M2[k]) for k in M2}) for wc in list(want): if lo <= wc <= hi: Mc, _ = taylor_step(self.funcs, M, w, wc - w, self.N, real) self.checkpoints[wc] = Mc want.remove(wc) M = M2 w = w + h self._ns += 1 return M, snap def _scan_eval(self, w): """rho during construction (linear scan of stored segments).""" for (a, b, w0, ser) in reversed(self.segs): if a <= w <= b: return self._horner(ser, w - w0) raise KeyError(f"rho: w={mp.nstr(w, 20)} not covered yet") def _eval(self, w): """rho at w from the covering stored local series (Horner).""" i = bisect.bisect_right(self._starts, w) - 1 for j in (i, i + 1, i - 1): if 0 <= j < len(self.segs): a, b, w0, ser = self.segs[j] if a <= w <= b: return self._horner(ser, w - w0) return None @staticmethod def _horner(ser, hh): acc = mp.mpc(0) for cn in reversed(ser): acc = acc * hh + cn return -mp.im(acc) # -- public --------------------------------------------------------------- def __call__(self, w): w = mp.mpf(w) key = mp.nstr(w, self.dps + 6) v = self._memo.get(key) if v is not None: return v v = self._eval_raw(w) self._memo[key] = v return v def _eval_raw(self, w): if w <= 1: return mp.mpf(0) u = w - 1 if u < U_MIN: return u * (-2 * mp.pi) * mp.log(u) # closed-form turn-on word if abs(w - 9) < V_MIN: return self.rho_9m if w < 9 else self.rho_9p if w > W_MAX: return mp.mpf(0) v = self._eval(w) if v is not None: return v # rounding-gap fallbacks at the guard edges (widths ~1 ulp of the guards) if u < 2 * U_MIN: return u * (-2 * mp.pi) * mp.log(u) if abs(w - 9) < 2 * V_MIN: return self.rho_9m if w < 9 else self.rho_9p if w > mp.mpf('0.9') * self._wcap: return mp.mpf(0) raise KeyError(f"rho: w={mp.nstr(w, 20)} not covered") # ====== C2. rho(w) in CLOSED FORM: the printed density (2026-09-11) ========== # The kite spectral density written out in closed form (the write-up's eqs. # rho_elem / rho_full / F23; m^2 = 1, w = l^2/m^2): # # rho(w) = rho_elem(w) 1 < w <= 9 # rho(w) = rho_elem(w) - R(w)/w, R(w) = int_9^w v F23(v) dv w > 9 # rho_elem(w) = -(pi/w) [ 2 ln(w-1) ln w + 3 Li2(1-w) ] # F23(v) = [ 2 Im S(v) - (v+3) Im Sdot(v) ] / ( v (v-1)^2 ), # x_pm = (sqrt v pm 1)^2 # Im S(v) = -(pi/v) int_4^{x_-} sqrt((x_- - s)(x_+ - s)) sqrt((s-4)/s) ds # Im Sdot(v) = +(pi/v) int_4^{x_-} (v-1+s)/sqrt((x_- - s)(x_+ - s)) # * sqrt((s-4)/s) ds # # rho_elem is the weight-two part (the bubble and (m,0,0)-sunrise pieces of # the kite system's inhomogeneity integrated in closed form; its threshold # expansion -2 pi (w-1) ln(w-1) + 3 pi (w-1) carries the turn-on constant # A1 = -2 pi as a theorem of the closed form); Im S is the d = 4 equal-mass # sunrise three-body cut (the Gamma_1(6) object) and Im Sdot the dotted- # sunrise cut, so above w = 9 the density is a one-fold over the sunrise cut. # THIS is the density the script evaluates by default at the density gate # points (w = 5, 12, 100) in every tier, gated against the AMFlow # solve_integrals references at the served floors; the Picard-Fuchs series # density (RhoPF above: the transport of the exact kite system, which the # dispersion quadrature consumes -- the value path of the integrals is not # touched by this section) is the cross-check, the agreement of the two # densities printed and gated at the rho self-check bar. --density series # restores the previous prints (the series density alone, the closed form # not evaluated). --sum-rule runs the spot gate (1/pi) int_1^W rho dw = # 6 zeta_3 (the large-l^2 normalisation of the kite) on the closed form. # NUMERICS (every difference formed without cancellation): the two # s-integrals are taken in two pieces, s = 4 cosh^2 xi on [4, s_mid] # (sqrt((s-4)/s) = tanh xi) and, on [s_mid, x_-] with s_mid = 2 sqrt(x_-), # s = s_mid + D cos^2 phi, sin^2 phi = (4 sqrt v / D) sinh^2 eta (so that # x_+ - s = 4 sqrt v cosh^2 eta), D = x_- - s_mid -- tanh-sinh on the raw # limits loses half the working digits in Im Sdot at the 1/sqrt(x_- - s) # endpoint; R(w) by Gauss-Legendre panels in y = ln v with 9, 12, 100 and the # cutoff among the panel ends (the node count per panel from the Bernstein # ellipse to the nearest singularity y = 0, as for the kernel device); the # inner one-folds at CLOSED_DPS_IN = DPS + 11 digits, the sums at DPS + 20. CLOSED_DPS_IN = DPS + 11 # inner one-fold working digits (the closed form's accuracy) CLOSED_DPS = DPS + 20 # outer sums / Li2 / the sum-rule integrals CLOSED_AMBIENT = max(90, CLOSED_DPS + 10) # parse / compare precision of the gates CLOSED_PANEL_ENDS = ('9', '9.25', '9.5', '10', '11', '12', '13', '16', '20', '30', '50', '100', '300', '1000', '1e4', '1e6', '1e9', '1e14', '1e22', '1e34', '1e50', '1e75', '1e110', '1e160', '1e240') CLOSED_ELEM_ENDS = ('1.5', '2', '3', '5', '9', '16', '30', '100', '1000', '1e5', '1e8', '1e12', '1e18', '1e26', '1e38', '1e57', '1e85', '1e128', '1e190') SUMRULE_LOG10_W_RANGE = (34, 240) # --sum-rule LOG10_W static range (default 50); within it the # next-omitted-tail rule of _sum_rule_tier refuses a cutoff # too low for the running --dps by name def rho_elem(w): """The weight-two closed form -(pi/w)[2 ln(w-1) ln w + 3 Li2(1-w)] at the ambient precision; 0 at and below the threshold w = 1.""" w = mp.mpf(w) if w <= 1: return mp.mpf(0) return -(mp.pi / w) * (2 * mp.log(w - 1) * mp.log(w) + 3 * mp.polylog(2, 1 - w)) def _sunrise_cut(v, dps, dotted): """Im S(v) (dotted=False) or Im Sdot(v) (dotted=True) for v > 9 from the printed one-folds, in two pieces with the cancellation-free substitutions of the block comment. Returns (value, tanh-sinh engine error estimate).""" with mp.workdps(dps): v = mp.mpf(v) sv = mp.sqrt(v) xm = (sv - 1) ** 2 # x_- g = 4 * sv # x_+ - x_- smid = 2 * mp.sqrt(xm) D = xm - smid ximax = mp.acosh(mp.sqrt(smid) / 2) def f1(xi): # s = 4 cosh^2 xi on [4, s_mid] ch, sh = mp.cosh(xi), mp.sinh(xi) s = 4 * ch * ch xms = xm - s if xms <= 0: return mp.mpf(0) xps = g + xms jac = mp.tanh(xi) * 8 * ch * sh if dotted: return (v - 1 + s) / mp.sqrt(xms * xps) * jac return mp.sqrt(xms * xps) * jac I1, e1 = mp.quad(f1, [0, ximax], error=True) etamax = mp.asinh(mp.sqrt(D / g)) def f2(eta): # s = s_mid + D cos^2 phi on [s_mid, x_-] sh, ch = mp.sinh(eta), mp.cosh(eta) S2 = (g / D) * sh * sh # sin^2 phi if S2 >= 1: return mp.mpf(0) s = smid + D * (1 - S2) r = mp.sqrt(((smid - 4) + D * (1 - S2)) / s) if dotted: return 2 * (v - 1 + s) * r return 2 * g * g * r * ch * ch * sh * sh I2, e2 = mp.quad(f2, [0, etamax], error=True) val = (mp.pi / v) * (I1 + I2) return (val if dotted else -val), abs(e1) + abs(e2) def F23(v, dps): """The sunrise-cut inhomogeneity F23(v) of the printed density (v > 9) at dps inner digits. Returns (value, engine error estimate).""" with mp.workdps(dps): v = mp.mpf(v) iS, eS = _sunrise_cut(v, dps, False) iSd, eSd = _sunrise_cut(v, dps, True) den = v * (v - 1) ** 2 return (2 * iS - (v + 3) * iSd) / den, (2 * eS + (v + 3) * eSd) / den class RhoClosed: """The closed-form density: rho_elem below w = 9, rho_elem - R(w)/w above, R(w) = int_9^w v F23(v) dv on Gauss-Legendre panels in y = ln v (F23 at dps_in inner digits at every node; the panel integrals of v F23 and, for the sum rule, of v F23 ln(W/v) from the same nodes). Panels are built on demand up to the requested w; a w strictly inside a panel gets one more rule on the part panel [panel start, w].""" def __init__(self, dps_in=None, dps_out=None, wmax=None): self.dps_in = dps_in or CLOSED_DPS_IN self.dps = dps_out or CLOSED_DPS ends = [mp.mpf(e) for e in CLOSED_PANEL_ENDS] if wmax is not None: wmax = mp.mpf(wmax) ends = [e for e in ends if e < wmax] + [wmax] self.ends = ends self.panels = [] # per panel: dict(a, b, n, nodes [(y, v, f, wt*h*v^2)], R) self.Rcum = {0: mp.mpf(0)} # index k -> R(ends[k]) self.n_eval = 0 self.err_in = mp.mpf(0) # engine estimates of the inner one-folds, propagated into R self.wall = 0.0 @staticmethod def n_nodes(ya, yb, dps_in): """Gauss-Legendre node count for the panel [ya, yb] in y = ln v: the Bernstein-ellipse rate to the nearest singularity of F23 in y (y = 0, i.e. v = 1) for dps_in + 4 digits, floor 12.""" c, L = (ya + yb) / 2, (yb - ya) / 2 r = c / L + mp.sqrt((c / L) ** 2 - 1) return max(12, int(mp.ceil((dps_in + 4) * mp.log(10) / (2 * mp.log(r)))) + 4) def _rule(self, a, b): """one Gauss-Legendre rule on [a, b] in y = ln v: returns (n, nodes, sum of v F23 dv).""" ya, yb = mp.log(a), mp.log(b) h, c = (yb - ya) / 2, (ya + yb) / 2 n = self.n_nodes(ya, yb, self.dps_in) xs, ws = _gl(n, self.dps) nodes = [] acc = mp.mpf(0) for x, wt in zip(xs, ws): y = c + h * x v = mp.exp(y) f, e = F23(v, self.dps_in) self.n_eval += 1 gw = wt * h * v * v # int v F23 dv = int v^2 F23 dy acc += gw * f self.err_in += abs(gw) * e nodes.append((y, v, f, gw)) return n, nodes, acc def _build_panel(self, k): """panel k = [ends[k], ends[k+1]]: F23 at the nodes, the R increment.""" t0 = time.time() with mp.workdps(self.dps): a, b = self.ends[k], self.ends[k + 1] n, nodes, pR = self._rule(a, b) self.panels.append({'a': a, 'b': b, 'n': n, 'nodes': nodes, 'R': pR}) self.Rcum[k + 1] = self.Rcum[k] + pR self.wall += time.time() - t0 def _extend(self, w): """build every panel whose end is <= w; returns the index j of the largest panel end ends[j] <= w (0 when w is below the first end).""" k = len(self.panels) while k + 1 < len(self.ends) and self.ends[k + 1] <= w: self._build_panel(k) k += 1 j = k while j > 0 and self.ends[j] > w: j -= 1 return j def R(self, w): """R(w) = int_9^w v F23(v) dv (0 at and below 9).""" w = mp.mpf(w) if w <= 9: return mp.mpf(0) if w > self.ends[-1]: raise ValueError(f"RhoClosed: w = {mp.nstr(w, 8)} lies beyond the last panel end " f"{mp.nstr(self.ends[-1], 6)}") j = self._extend(w) acc = self.Rcum[j] if self.ends[j] < w: # a part panel [ends[j], w]: one more rule t0 = time.time() with mp.workdps(self.dps): acc = acc + self._rule(self.ends[j], w)[2] self.wall += time.time() - t0 return acc def __call__(self, w): w = mp.mpf(w) with mp.workdps(self.dps): if w <= 9: return rho_elem(w) return rho_elem(w) - self.R(w) / w # -- the sum rule (1/pi) int_1^W rho dw = 6 zeta_3 ------------------------- def sum_rule(self, W): """The finite-cutoff Fubini form (exact): int_1^W rho dw = int_1^W rho_elem dw - int_9^W v F23(v) ln(W/v) dv, W a panel end (the constructor's wmax). The elementary integral: tanh- sinh on [1, 3/2], then Gauss-Legendre panels in ln w. Returns a dict of mpf: I_elem, J (the sunrise-cut part), S = (I_elem - J)/pi, six_zeta3, tail = the leading omitted tail (4 pi (ln W + 1) + 6 pi)/(pi W) of (1/pi) int_W^inf rho, from rho ~ (4 pi ln w + 6 pi)/w^2 at large w.""" W = mp.mpf(W) j = self._extend(W) if self.ends[j] != W: raise ValueError("RhoClosed.sum_rule: W must be a panel end (construct with wmax=W)") t0 = time.time() with mp.workdps(self.dps): Y = mp.log(W) J = mp.mpf(0) for P in self.panels[:j]: for (y, v, f, gw) in P['nodes']: J += gw * f * (Y - y) el_ends = [mp.mpf(e) for e in CLOSED_ELEM_ENDS] el_ends = [e for e in el_ends if e < W] + [W] I_el = mp.quad(rho_elem, [1, el_ends[0]]) n_el = 0 for a, b in zip(el_ends[:-1], el_ends[1:]): ya, yb = mp.log(a), mp.log(b) h, c = (yb - ya) / 2, (ya + yb) / 2 n = self.n_nodes(ya, yb, self.dps_in) xs, ws = _gl(n, self.dps) for x, wt in zip(xs, ws): wv = mp.exp(c + h * x) I_el += wt * h * wv * rho_elem(wv) n_el += 1 S = (I_el - J) / mp.pi tail = (4 * mp.pi * (Y + 1) + 6 * mp.pi) / (mp.pi * W) res = {'W': W, 'I_elem': I_el, 'J': J, 'S': S, 'six_zeta3': 6 * mp.zeta(3), 'tail': tail, 'n_elem_nodes': n_el, 'n_elem_panels': len(el_ends) - 1, 'wall_elem': time.time() - t0} return res def _closed_fail(msg): raise RuntimeError(f"[closed form] DENSITY GATE FAIL: {msg} -- the closed-form density does not " "reproduce its references; STOP") def closed_density_gates(series, series_label, bar12, bar_series, indent=' '): """The closed-form density evaluated NOW at the three density gate points w = 5, 12, 100 and gated: against the AMFlow solve_integrals references (RHO5_STR / RHO_W12_STR / RHO100_STR) at the served floors (GATE_FLOORS at the --dps target for w = 5 and 100; bar12, the rho self-check bar of the calling tier, for w = 12), and against the series density `series` ({'rho(5)': mpf, 'rho(12)': mpf, 'rho(100)': mpf}: the cached strings in the fast-start tier, the live transport's values in the live tiers) at bar_series. Prints one value line (min(50, CLOSED_DPS_IN - 10) digits) and one gate line per point; RAISES (rc 1) naming every failed comparison. Returns the RhoClosed object.""" t0 = time.time() pad = ' ' * len(indent) print(f"{indent}[closed form] the density in CLOSED FORM (rho_elem = -(pi/w)[2 ln(w-1) ln w + 3 Li2(1-w)] " "below w = 9;") print(f"{pad}rho_elem - R(w)/w above, R the one-fold over the equal-mass sunrise cut) evaluated NOW at the") print(f"{pad}density gate points; the cross-check is the series density ({series_label}):") rc = RhoClosed(wmax=100) vals = {name: rc(wv) for name, wv in (('rho(5)', 5), ('rho(12)', 12), ('rho(100)', 100))} refs = {'rho(5)': (RHO5_STR, 'AMFlow solve_integrals w=5', float(GATE_FLOORS['rho(5)'](QUAD_DPS))), 'rho(12)': (RHO_W12_STR, 'AMFlow solve_integrals w=12', float(bar12)), 'rho(100)': (RHO100_STR, 'AMFlow solve_integrals bnd_w100', float(GATE_FLOORS['rho(100)'](QUAD_DPS)))} failed = [] with mp.workdps(CLOSED_AMBIENT): for name in ('rho(5)', 'rho(12)', 'rho(100)'): val = vals[name] lit, what, floor = refs[name] d_ref = agree_digits(val, mp.mpf(lit)) d_ser = agree_digits(val, mp.mpf(series[name])) if name in series else None ok_ref = bool(d_ref >= floor) ok_ser = True if d_ser is None else bool(d_ser >= bar_series) print(f"{pad}[closed form] {name:8s} = {mp.nstr(val, min(50, CLOSED_DPS_IN - 10))}") line = (f"{pad} vs {what}: {mp.nstr(d_ref, 6)} d >= floor {floor:.2f} d -- " f"{'OK' if ok_ref else 'FAIL'}") if d_ser is not None: line += (f"; vs the series density: {mp.nstr(d_ser, 6)} d >= bar {float(bar_series):.0f} d -- " f"{'OK' if ok_ser else 'FAIL'}") print(line) if not ok_ref: failed.append(f"{name} vs {what} {mp.nstr(d_ref, 6)} d < floor {floor:.2f} d") if not ok_ser: failed.append(f"{name} vs the series density {mp.nstr(d_ser, 6)} d < bar {float(bar_series):.0f} d") err = mp.nstr(rc.err_in, 3) n_ser = sum(1 for k in vals if k in series) print(f"{pad}[closed form] density gate: {3 - sum(1 for f in failed if 'AMFlow' in f)}/3 points within the " f"AMFlow floors, {n_ser - sum(1 for f in failed if 'series' in f)}/{n_ser} within the series bar " f"(dps={QUAD_DPS}) -- RAISES on fail") print(f"{pad} ({rc.n_eval} F23 evaluations on {len(rc.panels)} panels over [9, 100], inner one-folds at " f"{rc.dps_in} digits, sums at {rc.dps}; the engine's estimate for the one-folds, propagated into R: {err}, " f"informational; {time.time()-t0:.1f}s)") if failed: _closed_fail('; '.join(failed)) return rc def _sum_rule_tier(): """--sum-rule [LOG10_W]: the spot gate (1/pi) int_1^inf rho(w) dw = 6 zeta_3 on the closed-form density, at the finite cutoff W = 10^LOG10_W in the Fubini form of RhoClosed.sum_rule (the elementary and the sunrise-cut parts integrated separately, so their large-w cancellation never happens numerically). GATE: the no-tail residual against the cutoff's own truncation CAPPED at the inner working precision (bar min(floor(-log10( tail / 6 zeta_3)) - 2, CLOSED_DPS_IN - 12), tail = the leading omitted tail (4 pi (ln W + 1) + 6 pi)/(pi W); a residual is never resolved below the inner one-folds' working precision, so the truncation bar alone -- 66 at W = 1e70, 95 at W = 1e100 against 66 inner digits at the default dps -- would fail correct bytes), the tail-restored residual against the inner working precision (bar CLOSED_DPS_IN - 12), and the density gate points w = 5, 12, 100 from the same grid against the AMFlow references at the served floors. ACCEPTED CUTOFFS: LOG10_W in SUMRULE_LOG10_W_RANGE and such that the NEXT omitted tail clears the tail-restored bar: that tail is below (4 ln W + 8)/W^2 (rho's next large-w term measured below 2 (4 pi ln w + 6 pi)/w^3: rho w^2/(4 pi ln w + 6 pi) - 1 = 1.73/w, 1.87/w, 1.94/w at w = 1e4, 1e9, 1e22), and the tier requires floor(-log10((4 ln W + 8)/(W^2 6 zeta_3))) - 2 >= CLOSED_DPS_IN - 12 -- the whole static range at the default dps; a cutoff too low for a raised --dps is refused by name (rc 2) with the smallest LOG10_W that passes. rc 0 PASS / 1 FAIL / 2 usage or refused cutoff.""" tag = SCRIPT.split('-')[0].upper() others = [n for n in ('point', 'fastcache_bank', 'boundary_recompute', 'checkpoint', 'resume') if getattr(_ARGS, n, None) not in (None, False)] if others: print(f"{SCRIPT}: --sum-rule is a tier of its own; it does not combine with " f"{', '.join('--' + n.replace('_', '-') for n in others)}", file=sys.stderr) raise SystemExit(2) if _ARGS.density != 'closed': print(f"{SCRIPT}: --sum-rule integrates the closed-form density; --density series turns it off", file=sys.stderr) raise SystemExit(2) lw = int(_ARGS.sum_rule) lo, hi = SUMRULE_LOG10_W_RANGE if not lo <= lw <= hi: print(f"{SCRIPT}: --sum-rule LOG10_W must lie in [{lo}, {hi}] (got {lw}), the tier's static range ({hi}: the " f"last panel end; below {lo} the cutoff's own truncation (4 ln W + 10)/W of the sum leaves the no-tail " "check under 32 digits); within the range a cutoff whose next omitted tail does not clear the " "tail-restored bar at the running --dps is refused by name", file=sys.stderr) raise SystemExit(2) W = mp.mpf(10) ** lw # the bars under the running precision: bar1 = the inner one-folds' working precision (the tail-restored bar and # the cap of the no-tail bar); bar2 = the agreement the NEXT omitted tail lets the tail-restored residual reach at # this cutoff, from the bound (4 ln W + 8)/W^2 on (1/pi) int_W^inf of rho's next large-w term (measured below # 2 (4 pi ln w + 6 pi)/w^3: rho w^2/(4 pi ln w + 6 pi) - 1 = 1.73/w, 1.87/w, 1.94/w at w = 1e4, 1e9, 1e22); a cutoff # with bar2 < bar1 cannot pass on correct bytes and is refused by name, with the smallest cutoff that can bar1 = CLOSED_DPS_IN - 12 with mp.workdps(30): z30 = 6 * mp.zeta(3) def next_tail(L): return (4 * L * mp.log(10) + 8) / mp.mpf(10) ** (2 * L) def next_bar(L): return int(mp.floor(-mp.log10(next_tail(L) / z30))) - 2 tail2, bar2 = next_tail(lw), next_bar(lw) if bar2 < bar1: lmin = next((L for L in range(lw + 1, hi + 1) if next_bar(L) >= bar1), None) print(f"{SCRIPT}: --sum-rule {lw} refused at --dps {QUAD_DPS}: the tail-restored residual is gated at " f"{bar1} d (CLOSED_DPS_IN - 12, the inner one-folds' working precision) and the next omitted tail at " f"W = 1e{lw}, below (4 ln W + 8)/W^2 = {mp.nstr(tail2, 3)}, lets it reach {bar2} d at most; " + (f"raise LOG10_W to {lmin} or more, or lower --dps" if lmin is not None else f"no cutoff up to LOG10_W = {hi} passes at this --dps; lower --dps"), file=sys.stderr) raise SystemExit(2) print(f"{tag} SUM-RULE tier -- the closed-form density's large-l^2 normalisation") print(" (1/pi) int_1^inf rho(w) dw = 6 zeta_3, checked at the finite cutoff W = 1e%d in the Fubini form" % lw) print(" int_1^W rho dw = int_1^W rho_elem dw - int_9^W v F23(v) ln(W/v) dv (exact; the elementary and the") print(" sunrise-cut parts integrated separately, so their large-w cancellation never happens numerically)") print(f" settings: QUAD_DPS={QUAD_DPS} DPS={DPS}; inner one-folds at CLOSED_DPS_IN={CLOSED_DPS_IN} digits, " f"sums at CLOSED_DPS={CLOSED_DPS};") print(f" bars: the tail-restored residual at CLOSED_DPS_IN - 12 = {bar1} d, the no-tail residual at the cutoff's own " f"truncation capped there; the next omitted tail at this cutoff, below (4 ln W + 8)/W^2 = {mp.nstr(tail2, 3)}, " f"allows {bar2} d;") print(" Gauss-Legendre panels in ln v / ln w, node counts from the Bernstein ellipse to ln = 0\n") print(f"[s1] [{time.time()-T0:6.1f}s] building the Gauss-Legendre panels of the sunrise-cut part on [9, 1e{lw}] " "and the elementary integral ...") rc = RhoClosed(wmax=W) res = rc.sum_rule(W) k = len(rc.panels) print(f" the sunrise-cut part: {k} panels ({', '.join(str(P['n']) for P in rc.panels)} nodes), " f"{rc.n_eval} F23 evaluations in {rc.wall:.1f}s;") print(f" the elementary part: tanh-sinh on [1, 3/2] + {res['n_elem_panels']} panels " f"({res['n_elem_nodes']} nodes) in {res['wall_elem']:.1f}s") with mp.workdps(CLOSED_AMBIENT): S, z, tail = res['S'], res['six_zeta3'], res['tail'] r0 = S - z # the residuals formed at the compare precision, printed below r1 = (S + tail) - z d0 = agree_digits(S, z) d1 = agree_digits(S + tail, z) bar_trunc = int(mp.floor(-mp.log10(tail / z))) - 2 # the cutoff's own truncation bar0 = min(bar_trunc, bar1) # capped at the inner working precision ok0, ok1 = bool(d0 >= bar0), bool(d1 >= bar1) nd = CLOSED_DPS_IN # the printed digits follow the inner precision print(f"[s2] [{time.time()-T0:6.1f}s] (1/pi) int_1^W rho dw = {mp.nstr(S, nd)}") print(f" 6 zeta_3 = {mp.nstr(z, nd)}") print(f" residual (no tail) = {mp.nstr(r0, 4)} [the omitted tail (1/pi) int_W^inf rho: " f"leading term (4 pi (ln W + 1) + 6 pi)/(pi W) = {mp.nstr(tail, 4)}") print(" from rho ~ (4 pi ln w + 6 pi)/w^2 (1 + c/w + ...), " f"c < 2 measured; the next term is below (4 ln W + 8)/W^2 = {mp.nstr(tail2, 3)}]") print(f" residual, tail added = {mp.nstr(r1, 4)}") print(f" [gate] no-tail agreement {mp.nstr(d0, 6)} d >= bar {bar0} d " f"(= min(floor(-log10(tail / 6 zeta_3)) - 2 = {bar_trunc}, CLOSED_DPS_IN - 12 = {bar1}): the cutoff's " f"own truncation, capped at the inner one-folds' working precision) -- {'OK' if ok0 else 'FAIL'}") print(f" [gate] tail-restored agreement {mp.nstr(d1, 6)} d >= bar {bar1} d " f"(= CLOSED_DPS_IN - 12: the inner one-folds' working precision; the next omitted tail allows {bar2} d " f"here) -- {'OK' if ok1 else 'FAIL'}") print(f"[s3] [{time.time()-T0:6.1f}s] the density gate points from the same grid (12 and 100 are panel ends; " "5 is rho_elem alone):") oks = [] for name, wv, lit, what, floor in ( ('rho(5)', 5, RHO5_STR, 'AMFlow solve_integrals w=5', float(GATE_FLOORS['rho(5)'](QUAD_DPS))), ('rho(12)', 12, RHO_W12_STR, 'AMFlow solve_integrals w=12', float(QUAD_DPS - 10)), ('rho(100)', 100, RHO100_STR, 'AMFlow solve_integrals bnd_w100', float(GATE_FLOORS['rho(100)'](QUAD_DPS)))): val = rc(wv) d = agree_digits(val, mp.mpf(lit)) ok = bool(d >= floor) oks.append(ok) print(f" {name:8s} = {mp.nstr(val, min(50, CLOSED_DPS_IN - 10))} vs {what}: {mp.nstr(d, 6)} d " f">= floor {floor:.2f} d -- {'OK' if ok else 'FAIL'}") ok_all = ok0 and ok1 and all(oks) lim = ("the cutoff's truncation" if bar_trunc <= bar1 else "the inner working precision, the cutoff's truncation below it") print(f"\ntotal wall time: {time.time()-T0:.1f}s") print(f"PASS (sum rule) -- (1/pi) int rho dw of the closed-form density reproduces 6 zeta_3 to " f"{int(mp.floor(d0))} digits at the cutoff W = 1e{lw} ({lim}) and to {int(mp.floor(d1))} " "digits with the leading tail restored; the density gate points within their floors" if ok_all else "FAIL (sum rule) -- see the [gate] / floor lines above") raise SystemExit(0 if ok_all else 1) # ====== D. closed-form checks of the non-elliptic kite masters =============== # AMFlow normalization: each loop carries measure d^dk/(i pi^{d/2}) e^{eps g}, # g = EulerGamma. With m^2 = 1: # TAD = -e^{eps g} Gamma(eps-1) (massive tadpole) # FMIX(w) = e^{eps g} Gamma(eps) \int_0^1 dx x^{-eps} (1-(1-x)w-i0)^{-eps} # so J[1,1,0,0,0] = TAD^2, J[1,1,0,1,0] = TAD*FMIX, J[1,1,0,1,1] = FMIX^2. # Checked LIVE against the transported Laurent blocks at the w=12 checkpoint # (small-eps probe; limited by the O(eps) truncation of the Laurent block). def _fmix(e, w): xstar = 1 - 1 / mp.mpf(w) def f_below(x): # 1-(1-x)w < 0: (A e^{-i pi})^{-eps} = A^{-eps} e^{i pi eps} A = (1 - x) * w - 1 return x ** (-e) * A ** (-e) * mp.exp(1j * mp.pi * e) def f_above(x): return x ** (-e) * (1 - (1 - x) * w) ** (-e) val = mp.quad(f_below, [0, xstar]) + mp.quad(f_above, [xstar, 1]) return mp.exp(e * mp.euler) * mp.gamma(e) * val def closed_masters_check(state, w, masters, eps_list=('1e-18', '1e-20')): old = mp.mp.dps mp.mp.dps = 70 rows = [] for name, idx in (("J[1,1,0,0,0] = TAD^2", (1, 1, 0, 0, 0)), ("J[1,1,0,1,0] = TAD*FMIX", (1, 1, 0, 1, 0)), ("J[1,1,0,1,1] = FMIX^2", (1, 1, 0, 1, 1))): loc = masters.index(idx) worst = mp.inf for es in eps_list: e = mp.mpf(es) tad = -mp.exp(e * mp.euler) * mp.gamma(e - 1) fmix = _fmix(e, w) closed = {"J[1,1,0,0,0] = TAD^2": tad * tad, "J[1,1,0,1,0] = TAD*FMIX": tad * fmix, "J[1,1,0,1,1] = FMIX^2": fmix * fmix}[name] laurent = state[-2][loc] / e ** 2 + state[-1][loc] / e + state[0][loc] worst = min(worst, agree_digits(laurent, closed)) rows.append((name, worst)) mp.mp.dps = old return rows # ==================================== main =================================== def main(): if _ARGS.sum_rule is not None: _sum_rule_tier() # the 6 zeta_3 spot gate on the closed form (2026-09-11); raises SystemExit _fc, _fc_note = _fc_eligible() if _fc_note: print(f"[fast-cache] {_fc_note}") if _fc is not None: _fc_fast_path(_fc) # raises SystemExit(0) on success mp.mp.dps = 200 ORACLE_KP = mp.mpf(ORACLE_KP_STR) ORACLE_K1 = mp.mpf(ORACLE_K1_STR) ORACLE_SE = mp.mpf(ORACLE_SE_STR) GT_AMF = mp.mpf(GT_AMF_STR) RHO_W12 = mp.mpf(RHO_W12_STR) AMF_BOX_M = mp.mpf(AMF_BOX_M_STR) AMF_BOX_DOT = mp.mpf(AMF_BOX_DOT_STR) mp.mp.dps = 90 A1 = -2 * mp.pi print("LBL3KP compliant final form -- parent IBP-partial-fractioned onto kite-family") print("masters (runtime Picard-Fuchs series rho) x closed dilog-box stub kernels,") print("one-fold in the insertion mass w. No node cache (round-1 cache DELETED).") print(f"settings: QUAD_DPS={QUAD_DPS} LEVEL={LEVEL} DPS={DPS} NORD={NORD}") print("guards (derived from QUAD_DPS, all bounds << gate target):") print(f" U_MIN={mp.nstr(U_MIN,3)} (dropped < {mp.nstr(10*U_MIN**2*mp.mpf('0.2')/mp.pi,3)})") print(f" V_MIN={mp.nstr(V_MIN,3)} (dropped < {mp.nstr(V_MIN**2*mp.fabs(mp.log(V_MIN))*mp.mpf('0.02')/mp.pi,3)})") print(f" W_MAX={mp.nstr(W_MAX,3)} (dropped < {mp.nstr(2*mp.mpf('0.57')*mp.log(W_MAX)/W_MAX**2,3)})\n") # --- [0] the exact rational connection + the single boundary seed -------- print(f"[0] [{time.time()-T0:6.1f}s] loading exact rational kite IBP connection ...") funcs, masters = load_kite_de() M0, d_seed, seed_dps = load_kite_boundary(masters) print(" 8 masters, eps-graded to eps^2; seed = w=5 Laurent vector " "(DERIVED, AMFlow-free:") print(f" lbl3kp-w5-derived.json, emitted at dps {seed_dps} by the " f"lbl3kp-w5-seed.py shim --") print(" vacuum closed forms + Frobenius at w=0 + exact-DE march; " "rerunnable here)") print(f" seed vs held-out retired AMFlow w=5 vector " f"(goal_digits=140): >= {mp.nstr(d_seed, 5)} d agreement") # --- [1] rho(w) as a PF-series sweep over (1, inf) ----------------------- print(f"[1] [{time.time()-T0:6.1f}s] building rho(w) by one adaptive local-series sweep ...") rho = RhoPF(funcs, M0) # held-out check: rho(12) vs the independent AMFlow w=12 run mp.mp.dps = DPS r12 = rho(mp.mpf(12)) d12 = agree_digits(r12, RHO_W12) r12_cap = _sig_digits(RHO_W12_STR) print(f" rho(12) live = {mp.nstr(r12, min(50, DPS-4))}") print(f" held-out oracle = {mp.nstr(RHO_W12, 50)}") print(f" MEASURED agreement = {mp.nstr(d12, 6)} d " f"(limited by DPS={DPS}; oracle string carries {r12_cap} digits;") print(" archived dps=130 record: 99.95 d)\n") # the CLOSED-FORM density at w = 5, 12, 100 (2026-09-11, section C2): # gated against the AMFlow references at the sibling's floors and # against the live series density above at the rho bar of the point tiers if _ARGS.density == 'closed': bar_c = min(30, QUAD_DPS - 10) closed_density_gates({'rho(5)': rho(mp.mpf(5)), 'rho(12)': r12, 'rho(100)': rho(mp.mpf(100))}, f"the live transport above, DPS={DPS}", bar12=bar_c, bar_series=bar_c) # --- [2] closed-form Gamma/2F1 checks of the non-elliptic masters -------- print(f"[2] [{time.time()-T0:6.1f}s] IBP-onto-kite-masters exhibit: transported masters vs") print(" Gamma/2F1 closed forms at the w=12 checkpoint (small-eps Laurent probe):") if mp.mpf(12) in rho.checkpoints: for name, d in closed_masters_check(rho.checkpoints[mp.mpf(12)], 12, masters): print(f" {name:26s} {mp.nstr(d, 5)} d") print(" (the remaining masters are the Gamma_1(6) sunrise block --") print(" the elliptic content -- and the m00 sunset + top one-fold)\n") else: print(" checkpoint missed -- skipped\n") # --- [3] kinematic point + the three stub kernels ------------------------ mp.mp.dps = 90 point_mode = _ARGS.point is not None if point_mode: S, T = mp.mpf(_ARGS.point[0]), mp.mpf(_ARGS.point[1]) else: S, T = mp.mpf(-1), mp.mpf(-1) / 3 M2 = mp.mpf(1) if not (S < 0 and T < 0 and -S - T < 4 * M2): raise SystemExit("DOMAIN: need s<0, t<0, u=-s-t<4*m2 " "(deep Euclidean, units m2=1)") print(f"[3] [{time.time()-T0:6.1f}s] kernels at (s,t)=({mp.nstr(S, 8)},{mp.nstr(T, 8)}): " "closed dilog box, runtime") print(" Chebyshev evaluation device + exact stub partial fractions:") K = BoxKernel(S, T, M2, QUAD_DPS) # kernel constants: K(1) direct vs Chebyshev fit; K'(1) Richardson FD vs fit kd = K.kd mp.mp.dps = kd + 20 Box_m_direct = K._direct(mp.mpf(1)) Box_dot_fd = box_dot_fd(K, dps=max(120, kd + 20)) d_m = agree_digits(K.tay1[0], Box_m_direct) d_d = agree_digits(K.tay1[1], Box_dot_fd) print(f" K(1) = {mp.nstr(Box_m_direct, 40)}") print(f" K'(1) = {mp.nstr(Box_dot_fd, 40)}") print(f" cross-route: fit-vs-direct K(1) {mp.nstr(d_m,5)} d, " f"fit-vs-FD K'(1) {mp.nstr(d_d,5)} d") if d_m < QUAD_DPS or d_d < min(QUAD_DPS, 60): raise SystemExit("kernel-constant cross-route check below target -- STOP") if not point_mode: mp.mp.dps = 160 da = agree_digits(Box_m_direct, AMF_BOX_M) db = agree_digits(Box_dot_fd, AMF_BOX_DOT) capm, capd = _sig_digits(AMF_BOX_M_STR), _sig_digits(AMF_BOX_DOT_STR) print(f" vs archived AMFlow nu=[1,1,1,1]: {mp.nstr(da,5)} d " f"(oracle string {capm} digits)") print(f" vs archived AMFlow nu=[2,1,1,1]: {mp.nstr(db,5)} d " f"(oracle string {capd} digits; FD route truncation ~65-70 d)") if da < 45 or db < 45: raise SystemExit("kernel constants vs archived AMFlow below 45 d -- STOP") mp.mp.dps = 90 # exact PF stub kernels, Taylor route inside |w-1|<0.3 with tay1 slices Box_m, Box_dot = K.tay1[0], K.tay1[1] KP = KernelPF(K, Box_m, Box_dot, npow=2) K1 = KernelPF(K, Box_m, Box_dot, npow=1) KP.seam_check("K_P") K1.seam_check("K_1") # --- [4] the three one-folds (shared memoized rho) ----------------------- panels = lbl3kp_panels(A1) kn = int(1.7 * QUAD_DPS) + 10 jobs = [("I_KP", "(nu3=2, parent)", KP, K.tay1[2:]), ("I_K1", "(nu3=1, stub) ", K1, K.tay1[1:]), ("I_SE", "(nu3=0, LBL3SE)", K, K.tay1)] results = {} for name, tag, kern, Kc in jobs: print(f"[4] [{time.time()-T0:6.1f}s] {name} {tag} = (1/pi) int rho K dw -- " "subtracted tanh-sinh ladder:") prev = None vals = {} for lev in range(3, LEVEL + 1): t0 = time.time() vals[lev] = disp_subtracted(rho, kern, panels, dps=QUAD_DPS, kernel_nmax=kn, maxdegree=lev, Kc_pre=Kc) mp.mp.dps = 90 step = mp.nstr(vals[lev] - prev, 4) if prev is not None else '-' print(f" L{lev}: {name} = {mp.nstr(vals[lev], 44)} " f"step={step} ({time.time()-t0:.1f}s)") prev = vals[lev] for k in sorted(vals): # 2026-09-06 (Q20e): fail CLOSED on a non-finite level if not mp.isfinite(vals[k]): # value, by name; the downstream step and raise RuntimeError( # gate tests are False for a NaN (unnamed) f"[refine] L{k}: {name} (the subtracted tanh-sinh level value) " "is not finite (nan/inf) -- fail closed") results[name] = vals # --- [5] gate (demo point) or self-consistency (other points) ------------ mp.mp.dps = 90 I_KP, I_K1, I_SE = (results[n][LEVEL] for n in ("I_KP", "I_K1", "I_SE")) if not point_mode: print(f"\n[5] [{time.time()-T0:6.1f}s] gate vs held-out oracles " "(d = -log10|f-oracle|/|oracle|, recomputed now):") rows = [ ("I_KP vs independent 9-prop AMFlow (Neville)", I_KP, ORACLE_KP, ORACLE_KP_STR, "the twenty-sample x_order-400/500 grid peel: order-convergence bar 68.73 relative (69.68 absolute), leave-one-out bar 74.45 relative (75.40 absolute)"), ("I_K1 vs independent 9-prop AMFlow (Neville)", I_K1, ORACLE_K1, ORACLE_K1_STR, "the twenty-sample x_order-400/500 grid peel: order-convergence bar 68.56 relative (69.23 absolute), leave-one-out bar 74.28 relative (74.95 absolute)"), ("I_SE vs independent 9-prop AMFlow (Neville)", I_SE, ORACLE_SE, ORACLE_SE_STR, "the twenty-sample x_order-400/500 grid peel: order-convergence bar 68.17 relative (68.45 absolute), leave-one-out bar 73.89 relative (74.17 absolute)"), ("I_SE vs prior 8-prop-family AMFlow eps-grid ", I_SE, GT_AMF, GT_AMF_STR, "the 42-digit artifact of the earlier 8-propagator-family eps-grid: its own ~40 d floor"), ] ok = True _fc_row_ds = [] thresh = 30.0 if QUAD_DPS >= 40 else QUAD_DPS - 5.0 if QUAD_DPS < 40: print(f" NOTE: QUAD_DPS={QUAD_DPS} < 40 targets fewer digits than the " f">=30d campaign gate;\n verdict threshold here is " f"{thresh:.0f} d (requested-precision consistency), the 30d gate " "runs at default --dps 40") for label, ours, oracle, lit, note in rows: d = float(agree_digits(ours, oracle)) _fc_row_ds.append(d) print(f" {label}") print(f" computed : {mp.nstr(ours, 44)}") print(f" oracle : {mp.nstr(oracle, 44)}") print(f" agreement: {d:.2f} d (oracle literal {_sig_digits(lit)} digits; {note})") ok = ok and d >= thresh print("\n (the three Neville rows read this script's own digits at the run's --dps") print(" against 60-digit literals whose grid bars sit at 68-75 d; the 8-prop-family") print(" row saturates at its 42-digit literal; dps growth also shows on rho(12) and") print(" the 138-digit kernel-constant checks above)") verdict = ("PASS -- all three integrals recomputed from the runtime PF rho + " "closed kernels agree with their held-out oracles to >= 30 digits" if ok else "FAIL -- see above") elif _POINT_TAG is not None: # the tagged tier: the level-step self-consistency and the rho(12) gate as at # any point, then every target vs the shipped record (bar from the record) print(f"\n[5] [{time.time()-T0:6.1f}s] results at the tagged point {_POINT_TAG} " f"(s,t)=({mp.nstr(S,8)},{mp.nstr(T,8)}):") ok = True for name in ("I_KP", "I_K1", "I_SE"): vals = results[name] step_rel = mp.fabs(vals[LEVEL] - vals[LEVEL - 1]) / max( mp.fabs(vals[LEVEL]), mp.mpf('1e-60')) print(f" {name} = {mp.nstr(vals[LEVEL], min(50, QUAD_DPS + 2))} " f"|L{LEVEL}-L{LEVEL-1}|/|L{LEVEL}| = {mp.nstr(step_rel, 3)}") ok = ok and step_rel < mp.mpf(10) ** (-(QUAD_DPS - 15)) ok = ok and float(d12) >= min(30, QUAD_DPS - 10) ok = point_tag_gate(_POINT_TAG, _POINT_REC, {'I_KP': I_KP, 'I_K1': I_K1, 'I_SE': I_SE}, ok) verdict = (f"PASS -- tagged point {_POINT_TAG}: every target within its bar of " "the shipped reference; rho(12) gate met, level steps converged" if ok else f"FAIL -- tagged point {_POINT_TAG}: see above") else: print(f"\n[5] [{time.time()-T0:6.1f}s] results at (s,t)=({mp.nstr(S,8)}," f"{mp.nstr(T,8)}) [no banked oracle at this point]:") ok = True for name in ("I_KP", "I_K1", "I_SE"): vals = results[name] step_rel = mp.fabs(vals[LEVEL] - vals[LEVEL - 1]) / max( mp.fabs(vals[LEVEL]), mp.mpf('1e-60')) print(f" {name} = {mp.nstr(vals[LEVEL], min(50, QUAD_DPS + 2))} " f"|L{LEVEL}-L{LEVEL-1}|/|L{LEVEL}| = {mp.nstr(step_rel, 3)}") ok = ok and step_rel < mp.mpf(10) ** (-(QUAD_DPS - 15)) ok = ok and float(d12) >= min(30, QUAD_DPS - 10) print(" (ladder steps are self-consistency; the held-out gates that run at") print(" ANY point are rho(12) [section 1] and the kernel cross-route checks)") verdict = ("PASS -- rho held-out gate met, kernel cross-routes passed, ladders " "converged" if ok else "FAIL -- see above") if _ARGS.fastcache_bank and not point_mode and ok: # banked_d derived from the STORED strings (what the fast path # re-checks) at the fast path's parse precision (90) -- see the # lbl3se sibling comment (its fast-start cache section). s_KP = mp.nstr(I_KP, QUAD_DPS + 8) s_K1 = mp.nstr(I_K1, QUAD_DPS + 8) s_SE = mp.nstr(I_SE, QUAD_DPS + 8) s_r12 = mp.nstr(r12, DPS + 10) with mp.workdps(90): bd = {'row0': float(agree_digits(mp.mpf(s_KP), mp.mpf(ORACLE_KP_STR))), 'row1': float(agree_digits(mp.mpf(s_K1), mp.mpf(ORACLE_K1_STR))), 'row2': float(agree_digits(mp.mpf(s_SE), mp.mpf(ORACLE_SE_STR))), 'row3': float(agree_digits(mp.mpf(s_SE), mp.mpf(GT_AMF_STR))), 'd12': float(agree_digits(mp.mpf(s_r12), mp.mpf(RHO_W12_STR)))} _fc_bank({'cached_dps': QUAD_DPS, 'I_KP': s_KP, 'I_K1': s_K1, 'I_SE': s_SE, 'rho12': s_r12, 'Box_m': mp.nstr(Box_m, 60), 'Box_dot': mp.nstr(Box_dot, 60), 'banked_d': bd}, time.time() - T0) elif _ARGS.fastcache_bank: print("[fast-cache] bank REFUSED: gate run did not PASS " "(or --point mode)") print(f"\ntotal wall time: {time.time()-T0:.1f}s") print(verdict) raise SystemExit(0 if ok else 1) if __name__ == '__main__': main()