#!/usr/bin/env python3 r"""lbl3e-master-evaluate.py -- LBL3E (the one-loop massive box with an equal-mass sunrise insertion on one rung; row 27 of the paper's table of integrals): standalone arbitrary-precision evaluator of the FULL elliptic master, assembled at runtime from the closed form written out in lbl3e-expression.md. elliptic-block master (named form, m^2 = 1, w = box-rung virtuality): g^(k)(w) = psi1(w) * [ Eichler_f3(tau(w)) + c1^(k) ] + c2^(k) * psi2(w) sector-15 representative master (the inherited sunrise; AMFlow convention J[1,1,1](d0=2) = -S111(2-2eps, w)), assembled through eps^2: S(2-2eps,w) = (psi1/pi)(w) * Gamma(1+eps)^2 e^{-eps I(f2;q)} * [ F1(eps;q) * B(eps) + F3(eps;q) ] F1 = 1 - (eps/2) I(1;q) + eps^2 I(1,f4;q) + ... F3 = I(1,f3;q) + eps I(1,f3,f2;q) + eps^2 [ I(1,f3,f2,f2;q) + I(1,f4,1,f3;q) ] + ... (arXiv:1704.08895 conventions; Eichler_f3(tau(w)) == I(1,f3;q_C(w)); the eps^0 row is exactly J^(0) = -(psi1/pi) [ B^(0) + Eichler_f3 ].) Everything is computed AT RUNTIME from q-series / Picard-Fuchs series / 2F1 eps-series: * psi1/pi two independent ways: Gamma_1(6) q-series 2*sqrt3*(e1+e2) at the Newton-inverted nome, AND the Picard-Fuchs Frobenius solution varpi0 of L_sun = t(t-1)(t-9) d^2 + (3t^2-20t+9) d + (t-3) Taylor-transported along the negative real axis; exact relation psi1/pi = (2/sqrt3)*varpi0. * psi2/pi = tau_C * psi1/pi, with tau_C from TWO independent routes: (a) Newton nome: tau = log(q_C)/(2 pi i), log q = ln|q| + i pi (q<0); (b) Picard-Fuchs SECOND (log-Frobenius) solution y2 = varpi0*ln t + g(t), recurrence for g derived from L_sun at runtime; tau = (y2/varpi0 - ln 9)/(2 pi i). * Eichler_f3 in the convention of the reference values = the one-fold cusp-regularized integral I(f3;q) = sum a_n/n q^n of the weight-3 form f3 = 36 sqrt3 (e1^3 - e1^2 e2 - 4 e1 e2^2 + 4 e2^3) (compared with the six 80-digit reference values at runtime); the sector-15 assembly uses the two-fold I(1,f3;q) = sum a_n/n^2 q^n. * Boundary tower B^(0..4) from the 2F1 form (finite-difference-free): F(eps) = 2F1(-2eps,-eps;1-eps; r3), r3 = e^{2 pi i/3} h(eps) = ( e^{i pi eps/3} F - e^{-i pi eps/3} conj-F ) / i B(eps) = (1/2) 3^{-eps} [ (3/2) h(eps)/eps^2 - pi (Gamma(1+2eps)/Gamma(1+eps)^2)/eps ] with the 2F1 eps-Taylor computed by the Pfaff-transformed term recurrence at u = r3/(r3-1) (|u| = 1/sqrt3) in eps-polynomial arithmetic (no finite differences). B^(0) = 2 Cl2(pi/3) = (3 sqrt3/2) L(chi_-3,2) is the classical positive control; B^(1),B^(2) are the 2F1 eps-derivative partners of the eps^{1,2} rows. WHAT THE DEFAULT RUN CHECKS (every comparison is made at runtime; the reference values are read from row27_data.json and row27_gate_amflow.json beside this script, sha256-pinned below and refused on any mismatch): G1 psi1/pi: Picard-Fuchs route vs q-series route (two codepaths) at the six reference points w = -1/4, -1/2, -1, -2, -3, -5, and vs the 80-digit reference values (the closure record's own cross-check: min 35.7 d, 92.8-116.4 d at 5/6 points). G2 psi2/tau: second-solution Picard-Fuchs route vs Newton-nome route, and psi2/pi vs the reference values. G3 Eichler_f3 vs the 80-digit cusp-regularized reference values. G4 FULL MASTER: the assembled J^(0),J^(1),J^(2) vs 301-digit ball-certified evaluations of the arXiv:1704.08895 literature representation at w = -1, -3, -5 (comparison only; these values enter no computation here). This reference lives in the same homogeneous-period frame as the construction. G5 B^(0) vs classical 2 Cl2(pi/3), and B^(0..4) vs 120-digit reference values. G6 INDEPENDENT AMFLOW VALUES: the recorded closed-form prediction of the row's elliptic content (the uncut equal-mass sunrise, AMFlow convention -S_111(4-2eps,w), eps^-2..eps^2; written to a sha256-stamped file, sha256 71acf39e..., before any AMFlow output existed) vs the recorded AMFlow values (blackbox numeric IBP, precision goals 60 and 40) at w = -7/3 and w = -4 -- two points kept out of the construction (they appear in no fit or anchor set; the construction's reference points are the six above). w = -1 is one of those reference points (and served the one-time convention identification), so it is printed as a control and excluded. Bar >= 30 d. NOTE: G6 compares two sets of RECORDED values carried in row27_gate_amflow.json -- the digits of agreement are computed here, but neither side is recomputed at runtime: the d = 4-2eps prediction is not assembled by this script, whose runtime assembly is the d = 2-2eps master compared in G4. WHAT IS NOT DONE HERE: * The box-dressed sectors 55a-e/63 have no pointwise (c1,c2) as posed: the pointwise named form is structurally excluded for the dressed sectors (the dressing is a one-fold over a moving fibre; it closes via the staged spectral one-fold assembly), and the shared elliptic content carries (c1,c2) = (B^(0),0). This script assembles the sector-15 representative master exactly ((c1,c2) = (B^(0),0), classical) and exposes assemble_master(w, c1, c2, dps) as the general w-pencil solution. * eps^3+ of the assembled master is not built (needs the longer word tower); B^(3),B^(4) ARE computed and compared. * Domain: real Euclidean w < 0 only (singular fibres w in {0,1,9}); points with |q_C| >= 0.92 are refused. Minkowski w needs an i0+ path, not here. * Polylog block (sectors 11/51/59) = textbook one-loop box x bubble, not evaluated here. * lbl3e-evaluate.py (same page) evaluates psi1/pi and B^(0) only; the two scripts are complementary and share no files. CERTIFIED TRUNCATION (fail-closed; every bound below is backed by a RAISE): the q-series truncation is certified by a PROVEN eta-quotient/Eisenstein majorant tail bound (derivation at _qtail_bound_log10); the Newton nome by a residual-based shift bound (non-convergence RAISES); the 2F1 tower by a trailing-window l1 ratio bound with the geometric form guarded at r < 0.95; the Picard-Fuchs transport by per-step trailing-window bounds; the eps^{-2,-1} pole cancellation and the Im residuals of B and of the assembled J are RAISING checks; in --point mode the --check two-precision agreement RAISES below dps-8 (in default mode --check reruns the comparison table at dps+60 under the same bars and prints both worst figures; that rerun must pass). Seed lengths are STARTING SEEDS only; escalation continues the same recurrences x1.7 to a cap and RAISES there. Healthy runs never escalate (measured headrooms 1e4-1e34 when the guards were calibrated, 2026-07-05). EXIT CODES: 0 every check passes; 1 a check fails (this is what --mutate must produce) or a certified bound is not reached; 2 usage error or --point domain refusal; 3 a pinned data file's sha256 does not match (refused before any computation); 4 a required file or mpmath is missing. USAGE python3 lbl3e-master-evaluate.py # default: dps 80, all checks (about ten seconds) python3 lbl3e-master-evaluate.py --dps 120 --point=-7/3 # any Euclidean w < 0 (use --point= for negatives) python3 lbl3e-master-evaluate.py --check # two-precision rule: rerun at dps+60 and diff python3 lbl3e-master-evaluate.py --quick # fewer comparison points python3 lbl3e-master-evaluate.py --mutate # control: c1^(0) = B^(0) -> B^(0)(1+1e-9); MUST exit nonzero Imports: mpmath + stdlib only. Data: row27_data.json, row27_gate_amflow.json (same directory; kept under their record names; sha256-pinned in this script). """ import argparse import hashlib import json import math import os import sys import time EXIT_PASS, EXIT_FAIL, EXIT_USAGE, EXIT_PIN, EXIT_MISSING = 0, 1, 2, 3, 4 try: import mpmath as mp except ImportError: print("MISSING: python3 module mpmath (pip install mpmath); nothing computed") sys.exit(EXIT_MISSING) HERE = os.path.dirname(os.path.abspath(__file__)) DATA_FILE = os.path.join(HERE, "row27_data.json") GATE_FILE = os.path.join(HERE, "row27_gate_amflow.json") DATA_SHA256 = "7d6b5990662e45b1530e8df51bdeacf60fead19d35e6e74ea419dc18a5fd7009" GATE_SHA256 = "951aaac840c4095b9d6f477c662de2783b3a0a2004fe3db38f4b52cb3977de45" # --------------------------------------------------------------------------- # Gamma_1(6) q-series machinery (the same construction as sunrise-evaluate.py # on the sunrise page, extended with f2, f4) # --------------------------------------------------------------------------- def chi_m3(n): r = n % 3 return 1 if r == 1 else (-1 if r == 2 else 0) def pmul(a, b, N): r = [mp.mpf(0)] * (N + 1) for i, ai in enumerate(a): if ai: bi = b[:N + 1 - i] for j, bj in enumerate(bi): if bj: r[i + j] += ai * bj return r def euler_sparse(N): sp, k = [(0, 1)], 1 while k * (3 * k - 1) // 2 <= N: s = -1 if k % 2 else 1 for e in (k * (3 * k - 1) // 2, k * (3 * k + 1) // 2): if e <= N: sp.append((e, s)) k += 1 return sorted(sp) def mul_sparse(a, sp, N): r = [mp.mpf(0)] * (N + 1) for e, s in sp: for i in range(0, N + 1 - e): r[i + e] += s * a[i] return r def div_sparse(a, sp, N): r = [mp.mpf(0)] * (N + 1) for m in range(N + 1): acc = a[m] for e, s in sp[1:]: if e <= m: acc -= s * r[m - e] r[m] = acc return r def build_series(N): """psi1/pi, f2, f3, f4, t(q) as length-(N+1) mpf coefficient lists.""" e1 = [mp.mpf(0)] * (N + 1) e1[0] = mp.mpf(1) / 6 for d in range(1, N + 1): c = chi_m3(d) if c: for m in range(d, N + 1, d): e1[m] += c e2 = [mp.mpf(0)] * (N + 1) for n in range(0, N // 2 + 1): e2[2 * n] = e1[n] s3 = mp.sqrt(3) psi1 = [2 * s3 * (x + y) for x, y in zip(e1, e2)] e11, e12, e22 = pmul(e1, e1, N), pmul(e1, e2, N), pmul(e2, e2, N) f2 = [-6 * (a + 6 * b - 4 * c) for a, b, c in zip(e11, e12, e22)] f3 = [36 * s3 * (a - b - 4 * c + 4 * d) for a, b, c, d in zip(pmul(e11, e1, N), pmul(e11, e2, N), pmul(e12, e2, N), pmul(e22, e2, N))] f4 = [324 * c for c in pmul(e11, e11, N)] # (3 sqrt2 e1)^4 = 324 e1^4 t = [mp.mpf(0)] * (N + 1) t[0] = mp.mpf(1) for d, r in [(6, 8), (1, 4), (2, -8), (3, -4)]: spd = [(d * e, s) for e, s in euler_sparse(N // d)] for _ in range(abs(r)): t = mul_sparse(t, spd, N) if r > 0 else div_sparse(t, spd, N) t = [mp.mpf(0)] + [9 * c for c in t[:-1]] return psi1, f2, f3, f4, t _SER_CACHE = {} def get_series(N): key = (mp.mp.dps, N) if key not in _SER_CACHE: _SER_CACHE[key] = build_series(N) return _SER_CACHE[key] def horner(c, q): acc = mp.mpf(0) for ck in reversed(c): acc = acc * q + ck return acc def _newton(q, tt, tser, dser, iters=60, raise_on_fail=False): """Newton solve of t(q) = tt. raise_on_fail=True (final solves) promotes silent loop exhaustion to a RuntimeError; warm-start tracking calls keep raise_on_fail=False (seeds only).""" tol = mp.mpf(10) ** (-mp.mp.dps + 5) dq = None for _ in range(iters): dq = (horner(tser, q) - tt) / horner(dser, q) q = q - dq if abs(dq) < tol: break else: if raise_on_fail: raise RuntimeError( f"_newton: nome iteration did NOT converge: |dq| = " f"{mp.nstr(abs(dq), 3)} >= tol {mp.nstr(tol, 3)} after {iters} " f"iterations at t = {mp.nstr(tt, 10)} (dps {mp.mp.dps})") return q def nome_for_t(tser, tval): """Euclidean cusp-connected branch: real negative nome. For t < -3 the nome is Newton-tracked along the real t-axis from t = -3 (no singular points on t < 0), avoiding spurious roots of the truncated series.""" dser = [k * tser[k] for k in range(1, len(tser))] if -3 <= tval < 0: return _newton(tval / 9, tval, tser, dser, raise_on_fail=True) q = _newton(mp.mpf(-3) / 9, mp.mpf(-3), tser, dser, raise_on_fail=True) steps = max(40, int(4 * abs(math.log(abs(float(tval)) / 3))) * 10) for k in range(1, steps + 1): u = mp.mpf(k) / steps q = _newton(q, -3 + (tval + 3) * u, tser, dser, iters=12) return _newton(q, tval, tser, dser, raise_on_fail=True) QMAX = 0.92 # --- certified-truncation guards (calibrated 2026-07-05 on healthy runs) ---- TAILWIN = 8 # trailing-window width of the tail certificates TAIL_GUARD_Q = 10 # q-series proven tail: bound < 10^-(dps + 10) TAIL_GUARD_STEP = 6 # PF Taylor step / Frobenius seed value tail: 10^-(digits+6) # (measured healthy 10^-(digits+10.0) -> 1e4 headroom) TAIL_GUARD_STEP_YP = 2 # derivative-series tail: 10^-(digits+2) # (measured healthy 10^-(digits+6.3) -> 1e4.3 headroom) TAIL_GUARD_2F1 = 12 # 2F1 tower value-level bound: 10^-(target dps + 12) POLE_MARGIN = 10 # boundary_B pole-cancellation RAISE at 10^-(ambient-10) # (measured healthy ~10^-ambient -> 1e10 headroom) BIM_MARGIN = 10 # boundary_B Im-residual RAISE at 10^-(ambient-10) # (measured healthy < 10^-2*ambient) JIM_MARGIN = 8 # assembled-J Im RAISE at 10^-(dps-8) relative # (measured healthy ~10^-(dps+15) -> 1e23 headroom) R2F1_RMAX = mp.mpf("0.95") # hyp2f1_r3_es: the geometric tail bound r/(1-r) # is VALID ONLY for r < 1; RAISE (fail-closed) if the # trailing-window sup ratio reaches this (measured # healthy r <= 0.62 at the smallest seed window -> # 1.5x margin; an unguarded r >= 1 would flip the sign # of r/(1-r) and silently ACCEPT a divergent tail) # eta-quotient exponents |r_d| of the built t-series (build_series: t/(9q) = # E(q^6)^8 E(q)^4 / (E(q^2)^8 E(q^3)^4), E = Euler product) _ETA_POWS = ((1, 4), (2, 8), (3, 4), (6, 8)) # K = (pi^2/6) sum_d |r_d|/d = (pi^2/6)(4 + 4 + 4/3 + 4/3): only used to PLACE # the candidate rho (any rho in (qa,1) yields a rigorous bound) _ETA_K = 17.6 _ETA_LOGG_CACHE = {} def _eta_majorant_log10(rho): """log10 G(rho), G(q) = prod_d E(q^d)^(-|r_d|): rigorous numeric sum of log G(rho) = sum_d |r_d| sum_j rho^(dj)/(j(1-rho^(dj))) with its geometric tail ADDED (so the returned value is an upper bound; the 1e-30 cutoff is rigorous BECAUSE the tail term is added). rho is quantized to a double (exact cache key; any rho in (qa,1) yields a valid bound).""" rho = mp.mpf(float(rho)) key = float(rho) if key in _ETA_LOGG_CACHE: return _ETA_LOGG_CACHE[key] tot = mp.mpf(0) for d, c in _ETA_POWS: rd = rho ** d s = mp.mpf(0) j = 1 while True: rdj = rho ** (d * j) term = rdj / (j * (1 - rdj)) s += term if term < mp.mpf("1e-30") or j > 200000: break j += 1 s += rho ** (d * (j + 1)) / ((j + 1) * (1 - rd) ** 2) # rigorous tail tot += c * s out = tot / mp.log(10) _ETA_LOGG_CACHE[key] = out return out def _qtail_bound_log10(M, qa): r"""PROVEN log10 upper bound on the truncation tail |sum_{n>M} c_n q^n| of ALL five built series (psi1, f2, f3, f4, t) at |q| = qa, cut at index M. (A naive sizing |c_n| <= top built coefficient would be FALSE for eta quotients, whose coefficients grow ~10^(1.5 sqrt n); this bound replaces it.) DERIVATION (each step elementary and rigorous): * t/(9q) = prod_{d in {1,2,3,6}} E(q^d)^{r_d}, r = {1:+4, 2:-8, 3:-4, 6:+8}. Coefficient-wise |[q^n](1-q^m)^{+-1}| <= [q^n](1-q^m)^{-1}, so |[q^n] t/(9q)| <= [q^n] G(q) with G = prod_d E(q^d)^{-|r_d|}, a generalized-partition series with coefficients >= 0. Positivity + Cauchy: [q^n] G <= G(rho) rho^{-n} for ANY rho in (0,1) (_eta_majorant_log10 evaluates log G(rho) rigorously). Minimizing over rho reproduces the measured growth law as the proven form |c_n| <= e^{2 sqrt(K n) + o(sqrt n)}, K = (pi^2/6) sum_d |r_d|/d. * Eisenstein block: |e1[n]| = |sum_{d|n} chi_-3(d)| <= sigma_0(n), and e2[n] in {0, e1[n/2]} with sigma_0(n/2) <= sigma_0(n), so e1 and e2 are coefficient-dominated by m(q) = 1/6 + sum_n sigma_0(n) q^n with m(rho) <= 1/6 + rho/(1-rho)^2. Summing |integer coefficients| of the defining monomials: |psi1_n| <= 4 sqrt3 m(rho) rho^{-n}, |f2_n| <= 66 m(rho)^2 rho^{-n}, |f3_n| <= 360 sqrt3 m(rho)^3 rho^{-n}, |f4_n| <= 324 m(rho)^4 rho^{-n}. * Tail sum at |q| = qa < rho: |sum_{n>M} c_n q^n| <= C(rho) (qa/rho)^{M+1} / (1 - qa/rho) per series (t-series shifted by one through its 9q prefactor); the returned bound is max(eta-block, Eisenstein-block) + log10(2) >= their sum. Every candidate rho gives a valid bound, so the min over {1 - sqrt(K/(M+1)), (1+qa)/2, (3+qa)/4} is valid. * Inequalities also checked empirically to N=800 at rho in {0.5, 0.7, 0.85} (worst |c_n|/bound = 0.12, all <= 1). """ best = mp.inf r_opt = 1 - mp.sqrt(mp.mpf(_ETA_K) / (M + 1)) for rho in (r_opt, (1 + qa) / 2, (3 + qa) / 4): rho = mp.mpf(float(rho)) # double-quantized: exact cache key if not (qa < rho < 1): continue c_t, c_e, s = _qtail_consts(rho, qa) best = min(best, max(c_t + M * s, c_e + (M + 1) * s) + mp.log10(2)) return best def _qtail_consts(rho, qa): """(c_t, c_e, s) with lg_tail_t = c_t + M*s, lg_tail_e = c_e + (M+1)*s (both linear in the cut index M at fixed rho; s = log10(qa/rho) < 0).""" s = mp.log10(qa / rho) mrho = mp.mpf(1) / 6 + rho / (1 - rho) ** 2 c_t = (mp.log10(9 * qa) + _eta_majorant_log10(rho) - mp.log10(1 - qa / rho)) C_e = (4 * mp.sqrt(3) * mrho + 66 * mrho ** 2 + 360 * mp.sqrt(3) * mrho ** 3 + 324 * mrho ** 4) c_e = mp.log10(C_e) - mp.log10(1 - qa / rho) return c_t, c_e, s def _q_required_M(qa, tol_exp, ncap): """Smallest-M ESTIMATE with the proven bound < tol_exp, exploiting that the bound is linear in M at fixed rho (closed-form solve per candidate; the adaptive candidate rho = 1-sqrt(K/(M+1)) is fixed-point iterated). Rigor is not needed here: sized_series VERIFIES the returned length with _qtail_bound_log10 and escalates if short. None if ncap insufficient.""" best = None guesses = [(1 + qa) / 2, (3 + qa) / 4] Mg = 400 for _ in range(3): # fixed-point on the adaptive rho r_opt = 1 - mp.sqrt(mp.mpf(_ETA_K) / (Mg + 1)) Ms = [] for rho in guesses + [r_opt]: rho = mp.mpf(float(rho)) if not (qa < rho < 1): continue c_t, c_e, s = _qtail_consts(rho, qa) need = mp.mpf(tol_exp) - mp.log10(2) M1 = (need - c_t) / s # lg_t < tol M2 = (need - c_e) / s - 1 # lg_e < tol Ms.append(int(mp.ceil(max(M1, M2))) + 1) if not Ms: return None Mg = max(8, min(Ms)) best = Mg if best > ncap or _qtail_bound_log10(min(best, ncap), qa) >= tol_exp: if _qtail_bound_log10(ncap, qa) >= tol_exp: return None return ncap return best def sized_series(tval, dps): """(psi1,f2,f3,f4,tser,q,cert): q-series sized by the PROVEN eta-majorant tail bound (_qtail_bound_log10), certified < 10^-(dps+TAIL_GUARD_Q) at the returned length. N0 = max(400, 2.2*dps) survives as the STARTING SEED only (speed, never answer); escalation rebuilds the SAME recurrences (build_series) at x1.7 up to NCAP = 64*N0; RuntimeError at cap. Also certifies the Newton nome: raise_on_fail final solves + a residual-based nome-shift bound (|t(q)-t| at +10 dps + series tail, over |t'(q)|).""" N0 = max(400, int(2.2 * dps)) NCAP = 64 * N0 tol_exp = -(dps + TAIL_GUARD_Q) psi1, f2, f3, f4, tser = get_series(N0) q = nome_for_t(tser, tval) qa = abs(q) if qa >= QMAX: raise ValueError(f"|q_C(w)| = {float(qa):.3f} >= {QMAX}: too close to a " "cusp (w -> 0^-, -inf); outside the series domain") Mreq = _q_required_M(qa, tol_exp, NCAP) if Mreq is None: raise RuntimeError( f"sized_series: proven eta-majorant tail bound cannot reach tol " f"1e{tol_exp} within cap at t = {mp.nstr(tval, 10)}: |q| = " f"{float(qa):.4f}, bound(NCAP) = 1e{float(_qtail_bound_log10(NCAP, qa)):.1f}, " f"NCAP = {NCAP} (= 64 x seed {N0})") N = max(N0, Mreq) bexp = None for _ in range(8): if len(tser) - 1 < N: # EXACT continuation: build_series re-runs the same recurrences # to the longer length; re-polish the nome on the longer t-series psi1, f2, f3, f4, tser = get_series(N) dser = [k * tser[k] for k in range(1, len(tser))] q = _newton(q, tval, tser, dser, raise_on_fail=True) qa = abs(q) bexp = _qtail_bound_log10(len(tser) - 1, qa) if bexp < tol_exp: break if N >= NCAP: raise RuntimeError( f"sized_series: certified tail bound 1e{float(bexp):.1f} >= tol " f"1e{tol_exp} at t = {mp.nstr(tval, 10)}, |q| = {float(qa):.4f}, " f"N = {len(tser) - 1} (cap {NCAP})") N = min(NCAP, int(N * 1.7) + 1) else: raise RuntimeError( f"sized_series: escalation loop exhausted at t = {mp.nstr(tval, 10)}: " f"bound 1e{float(bexp):.1f} >= tol 1e{tol_exp}, N = {len(tser) - 1} " f"(cap {NCAP})") # certified nome-shift bound: residual re-evaluated at +10 dps so the # horner roundoff floor sits ~10 digits below the shift tolerance with mp.workdps(mp.mp.dps + 10): dser = [k * tser[k] for k in range(1, len(tser))] resid = abs(horner(tser, q) - tval) tprime = abs(horner(dser, q)) shift = (resid + mp.mpf(10) ** bexp) / tprime shift_tol = mp.mpf(10) ** (-(dps + 2)) if not (shift < shift_tol): raise RuntimeError( f"sized_series: certified nome-shift bound {mp.nstr(shift, 3)} >= " f"tol {mp.nstr(shift_tol, 3)} at t = {mp.nstr(tval, 10)} (Newton " f"residual {mp.nstr(resid, 3)}, series tail 1e{float(bexp):.1f}, " f"|t'(q)| = {mp.nstr(tprime, 3)})") cert = {"q_tail_log10": float(bexp), "q_tail_tol_log10": tol_exp, "N": len(tser) - 1, "nome_shift": shift, "nome_shift_tol": shift_tol} return psi1, f2, f3, f4, tser, q, cert # --------------------------------------------------------------------------- # Picard-Fuchs Frobenius machinery (the same construction as lbl3e-evaluate.py # on this page, extended with the second/log solution): L_sun = t(t-1)(t-9) y'' + (3t^2-20t+9) y' + (t-3) y # --------------------------------------------------------------------------- def _frob_extend(a, b, upto): """EXACT continuation of the paired Frobenius recurrences to index upto.""" for m in range(len(a) - 1, upto): am1 = a[m - 1] if m >= 1 else mp.mpf(0) a.append(((10 * m * m + 10 * m + 3) * a[m] - m * m * am1) / (9 * (m + 1) ** 2)) bm1 = b[m - 1] if m >= 1 else mp.mpf(0) # source S(t) = -[2(t^2-10t+9) y1' + (2t-10) y1], coefficient s_m: s_m = -(18 * (m + 1) * a[m + 1] - (20 * m + 10) * a[m] + 2 * m * am1) b.append(((10 * m * m + 10 * m + 3) * b[m] - m * m * bm1 + s_m) / (9 * (m + 1) ** 2)) def frob_pair_at_zero(t, digits): """(y1, y1', yL, yL', err_val, err_der) at t in (-1/4, 0): y1 = varpi0 Frobenius series, yL = y1*ln|t| + g with L_sun[y1 ln t + g] = 0, g(0) = 0 (MUM log solution; the +i pi y1 piece of the t<0 branch is added analytically later). nmax = f(digits) is a STARTING SEED only; the truncation is certified by the trailing-TAILWIN-window bound max(|a_n t^n|, |b_n t^n|) * r/(1-r) with r = |t| certified by the Frobenius radius (regular singular points of L_sun at {0,1,9} => the a- and b-series both converge on |t| < 1; here |t| <= 1/8). Escalation = EXACT continuation of the same recurrences x1.7 up to 16 x seed; RuntimeError at cap. Measured healthy bound 1e-(digits+31) at t=-1/8 vs tol 1e-(digits+6).""" nmax = int(digits / (-math.log10(abs(float(t))))) + 40 ncap = 16 * nmax a = [mp.mpf(1)] b = [mp.mpf(0)] r = abs(t) fac = r / (1 - r) tol_v = mp.mpf(10) ** (-(digits + TAIL_GUARD_STEP)) tol_d = mp.mpf(10) ** (-(digits + TAIL_GUARD_STEP_YP)) target = nmax while True: _frob_extend(a, b, target) M = len(a) - 1 wmax = max(max(abs(a[n]), abs(b[n])) * r ** n for n in range(max(0, M - (TAILWIN - 1)), M + 1)) e_val = wmax * fac e_der = (wmax / r) * (M * fac + r / (1 - r) ** 2) if e_val < tol_v and e_der < tol_d: break if target >= ncap: raise RuntimeError( f"frob_pair_at_zero: certified tail bound not reached at cap: " f"t = {mp.nstr(t, 10)}, bound_val = {mp.nstr(e_val, 3)} (tol " f"{mp.nstr(tol_v, 3)}), bound_der = {mp.nstr(e_der, 3)} (tol " f"{mp.nstr(tol_d, 3)}), N = {M} (cap {ncap})") target = min(ncap, int(target * 1.7) + 1) y1 = mp.mpf(0); y1p = mp.mpf(0); g = mp.mpf(0); gp = mp.mpf(0) for n in range(len(a) - 1, -1, -1): y1 = y1 * t + a[n] g = g * t + b[n] for n in range(len(a) - 1, 0, -1): y1p = y1p * t + n * a[n] gp = gp * t + n * b[n] lt = mp.log(abs(t)) # error transport into (yL, yL'): |lt|, 1/|t| are the exact prefactors e_L = e_val * (abs(lt) + 1) e_Lp = e_der * (abs(lt) + 1) + e_val / abs(t) return y1, y1p, y1 * lt + g, y1p * lt + y1 / t + gp, \ max(e_val, e_L), max(e_der, e_Lp) def _taylor_coeffs_extend(p, qq, r, c, upto): """EXACT continuation of the local L_sun Taylor recurrence: extends c (c[0], c[1] given) so that len(c) == upto.""" for n in range(len(c) - 2, upto - 2): s = mp.mpf(0) for j in range(1, 4): k = n - j + 2 if 0 <= k < len(c): s += p[j] * k * (k - 1) * c[k] for j in range(3): k = n - j + 1 if 0 <= k < len(c): s += qq[j] * k * c[k] for j in range(2): k = n - j if 0 <= k < len(c): s += r[j] * c[k] c.append(-s / (p[0] * (n + 2) * (n + 1))) def taylor_step(t0, y, yp, h, nterm, rloc, tol_y, tol_yp): """One CERTIFIED Taylor step of L_sun from t0 to t0+h. nterm is the STARTING SEED; the truncation is certified by the trailing-TAILWIN-window bound max|c_n h^n| * r/(1-r), r = |h|/rloc certified by the step rule (|h| <= 0.45 * rloc, rloc = distance to the nearest solution singularity in {0,1,9}); the derivative tail carries the exact n-weighted geometric factor sum_{n>M} n r^n <= M r/(1-r) + r/(1-r)^2 on the same envelope. Escalation = EXACT continuation of the same recurrence (the c-list is extended in place) x1.7 up to 16 x seed; RuntimeError at cap. Returns (ynew, ypnew, bound_y, bound_yp). Caveat: the trailing-window max is the geometric envelope constant; polynomial / log-singularity factors are absorbed by the guard digits + escalation.""" p = [t0 ** 3 - 10 * t0 ** 2 + 9 * t0, 3 * t0 ** 2 - 20 * t0 + 9, 3 * t0 - 10, mp.mpf(1)] qq = [3 * t0 ** 2 - 20 * t0 + 9, 6 * t0 - 20, mp.mpf(3)] r = [t0 - 3, mp.mpf(1)] rr = abs(h) / rloc assert rr < 1 fac = rr / (1 - rr) ncap = 16 * nterm target = nterm c = [y, yp] ha = abs(h) while True: _taylor_coeffs_extend(p, qq, r, c, target) M = len(c) - 1 wmax = max(abs(c[n]) * ha ** n for n in range(max(0, M - (TAILWIN - 1)), M + 1)) bound_y = wmax * fac bound_yp = (wmax / ha) * (M * fac + rr / (1 - rr) ** 2) if bound_y < tol_y and bound_yp < tol_yp: break if target >= ncap: raise RuntimeError( f"taylor_step: certified tail bound not reached at cap: t0 = " f"{mp.nstr(t0, 10)}, h = {mp.nstr(h, 6)}, bound_y = " f"{mp.nstr(bound_y, 3)} (tol {mp.nstr(tol_y, 3)}), bound_yp = " f"{mp.nstr(bound_yp, 3)} (tol {mp.nstr(tol_yp, 3)}), N = {M + 1} " f"(cap {ncap})") target = min(ncap, int(target * 1.7) + 1) ynew = mp.mpf(0); ypnew = mp.mpf(0) for n in range(len(c) - 1, -1, -1): ynew = ynew * h + c[n] for n in range(len(c) - 1, 0, -1): ypnew = ypnew * h + n * c[n] return ynew, ypnew, bound_y, bound_yp def pf_transport(t_target, digits): """(varpi0, yL, err_y1, err_yL) at real t_target < 0 by CERTIFIED Taylor continuation on the negative axis. The old nterm = f(digits) formula is the per-step seed only; every step is checked by its certified trailing- window tail bound (see taylor_step). err_* accumulate the per-step certified bounds at value level (first-order (c0,c1) error transport: E += |h|*Ep + bound_y, Ep += bound_yp; higher-order flow sensitivity is absorbed by the guard digits -- the stated caveat).""" t_target = mp.mpf(t_target) t0 = t_target if t_target > mp.mpf(-1) / 8 else mp.mpf(-1) / 8 y1, y1p, yL, yLp, E, Ep = frob_pair_at_zero(t0, digits) EL, ELp = E, Ep tol_y = mp.mpf(10) ** (-(digits + TAIL_GUARD_STEP)) tol_yp = mp.mpf(10) ** (-(digits + TAIL_GUARD_STEP_YP)) eps_land = mp.mpf(10) ** (-digits) while t0 > t_target + eps_land: dist = min(abs(t0), abs(t0 - 1), abs(t0 - 9)) h = -min(mp.mpf("0.45") * dist, t0 - t_target) # solution-radius: yL has a log singularity at t=0, so include |t0| rloc = min(abs(t0 - 1), abs(t0 - 9), abs(t0)) ratio = float(abs(h) / rloc) nterm = int(digits / (-math.log10(min(ratio, 0.45)))) + 30 y1, y1p, by, byp = taylor_step(t0, y1, y1p, h, nterm, rloc, tol_y, tol_yp) E, Ep = E + abs(h) * Ep + by, Ep + byp yL, yLp, bL, bLp = taylor_step(t0, yL, yLp, h, nterm, rloc, tol_y, tol_yp) EL, ELp = EL + abs(h) * ELp + bL, ELp + bLp t0 += h return y1, yL, E, EL # --------------------------------------------------------------------------- # Iterated q-integrals with tangential-base-point (cusp) regularization. # Value V(q) = sum_j F_j(q) L^j, L = log q = ln|q| + i pi (Euclidean branch, # q < 0; equivalently L = 2 pi i tau_C). Constant terms integrate to L-powers. # --------------------------------------------------------------------------- def integrate_kernel(f, V, N): """W = int_0^q f V dq'/q', V = {j: coeff-list}, tangential-base regularized.""" out = {} for j, Fj in V.items(): U = Fj if f is None else pmul(f, Fj, N) # f = None means kernel 1 if U[0]: lvl = out.setdefault(j + 1, [mp.mpf(0)] * (N + 1)) lvl[0] += U[0] / (j + 1) fac = 1 for i in range(0, j + 1): if i > 0: fac *= -(j - i + 1) lvl = out.setdefault(j - i, [mp.mpf(0)] * (N + 1)) for n in range(1, N + 1): if U[n]: lvl[n] += U[n] * fac / mp.mpf(n) ** (i + 1) return out class WordTable: """Iterated integrals I(letters...; q) with suffix caching. Letters: 'one', 'f2', 'f3', 'f4' (q-series from build_series).""" def __init__(self, f2, f3, f4, N): self.k = {"one": None, "f2": f2, "f3": f3, "f4": f4} self.N = N self.cache = {(): {0: [mp.mpf(1)] + [mp.mpf(0)] * N}} def logq_series(self, word): word = tuple(word) if word not in self.cache: inner = self.logq_series(word[1:]) self.cache[word] = integrate_kernel(self.k[word[0]], inner, self.N) return self.cache[word] def value(self, word, q, L): V = self.logq_series(word) acc = mp.mpc(0) for j, Fj in V.items(): acc += horner(Fj, q) * L ** j return acc # --------------------------------------------------------------------------- # eps-series helpers (truncated Taylor in eps, complex coefficients) # --------------------------------------------------------------------------- def es_mul(a, b, K): r = [mp.mpc(0)] * K for i in range(K): if a[i]: for j in range(K - i): r[i + j] += a[i] * b[j] return r def es_inv(a, K): r = [mp.mpc(0)] * K r[0] = 1 / a[0] for n in range(1, K): s = mp.mpc(0) for j in range(1, n + 1): s += a[j] * r[n - j] r[n] = -s / a[0] return r def es_exp(a, K): assert a[0] == 0 r = [mp.mpc(0)] * K r[0] = mp.mpc(1) for n in range(1, K): # r' = a' r => n r_n = sum_{j>=1} j a_j r_{n-j} s = mp.mpc(0) for j in range(1, n + 1): s += j * a[j] * r[n - j] r[n] = s / n return r def gamma1p_es(scale, K): """Gamma(1 + scale*eps) as an eps-series (lnGamma Taylor, no finite diff).""" ln = [mp.mpc(0)] * K if K > 1: ln[1] = -mp.euler * scale for k in range(2, K): ln[k] = (-1) ** k * mp.zeta(k) / k * scale ** k return es_exp(ln, K) # --------------------------------------------------------------------------- # Boundary tower B^(k) from the printed 2F1 form (Pfaff series at |u|=1/sqrt3) # --------------------------------------------------------------------------- def hyp2f1_r3_es(K, tol_abs=None): """eps-Taylor (length K) of 2F1(-2eps,-eps;1-eps;r3), r3 = e^{2 pi i/3}, via Pfaff: (1-r3)^eps * 2F1(1+eps,-eps;1-eps;u), u = r3/(r3-1). Term recurrence t_{n+1} = t_n (n+1+eps)(n-eps)/((n+1-eps)(n+1)) * u. CERTIFIED tail bound (returns (F, tail_bound)): the truncated eps-polynomial arithmetic is per-coefficient EXACT (order-i output depends only on orders <= i of the factors), and the l1 norm ||t||_1 = sum_i |t[i]| is submultiplicative under truncated products, so the term recurrence gives the certified ratio ||t_{n+1}||_1 <= |u| (n+2)/n * ||t_n||_1 (|u| = 1/sqrt3; ||n+1+eps||_1 = n+2, ||n-eps||_1 = n+1, ||1/(n+1-eps)||_1 <= 1/n; checked empirically over 400 terms, worst ratio quotient 0.9946). Per-coefficient tail: |sum_{n>M} t_n[i]| <= max(||t_n||_1, trailing TAILWIN) * r/(1-r), r = sup over the trailing window of |u|(k+2)/k (attained at the window start; decreasing in k => dominates the whole tail). The geometric form r/(1-r) is VALID ONLY for r < 1: guarded FAIL-CLOSED at R2F1_RMAX = 0.95 (RAISE with named diagnostics; an unguarded r >= 1 would give a negative "bound" that accepts divergent tails). The prefactor (1-r3)^eps multiplies the bound by ||pref||_1 <= e^{|log(1-r3)|}. The M = f(dps) formula is the STARTING SEED only; escalation continues the SAME recurrence x1.7 to 64 x seed; RuntimeError at cap. tol_abs = None sets the bar at 10^-(ambient dps).""" r3 = mp.exp(2j * mp.pi / 3) u = r3 / (r3 - 1) ua = abs(u) lg = mp.log(1 - r3) # = ln sqrt3 - i pi/6 (principal) pfac = mp.e ** abs(lg) # >= ||pref||_1 if tol_abs is None: tol_abs = mp.mpf(10) ** (-mp.mp.dps) M0 = int(mp.mp.dps / (-mp.log10(ua))) + 40 mcap = 64 * M0 t = [mp.mpc(0)] * K t[0] = mp.mpc(1) f = list(t) n = 0 target = M0 win = [] while True: while n < target: t1 = [(n + 1) * t[i] + (t[i - 1] if i else 0) for i in range(K)] # *(n+1+eps) t1 = [n * t1[i] - (t1[i - 1] if i else 0) for i in range(K)] # *(n-eps) inv = [mp.mpc(1) / (n + 1)] for _ in range(K - 1): inv.append(inv[-1] / (n + 1)) # 1/(n+1-eps) t = es_mul(t1, inv, K) t = [c * u / (n + 1) for c in t] f = [a + b for a, b in zip(f, t)] n += 1 win.append(sum(abs(c) for c in t)) if len(win) > TAILWIN: win.pop(0) mstart = n - (len(win) - 1) # first term index inside the window # sup of the certified l1 ratio over the trailing window AND the whole # tail: r_k = |u|(k+2)/k is strictly decreasing in k, so the window # sup (taken explicitly) dominates every tail ratio. if not win or mstart < 1: r = mp.inf # empty/degenerate window (k <= 0 ratio # undefined) -> no certificate else: r = max(ua * (k + 2) / k for k in range(mstart, n + 1)) if not (r < R2F1_RMAX): # FAIL-CLOSED guard: # for r >= 1 the geometric tail sum diverges and r/(1-r) goes # NEGATIVE -- the old unguarded formula could return a negative # "bound" < tol and silently certify a divergent tail. Escalate # while the ratio can still shrink (r_k -> |u| as k grows); RAISE # at the cap, or immediately when even the asymptotic ratio |u| # is >= R2F1_RMAX (escalation can never repair that). if ua >= R2F1_RMAX or target >= mcap: raise RuntimeError( f"hyp2f1_r3_es: trailing-window sup ratio r = " f"{mp.nstr(r, 6)} >= R2F1_RMAX = {mp.nstr(R2F1_RMAX, 3)} " f"at window start m = {mstart} (|u| = {mp.nstr(ua, 6)}, " f"n = {n}, seed {M0}, cap {mcap}, dps {mp.mp.dps}): " f"geometric tail bound r/(1-r) is INVALID -- refusing to " f"certify (fail-closed)") target = min(mcap, int(target * 1.7) + 1) continue bound = max(win) * r / (1 - r) * pfac if bound < tol_abs: break if target >= mcap: raise RuntimeError( f"hyp2f1_r3_es: certified tail bound {mp.nstr(bound, 3)} >= " f"tol {mp.nstr(tol_abs, 3)} at M = {n} (cap {mcap}, seed {M0}, " f"dps {mp.mp.dps})") target = min(mcap, int(target * 1.7) + 1) pref = es_exp([mp.mpc(0)] + [lg] + [mp.mpc(0)] * (K - 2), K) return es_mul(pref, f, K), bound def boundary_B(nB, tol_dps=None): """B^(0..nB-1) + pole-cancellation residuals (must vanish) from B(eps) = (1/2) 3^{-eps} [ (3/2) h/eps^2 - pi (Gamma(1+2eps)/Gamma(1+eps)^2)/eps ]. Returns (Bre, polres, imres, cert). The 2F1 series tail is certified in hyp2f1_r3_es and propagated to B level by the l1 norms of the EXACT prefactors: h picks up 2 e^{pi/3} (||e^{+-i pi eps/3}||_1 <= e^{pi/3}), the bracket 3/2, 3^{-eps} a factor ||p3||_1 <= 3, the final /2 one half => bound_B = (9/2) e^{pi/3} * bound_F, checked at 10^-(tol_dps) (callers pass target dps + TAIL_GUARD_2F1; measured healthy bound_B ~ 10^-(ambient+6.3) => 1e9-1e34 headroom, never fires on healthy runs). RAISING checks: the eps^{-2,-1} pole-cancellation residual and the Im residual must sit at the roundoff floor -- measured healthy polres ~ 10^-ambient, imres < 10^-2*ambient -- RAISE above 10^-(ambient-10).""" K = nB + 4 D = mp.mp.dps if tol_dps is None: tol_dps = D fac_B = mp.mpf(9) / 2 * mp.e ** (mp.pi / 3) F, bF = hyp2f1_r3_es(K, tol_abs=mp.mpf(10) ** (-tol_dps) / fac_B) bound_B = bF * fac_B Fc = [mp.conj(c) for c in F] ip3 = 1j * mp.pi / 3 ep = es_exp([mp.mpc(0)] + [ip3] + [mp.mpc(0)] * (K - 2), K) em = es_exp([mp.mpc(0)] + [-ip3] + [mp.mpc(0)] * (K - 2), K) h = [(a - b) / 1j for a, b in zip(es_mul(ep, F, K), es_mul(em, Fc, K))] G1 = gamma1p_es(1, K) G2 = gamma1p_es(2, K) gg = es_mul(G2, es_inv(es_mul(G1, G1, K), K), K) # bracket_k = (3/2) h_{k+2} - pi*gg_{k+1} (k >= -2; poles must cancel) res_m2 = mp.mpf(3) / 2 * h[0] res_m1 = mp.mpf(3) / 2 * h[1] - mp.pi * gg[0] br = [mp.mpf(3) / 2 * h[k + 2] - mp.pi * gg[k + 1] for k in range(K - 2)] p3 = es_exp([mp.mpc(0)] + [-mp.log(3)] + [mp.mpc(0)] * (K - 2), K) B = es_mul(p3, br, K - 2) Bre = [mp.re(c) / 2 for c in B[:nB]] imax = max(abs(mp.im(c)) for c in B[:nB]) polres = max(abs(res_m2), abs(res_m1)) thresh = mp.mpf(10) ** (-(D - POLE_MARGIN)) if not (polres < thresh): raise RuntimeError( f"boundary_B: eps^(-2,-1) pole-cancellation residual " f"{mp.nstr(polres, 3)} >= 10^-(dps-{POLE_MARGIN}) = " f"{mp.nstr(thresh, 3)} at ambient dps {D} (must sit at the " f"roundoff floor ~1e-{D}; healthy headroom 1e{POLE_MARGIN})") thresh_im = mp.mpf(10) ** (-(D - BIM_MARGIN)) if not (imax < thresh_im): raise RuntimeError( f"boundary_B: Im residual {mp.nstr(imax, 3)} >= " f"10^-(dps-{BIM_MARGIN}) = {mp.nstr(thresh_im, 3)} at ambient " f"dps {D} (B tower must be real to roundoff)") cert = {"B_tail_bound": bound_B, "B_tail_tol": mp.mpf(10) ** (-tol_dps)} return Bre, polres, imax, cert # --------------------------------------------------------------------------- # Full master assembly # --------------------------------------------------------------------------- def parse_pt(w): if isinstance(w, str) and "/" in w: n, d = w.split("/") return mp.mpf(n.strip()) / mp.mpf(d.strip()) return mp.mpf(w) def master_package(w, dps, B=None, eps_orders=3): """All named ingredients + assembled sector-15 master at Euclidean w < 0. Returns dict; work is done at dps+15 with series sized to dps+10.""" with mp.workdps(dps + 15): t = parse_pt(w) if not (t < 0): raise ValueError("Euclidean domain only: need real w < 0 (m^2=1); " "singular fibres at w in {0,1,9}; Minkowski w " "needs an i0+ prescription not implemented here") wd = dps + 10 psi1s, f2s, f3s, f4s, tser, q, qcert = sized_series(t, wd) N = len(tser) - 1 L = mp.log(abs(q)) + 1j * mp.pi # = 2 pi i tau_C, Euclidean branch tau = L / (2j * mp.pi) psi1_q = horner(psi1s, q) # psi1/pi, q-series route # PF routes (independent codepath), with certified accumulated bounds vp0, yL, e_vp0, e_yL = pf_transport(t, wd) psi1_pf = 2 / mp.sqrt(3) * vp0 tau_pf = mp.mpf(1) / 2 - 1j * (yL / vp0 - mp.log(9)) / (2 * mp.pi) # value-level certified bounds via exact prefactors (quotient rule # with the certified denominators) e_psi1_pf = 2 / mp.sqrt(3) * e_vp0 e_tau_pf = ((e_yL + abs(yL / vp0) * e_vp0) / (abs(vp0) - e_vp0) / (2 * mp.pi)) psi2_q = tau * psi1_q # psi2/pi # word integrals W = WordTable(f2s, f3s, f4s, N) eich1 = W.value(("f3",), q, L) # reference-value convention, real eich = W.value(("one", "f3"), q, L) # two-fold, enters the assembly K = eps_orders bcert = None if B is None: B, polres, imres, bcert = boundary_B(3, tol_dps=dps + TAIL_GUARD_2F1) F1 = [mp.mpc(1), -W.value(("one",), q, L) / 2, W.value(("one", "f4"), q, L)][:K] F3 = [eich, W.value(("one", "f3", "f2"), q, L), W.value(("one", "f3", "f2", "f2"), q, L) + W.value(("one", "f4", "one", "f3"), q, L)][:K] G1 = gamma1p_es(1, K) P = es_mul(es_mul(G1, G1, K), es_exp([mp.mpc(0)] + [-W.value(("f2",), q, L)] + [mp.mpc(0)] * (K - 2), K), K) Bs = [mp.mpc(b) for b in B[:K]] R = [a + b for a, b in zip(es_mul(F1, Bs, K), F3)] S = es_mul(P, R, K) J = [-psi1_q * c for c in S] # AMFlow convention J = -S111 # RAISING check: the assembled J # must be real up to roundoff on the Euclidean branch. Measured # healthy residual ~10^-(dps+15) (ambient roundoff), threshold # 10^-(dps-JIM_MARGIN) relative => ~1e23 headroom. jim = max(abs(mp.im(c)) for c in J) jscale = max(max(abs(mp.re(c)) for c in J), mp.mpf(1)) jthresh = jscale * mp.mpf(10) ** (-(dps - JIM_MARGIN)) if not (jim < jthresh): raise RuntimeError( f"master_package: assembled-J Im residual {mp.nstr(jim, 3)} >= " f"10^-(dps-{JIM_MARGIN}) x scale = {mp.nstr(jthresh, 3)} at " f"w = {mp.nstr(t, 10)}, dps {dps} (Euclidean-branch/assembly " f"inconsistency)") cert = {"q_tail_log10": qcert["q_tail_log10"], "q_tail_tol_log10": qcert["q_tail_tol_log10"], "N": qcert["N"], "nome_shift": qcert["nome_shift"], "nome_shift_tol": qcert["nome_shift_tol"], "e_psi1_pf": e_psi1_pf, "e_tau_pf": e_tau_pf, "B_tail_bound": (bcert or {}).get("B_tail_bound"), "B_tail_tol": (bcert or {}).get("B_tail_tol")} return {"w": t, "q": q, "tau": tau, "tau_pf": tau_pf, "psi1_over_pi_q": psi1_q, "psi1_over_pi_pf": psi1_pf, "psi2_over_pi": psi2_q, "eichler_f3": mp.re(eich1), "I_1f3": mp.re(eich), "J": J, "J_im_residual": jim, "B": B, "cert": cert} def assemble_master(w, c1, c2, dps): """General named form of the elliptic block (eps^0 row, in pi=stripped period units): g(w) = (psi1/pi)(w) * [Eichler_f3(tau(w)) + c1] + c2 * (psi2/pi)(w). General w-pencil variation-of-parameters solution. For the box-dressed sectors 55a-e/63 no per-master (c1,c2) exist as posed (pointwise named form structurally excluded -- moving-fibre one-fold; the dressing closes via the staged spectral assembly). Sector 15 is (c1,c2) = (B^(0), 0) with an overall minus (J = -S).""" pk = master_package(w, dps) return pk["psi1_over_pi_q"] * (pk["eichler_f3"] + c1) + c2 * pk["psi2_over_pi"] # --------------------------------------------------------------------------- # Comparisons G1-G6 (reference values read through the sha256 pins) # --------------------------------------------------------------------------- def _pinned_bytes(path, pin): """Read a companion file and refuse it unless its sha256 matches the pin.""" name = os.path.basename(path) if not os.path.exists(path): print(f"MISSING: {name} must sit beside this script (download it from the " f"same page); nothing computed") sys.exit(EXIT_MISSING) raw = open(path, "rb").read() sha = hashlib.sha256(raw).hexdigest() if sha != pin: print(f"REFUSED: {name} sha256 {sha} ({len(raw)} bytes) does not match the " f"pin {pin} -- the file was altered or is not the released version; " f"nothing computed") sys.exit(EXIT_PIN) return raw def load_data(): """row27_data.json through its sha256 pin (refused on mismatch, exit 3).""" return json.loads(_pinned_bytes(DATA_FILE, DATA_SHA256)) def load_gate_amflow(): """row27_gate_amflow.json through its sha256 pin (refused on mismatch, exit 3).""" return json.loads(_pinned_bytes(GATE_FILE, GATE_SHA256)) def agree_d(a, b): d = abs(a - b) if d == 0: return float(mp.mp.dps) m = max(abs(a), abs(b)) return float(-mp.log10(d / m)) def _pnum(s): """Parse a recorded midpoint string 'x' or 'x + y im'.""" if " + " in s and s.endswith("im"): re_s, im_s = s.split(" + ") return mp.mpc(mp.mpf(re_s), mp.mpf(im_s[:-2])) return mp.mpf(s) def _ball_mid(s): """Midpoint of an Arb-style ball string '[x +/- r]' (or a plain number).""" s = str(s).strip() if s.startswith("["): return mp.mpf(s[1:-1].split("+/-")[0].strip()) return mp.mpf(s) def gate_independent_amflow(g): """[G6] independent AMFlow values vs the locked closed-form prediction, both RECORDED in row27_gate_amflow.json (the prediction was written to a sha256-stamped file before any AMFlow output existed). The digits of agreement are computed here from the recorded strings; neither side is recomputed. Compared: w = -7/3, -4 only (kept out of the construction: no fit or anchor set contains them); w = -1 is a construction reference point -> control only.""" orders = ["-2", "-1", "0", "1", "2"] lock = g["locked_prediction"] print("\n[G6] INDEPENDENT AMFlow values (blackbox numeric IBP, uncut equal-mass " "sunrise, -S_111(4-2eps,w))\n vs the LOCKED closed-form prediction " f"(lock sha256 {lock['sha256'][:16]}..., locked {lock['locked_at']},\n" " before any AMFlow output; both sides recorded in " "row27_gate_amflow.json, digits of agreement computed here):") ok = True heldout_min, g60_heldout_min = mp.inf, mp.inf with mp.workdps(200): for goal, label in (("goal60", "goal 60"), ("goal40", "goal 40 (2nd prec)")): for pt in ("-1", "-7/3", "-4"): ds = [agree_d(_ball_mid(g["amflow_raw_outputs"][goal][pt][k]["re"]), mp.mpf(g["locked_values_stringD150"][pt][k])) for k in orders] dmin = min(ds) role = ("control: construction reference point (+ convention ID), " "EXCLUDED" if pt == "-1" else "KEPT-OUT comparison point") print(f" {label:18s} w={pt:>4}: " + " ".join(f"eps^{o} {d:6.1f}" for o, d in zip(orders, ds)) + f" min {dmin:5.1f} d [{role}]") if pt != "-1": heldout_min = min(heldout_min, dmin) if goal == "goal60": g60_heldout_min = min(g60_heldout_min, dmin) ok &= dmin > 30 rec = g["gate_digits"]["headline_min_over_orders_heldout"] consistent = abs(float(g60_heldout_min) - float(rec)) < 0.15 ok &= consistent print(f" kept-out min: {float(g60_heldout_min):.1f} d (goal 60) / " f"{float(heldout_min):.1f} d (both precisions); bar 30 d; " f"recorded headline {rec} d " f"[{'consistent' if consistent else 'MISMATCH'}]") return ok, float(g60_heldout_min) def run_gates(dps, data, quick=False): ok = True certs = [] with mp.workdps(dps + 40): B, polres, imres, bcert = boundary_B(5, tol_dps=dps + TAIL_GUARD_2F1) print("[G5] boundary tower from the printed 2F1 eps-form " "(Pfaff series, finite-difference-free):") with mp.workdps(dps + 40): cl2 = 2 * mp.im(mp.polylog(2, mp.exp(1j * mp.pi / 3))) d_cl = agree_d(B[0], cl2) print(f" B^(0) vs classical 2*Cl2(pi/3): {d_cl:7.1f} d " f"(recorded PSLQ residual 110.0 d)") for k in range(5): d_b = agree_d(B[k], mp.mpf(data["B"][str(k)])) print(f" B^({k}) vs 120-digit reference value: {d_b:7.1f} d") ok &= d_b > min(30, dps - 12) print(f" eps^(-2,-1) pole cancellation residual: {mp.nstr(polres, 3)} " f"(must be ~0); max Im residual {mp.nstr(imres, 3)}") ok &= float(polres) < 10 ** (-(dps + 20)) print(f" [floor] reference-precision detection floor: the B^(k) " f"reference values are 120-digit strings -- this table cannot " f"detect a discrepancy below ~1e-120 (agreement saturates at " f"min(ambient {dps + 40}, 120) d); honest floor of the " f"comparison, not a certificate beyond it") pts = ["-1", "-5"] if quick else ["-1", "-3", "-5"] print("\n[G4] FULL MASTER J[1,1,1](d0=2) = -S111(2-2eps,w), eps^0..2, vs the " "certified reference values\n (301-digit literature-representation " "strings, comparison only; HOMOGENEOUS-PERIOD-FRAME reference -- supplemented by the") print(" independent AMFlow leg [G6] below):") worst = mp.inf for p in pts: pk = master_package(p, dps, B=B) certs.append(pk["cert"]) orc = data["oracle_J111_d2"][p] ds = [] with mp.workdps(dps + 15): for k in range(3): d_k = agree_d(mp.re(pk["J"][k]), mp.mpf(orc[str(k)]["re"])) ds.append(d_k) worst = min(worst, d_k) print(f" w={p:>3}: eps^0 {ds[0]:6.1f} d eps^1 {ds[1]:6.1f} d " f"eps^2 {ds[2]:6.1f} d (cap min(dps={dps}, 301); " f"Im residual {mp.nstr(pk['J_im_residual'], 3)})") ok &= min(ds) > min(30, dps - 12) print(f" [floor] reference-precision detection floor: the reference " f"strings carry 301 digits -- this comparison cannot detect a " f"discrepancy below ~1e-301 (agreement caps at min(dps+15 = " f"{dps + 15}, 301) d); honest floor of the comparison, not " f"a certificate beyond it") g6ok, g6min = gate_independent_amflow(load_gate_amflow()) ok &= g6ok print("\n[G1-G3] period / tau / Eichler cross-checks at the six reference points:") apts = ["-1/4", "-2"] if quick else ["-1/4", "-1/2", "-1", "-2", "-3", "-5"] for p in apts: pk = master_package(p, dps, B=B) certs.append(pk["cert"]) anch = data["anchors"][p] with mp.workdps(dps + 15): d_route = agree_d(pk["psi1_over_pi_q"], pk["psi1_over_pi_pf"]) d_bank = agree_d(pk["psi1_over_pi_q"], mp.mpf(anch["psi1_over_pi"])) d_tau = agree_d(pk["tau"], pk["tau_pf"]) d_psi2 = agree_d(pk["psi2_over_pi"], _pnum(anch["psi2_over_pi"])) d_ei = agree_d(pk["eichler_f3"], mp.mpf(data["eichler_f3"][p])) print(f" w={p:>4}: psi1 PF-vs-qseries {d_route:6.1f} d | psi1 vs ref " f"{d_bank:5.1f} d | tau 2nd-sol-vs-nome {d_tau:6.1f} d | " f"psi2 vs ref {d_psi2:5.1f} d | Eichler vs ref {d_ei:5.1f} d") ok &= d_route > min(30, dps - 12) and d_tau > min(30, dps - 12) print(" (reference-value comparisons cap at ~80 d -- the stored strings; " "the psi2 reference at w=-5 is itself") print(" limited by the reference run's NQ=500 q-convergence (the same " "family as its recorded 35.7 d floor);") print(f" closure record: min " f"{data['recorded_gate_digits']['homog_cross_min_d']} d homogeneous cross-check / " f"{data['recorded_gate_digits']['ell_top_loo_min_d']} d elliptic-top leave-one-out " f"at PREC=400)") # certified-bound summary (fail-closed -- every line below is # backed by a RAISE inside the corresponding routine) wq = max(c["q_tail_log10"] for c in certs) wsh = max(c["nome_shift"] for c in certs) wpf = max(c["e_psi1_pf"] for c in certs) wtau = max(c["e_tau_pf"] for c in certs) print("\n[certified] fail-closed runtime truncation bounds (worst over " "comparison points; RAISE on non-convergence):") print(f" q-series tail (PROVEN eta-majorant) <= 1e{wq:.1f} " f"(tol 1e{certs[0]['q_tail_tol_log10']}, N = " f"{max(c['N'] for c in certs)}); nome-shift <= {mp.nstr(wsh, 3)} " f"(tol {mp.nstr(certs[0]['nome_shift_tol'], 3)})") print(f" PF-transport accumulated: psi1 <= {mp.nstr(wpf, 3)}, tau <= " f"{mp.nstr(wtau, 3)} (per-step tol 1e-{dps + 10 + TAIL_GUARD_STEP}); " f"2F1 tower B-tail <= {mp.nstr(bcert['B_tail_bound'], 3)} " f"(tol {mp.nstr(bcert['B_tail_tol'], 3)})") print(f" raising checks armed: 2F1 pole/Im residuals (10^-(ambient-" f"{POLE_MARGIN})), assembled-J Im (10^-(dps-{JIM_MARGIN}) rel), " f"Newton non-convergence, --check agreement (dps-8)") return ok, float(worst), g6min def _mutated_boundary_B(exact): """--mutate control: a wrapper over boundary_B that multiplies c1^(0) = B^(0) by (1 + 1e-9) before it enters the assembly (the eps^0 row is -(psi1/pi)[B^(0) + Eichler_f3]). Installed from main() only; every numeric function is untouched. G4 and G5 must then fail and the run must exit 1. G6 compares recorded values and is unaffected by construction.""" def mutated(nB, tol_dps=None): Bre, polres, imres, cert = exact(nB, tol_dps) Bre[0] = Bre[0] * (1 + mp.mpf("1e-9")) return Bre, polres, imres, cert return mutated def main(): ap = argparse.ArgumentParser( description="Row 27 (LBL3E) full elliptic master, standalone (mpmath). " "Domain: real Euclidean w < 0, m^2 = 1.") ap.add_argument("--dps", type=int, default=80, help="target digits (default 80; minimum 30)") ap.add_argument("--point", action="append", default=None, help="Euclidean w < 0 (repeatable; use the --point=-7/3 form for negatives)") ap.add_argument("--check", action="store_true", help="two-precision rule: rerun at dps+60 and diff (default mode: the " "comparison table at dps+60 with fewer points; --point mode: the " "assembled J)") ap.add_argument("--quick", action="store_true", help="fewer comparison points") ap.add_argument("--mutate", action="store_true", help="control: multiply c1^(0) = B^(0) by (1 + 1e-9) before the assembly; " "the run must exit nonzero") args = ap.parse_args() if args.dps < 30: ap.error("--dps must be >= 30 (the comparison bar is 30 digits)") # argparse exits 2 # fail-closed: both data files are read through their sha256 pins BEFORE any computation data = load_data() gate = load_gate_amflow() # run_gates re-reads (and re-checks) it as well print(f"[pins] row27_data.json sha256 {DATA_SHA256[:16]}... and row27_gate_amflow.json " f"sha256 {GATE_SHA256[:16]}... match ({len(gate['locked_values_stringD150'])} " f"recorded AMFlow points)") if args.mutate: globals()["boundary_B"] = _mutated_boundary_B(boundary_B) print("MUTATE CONTROL: c1^(0) = B^(0) is multiplied by (1 + 1e-9) before the " "assembly; G4 and G5 must FAIL and this run must exit nonzero") mp.mp.dps = args.dps + 15 t0 = time.perf_counter() rc = EXIT_PASS try: if args.point: for p in args.point: try: pk = master_package(p, args.dps) except ValueError as ex: print(f"w = {p}: DOMAIN REFUSAL -- {ex}") rc = EXIT_USAGE continue print(f"w = {p} (dps {args.dps}; |q_C| = {float(abs(pk['q'])):.4f})") print(f" psi1/pi (q-series) = {mp.nstr(pk['psi1_over_pi_q'], args.dps)}") print(f" psi1/pi (PF route) = {mp.nstr(pk['psi1_over_pi_pf'], args.dps)}" f" [cross: {agree_d(pk['psi1_over_pi_q'], pk['psi1_over_pi_pf']):.1f} d]") print(f" tau_C = {mp.nstr(pk['tau'], 30)}" f" [2nd-solution route: {agree_d(pk['tau'], pk['tau_pf']):.1f} d]") print(f" psi2/pi = {mp.nstr(pk['psi2_over_pi'], 30)}") print(f" Eichler_f3(tau(w)) = {mp.nstr(pk['eichler_f3'], args.dps)}" " [one-fold I(f3;q) convention of the reference values]") print(f" I(1,f3;q) = {mp.nstr(pk['I_1f3'], args.dps)}" " [two-fold; enters the eps^0 assembly]") for k in range(3): print(f" J^({k})(w) = {mp.nstr(mp.re(pk['J'][k]), args.dps)}") print(f" (Im residual {mp.nstr(pk['J_im_residual'], 3)}; sector-15 " "representative master, (c1,c2) = (B^(0),0); the box-dressed " "sectors have no pointwise (c1,c2) -- see assemble_master(w,c1,c2,dps))") ct = pk["cert"] print(f" [certified] q-series tail <= 1e{ct['q_tail_log10']:.1f} " f"(PROVEN eta-majorant, N = {ct['N']}, tol " f"1e{ct['q_tail_tol_log10']}); nome-shift <= " f"{mp.nstr(ct['nome_shift'], 3)}") print(f" [certified] PF-transport: psi1 <= " f"{mp.nstr(ct['e_psi1_pf'], 3)}, tau <= " f"{mp.nstr(ct['e_tau_pf'], 3)}; 2F1 tower B-tail <= " f"{mp.nstr(ct['B_tail_bound'], 3)} (all fail-closed: RAISE " f"on non-convergence)") if args.check: pk2 = master_package(p, args.dps + 60) d = min(agree_d(mp.re(pk['J'][k]), mp.re(pk2['J'][k])) for k in range(3)) print(f" --check: dps+60 rerun agrees to {d:.1f} d (>= {args.dps - 8} required)") if d < args.dps - 8: raise RuntimeError( f"--check FAILED at w = {p}: dps+60 rerun agrees to " f"only {d:.1f} d < required {args.dps - 8} d " f"(dps {args.dps} vs {args.dps + 60})") else: ok, worst, g6min = run_gates(args.dps, data, quick=args.quick) gate_fail = not ok print(f"\nGATE SUMMARY: {'PASS' if ok else 'FAIL'} " f"(worst full-master reference agreement {worst:.1f} d at dps {args.dps}; " f"independent AMFlow leg kept-out min {g6min:.1f} d)") if args.check: print("\n--check: rerunning the full-master comparison at dps+60 ...") ok2, worst2, g6min2 = run_gates(args.dps + 60, data, quick=True) gate_fail = gate_fail or not ok2 print(f"--check summary: {'PASS' if ok2 else 'FAIL'} " f"(worst {worst2:.1f} d at dps {args.dps + 60}; digits must GROW; " f"independent leg min {g6min2:.1f} d, dps-independent recorded strings)") if gate_fail: rc = EXIT_FAIL except RuntimeError as ex: # fail-closed: a certified bound not reached, or a --check disagreement print(f"\nFAIL (exit {EXIT_FAIL}): {ex}") rc = EXIT_FAIL print(f"\nmeasured wall time: {time.perf_counter() - t0:.2f} s") if rc == EXIT_PASS and not args.point: print("RESULT: PASS (every comparison above its bar; exit 0)") elif rc == EXIT_FAIL: print("RESULT: FAIL (exit 1)") return rc if __name__ == "__main__": sys.exit(main())