#!/usr/bin/env python3 r"""sunrise-row12-evaluate.py -- two-loop sunrise, unequal masses, eps^1 correction: STANDALONE evaluator (row 12 of the paper's table of integrals). Two-loop sunrise S_111(2-2eps, t), one external scale t = p^2 (Euclidean t < 0), three massive propagators with squared masses (m1^2, m2^2, m3^2) in {(1,1,2), (1,2,3), (1,1,4)}, mu = 1. Convention J_111 = -S_111 (AW sign); J^(n) = eps^n Taylor coefficient. CLOSED FORM (zero free parameters; agrees with independent AMFlow reference values to 288-308 digits at --full --dps 300, reference-capped at every point, for all three masses -- the data file marks none of them as held out of a fit): E^(1)(t) := J^(1)(t) / psihat1(t) = -J_4^(3)(t) + b * J_4^(2)(t), b = 2*gamma_E + 2*log(m3^2) (exact frame constant, NOT a fit) * psihat1 = |psi1_F|/pi, the Adams-Weinzierl Feynman-curve holomorphic period, computed here from the BMSW curve roots (elliptic K). * J_4^(2) = the weight-2 eps^0 master (the object of sunrise-row09/10/11-evaluate.py): C_{4,2} plus the depth-2 Kronecker-eMPL words I(omega_0, omega_3(z_j, N tau_C)); E^(0) = -J_4^(2) (positive control in the gate table). * J_4^(3) = the Bogner-Mueller-Stach-Weinzierl [BMSW, arXiv:1907.01251] weight-3 depth<=3 Kronecker-eMPL polynomial: 101 terms (6 depth-1, 8 depth-2, 86 depth-3, plus the boundary constant C_{4,3}) over a 14-letter kernel alphabet {omega_0, omega_2(z_j, N tau_C), omega_3(z_j, N tau_C), eta_2 : j=1..3, N=1,2}, 13 words carrying the second-kind kernel eta_2 = (e_2(tau)-2e_2(2tau)) dtau/(2 pi i). * C_{4,3} boundary constant in CLOSED FORM via the Nielsen reduction Li_{2,1}(x,1) = S_{1,2}(x) = -Li_3(1-x) + ln(1-x) Li_2(1-x) + (1/2) ln x ln^2(1-x) + zeta(3) (Koelbig 1986) -- exact on |x|=1, where the marked-point exponentials live for Euclidean t. (The naive nested sum converges like log(n)/n^2 there; an under-converged C_{4,3} was the source of an earlier uniform ~1e-6 residual.) RUNTIME MACHINERY (all in this file, mpmath only): 1. BMSW Feynman curve E_{3,F}: roots e_iF(t), modulus k_F^2, tau_F, period psi1_F, and the three Abel-Jacobi marked points z_1,z_2,z_3 (z1+z2+z3 = 1) -- elliptic integrals. 2. Kernel nome q_C = -exp(i pi tau_F) = exp(2 pi i tau_C), tau_C = (tau_F+1)/2; negative real in the documented domain. 3. Tau-iterated-integral engine: each one-form pulled back to C(qbar) dqbar/qbar, iterated integrals built inner-first with the shuffle / tangential-base-point regularisation at the cusp (formal letter L = log qbar carried as a polynomial tower; L^0 part kept). Engine validations: shuffle product to 1e-64, independent kernel rederivation to 99 digits, depth-3 shuffle with divergent letters in every slot to 52 digits. 4. Boundary blocks in closed form: L1, L2, L_Delta (marked-point logs), C_{4,2} (elliptic dilogarithm), C_{4,3} (Nielsen S_{1,2} + Li_3 + zeta(3)). DATA FILE sunrise-row12-data.json (same directory; sha256-pinned below, REFUSED on mismatch): * J43_terms / J42_terms: EXACT rational coefficients + word lists, parsed from the BMSW supplementary material (supplementary_material.mpl of arXiv:1907.01251). Published closed-form symbolic data, not a numeric cache. * gates: independent AMFlow reference values (d = 2-2eps, 288-308 digits) of the (1,1,1) master at eps^0 and eps^1 -- comparison only, never in the computation; no point is marked as held out of any fit (every one is a prediction of the zero-parameter closure); each point carries the provenance of its AMFlow value (goal digits, eps order, threads, both ball radii, the sha256 of the input configuration and of the output file, wall, landed stamp; the four points re-run at goal 280 also carry their goal-250 twin: output sha256 and the pair agreement). * paper_claims: the digit counts the paper states for this row -- 288-308 over the nine points at --full --dps 300; per mass (1,1,2) 288/288, (1,2,3) 288/308, (1,1,4) 288/288 (worst/best) -- gated at that tier (see below). * paper_literals: the digits the paper prints for E^(1) (trailing digits truncated). STATED GAPS (honest scope): * The 101-word J_4^(3) polynomial is INHERITED from the published BMSW supplementary material (arXiv:1907.01251), its 101 coefficients re-fit from the reference grid alone, every coefficient free, with 4 held-out points at >= 104 digits (sunrise-row12-j43-refit.json beside this file, sha-pinned), identical to the published values; the full 14-letter alphabet is not closed from scratch (see that file); the contribution of this work at eps^1 is the zero-parameter closure E^(1) = -J_4^(3) + b J_4^(2) (sign, frame constant b in closed form, corrected C_{4,3} evaluation) validated to 288-308 digits (--full --dps 300) against those independent reference values. * Gate reference values are AMFlow numerics (declared above); the evaluation itself never touches them. DOMAIN: Euclidean t < 0 (below the soft point t=0 and all thresholds; q_C negative real, marked points real). Arbitrary rational/decimal t < 0 via --point. |q_C| >= 0.8995 or non-Euclidean t is refused explicitly (--point prints REFUSED and exits 2). Physical t > 0 is NOT wired (analytic continuation of the q-series frame not vendored). RUNTIME CERTIFICATION (every evaluation, fail-closed): * FRAME-AGREEMENT DOUBLING GATE: the seed depth Nq is only a starting point; the full word polynomial is re-evaluated on a second frame at 2*Nq (exact continuation of the same q-series recurrences) and the E^(0) and E^(1) two-depth agreements must beat tol = 10^-(working dps - 2) relative (8 guard digits past the reported dps), or Nq doubles up to a 64x-seed cap and then the run REFUSES to print a value. * CLOSED-FORM BLOCK TRIPWIRE: every evaluate_point (i) re-derives b via the digamma route -2*psi(0,1) + 4*log(m3), (ii) recomputes C_{4,2}/C_{4,3} via INDEPENDENT quadrature references (Euler integral representations of Li_2, Li_3 and Nielsen S_{1,2}, tanh-sinh at working dps + 10), and (iii) recomputes every used block at +10 dps (L1/L2 recompute-only -- no independent route wired for the plain logs; disclosed). Any mismatch > 10^-(working dps - 6) relative raises BEFORE any value is printed. A 1e-30 perturbation of b or of C_{4,3} therefore cannot surface as a silent digit drop in the gate table. * TWO-PRECISION SELF-CHECK (--point): the two certifications above run at ONE working precision, so a loss of working precision inside the evaluation (rounding in the curve or block arithmetic at that dps) passes both. Every --point evaluation is therefore repeated at dps + 60, both values are printed with their agreement, and the run FAILS by name (exit 1) when they agree to fewer than dps - 12 digits. Specimen: (1,1,4) at t = -2 -- the value at dps 60 is good to about 35 digits while the frame-agreement gate and the block tripwire both PASS; the self-check names it (at dps 150 the pair still agrees to fewer than the bar, at dps 250 it clears). Not a table point: the nine table points are read directly against their references. WHAT THE DEFAULT RUN CHECKS (bare invocation; seconds): (1,1,2) at t = -5, dps 30: E^(1) vs the 288-digit independent AMFlow reference (bar: >= 30 digits, capped by min(reference digits, dps + 10)), E^(0) as positive control, and the digits the paper prints, E^(1) = 12.166814... (8 digits, trailing digits truncated), matched to within one unit in the last printed place. --full: the whole 3 masses x 3 points table at dps 60 (minutes rather than seconds). --full --dps 300: the nine-point table at the paper's counts (288-308 digits, every point read to its reference's own digit count; measured 15:42.95 wall, 186752 kB; longer on a loaded host) -- the tier at which the paper_claims lines are gated (PASS iff the measured worst >= the stated worst - 1, per mass and over the nine; at any lower tier they print SKIP by name with that run's cap). --dps 120 (the quick point at dps 120): 129 digits against the 288-digit reference at (1,1,2) t = -5 (cap 130) at a longer wall, which the script prints; the paper's own counts need the --full --dps 300 tier above. EXIT CODES: 0 every gate passes; 1 a gate fails (this is what --mutate must produce), a fail-closed certification refusal, a paper-claim FAIL at its tier or a --point self-check FAIL; 2 usage error or --point domain refusal; 3 the data file's or the re-fit file's sha256 does not match its pin (refused before any computation); 4 the data file, the re-fit file or mpmath is missing. USAGE python3 sunrise-row12-evaluate.py # quick default (see above) python3 sunrise-row12-evaluate.py --full # 3 masses x 3 points, dps 60 python3 sunrise-row12-evaluate.py --full --dps 300 # the paper's row-12 counts (minutes): # 288-308 digits, paper_claims gated python3 sunrise-row12-evaluate.py --dps 120 # quick point at dps 120: 129 digits # (slower) python3 sunrise-row12-evaluate.py --full --mass 123 # one mass, all its points python3 sunrise-row12-evaluate.py --mass 123 --point=-7/2 # arbitrary Euclidean t<0; carries the # two-precision self-check (dps, dps+60) # (use --point= for negatives) python3 sunrise-row12-evaluate.py --full --check # rerun the table at dps+60, diff python3 sunrise-row12-evaluate.py --mutate # control: first J_4^(3) coefficient # 6 -> 6(1+1e-12); MUST exit nonzero mp.dps is set inside main() after argparse (module-level mpf footgun avoided). Dependency: python3 + mpmath (pip install mpmath). """ import argparse import hashlib import json import os import re import sys import time from fractions import Fraction 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) DATA_FILE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "sunrise-row12-data.json") DATA_SHA256 = "478027348a582a925c379aaf80e43958fc56c62a5e8df94f8b0a14264c1b48f6" REFIT_FILE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "sunrise-row12-j43-refit.json") REFIT_SHA256 = "b7a8fbc6a3ce87e2edf7f880821f237108e3e99528f4dc289847a215a3c51b54" def _pinned_bytes(path, pin): """Read the data 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 same page); " f"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 pin " f"{pin} -- the file was altered or is not the released version; nothing computed") sys.exit(EXIT_PIN) return raw # --------------------------------------------------------------------------- # helpers (no module-level mpf constants!) # --------------------------------------------------------------------------- def _tpi(): return 2 * mp.pi * mp.mpc(0, 1) def _csqrt(x): return mp.sqrt(mp.mpc(x)) def _rat(s): """Exact rational string -> mpf at CURRENT dps.""" fr = Fraction(s) return mp.mpf(fr.numerator) / mp.mpf(fr.denominator) # --------------------------------------------------------------------------- # 1. BMSW Feynman curve E_{3,F} (mass-generic; port of frameF.fcurve) # --------------------------------------------------------------------------- def fcurve(t, m1, m2, m3, mu=1): """Curve data at t: roots e1F,e2F,e3F, kF2, kFp2, tauF, psi1F, marked points z1..z3. m1,m2,m3 are MASSES (not squared); BMSW def_roots / points_on_E_hat / def_coordinate_torus.""" t = mp.mpc(t) mu = mp.mpf(mu) m1s, m2s, m3s = m1 * m1, m2 * m2, m3 * m3 M100 = m1s + m2s + m3s mu1 = -m1 + m2 + m3 mu2 = m1 - m2 + m3 mu3 = m1 + m2 - m3 mu4 = m1 + m2 + m3 Delta = mu1 * mu2 * mu3 * mu4 mu4p = mu ** 4 rad = 3 * (_csqrt(mu1 * mu1 - t) * _csqrt(mu2 * mu2 - t) * _csqrt(mu3 * mu3 - t) * _csqrt(mu4 * mu4 - t)) e1F = (-t * t + 2 * M100 * t + Delta + rad) / (24 * mu4p) e2F = (-t * t + 2 * M100 * t + Delta - rad) / (24 * mu4p) e3F = (2 * t * t - 4 * M100 * t - 2 * Delta) / (24 * mu4p) Z1F = e3F - e2F Z2F = e1F - e3F Z3F = e1F - e2F kF2 = Z1F / Z3F kFp2 = -Z1F / Z2F tauF = mp.mpc(0, 1) * mp.ellipk(1 - kF2) / mp.ellipk(kF2) psi1F = 2 / _csqrt(Z3F) * mp.ellipk(kF2) Kp = mp.ellipk(kFp2) def xhat(mi2, mj2): return e3F + mi2 * mj2 / mu4p def zF(xjk): up = _csqrt((e1F - e3F) / (xjk - e3F)) return mp.ellipf(mp.asin(up), kFp2) / (2 * Kp) z1F = zF(xhat(m2s, m3s)) z2F = zF(xhat(m3s, m1s)) z3F = zF(xhat(m1s, m2s)) # Abel-Jacobi BRANCH REPAIR: the lattice relation z1+z2+z3 = 1 must hold exactly # (AGW: z3 = 1 - z1 - z2). The ellipf/asin principal branch reflects z3 -> 1 - z3 # once z3 crosses the half-period (e.g. (1,1,4) for -2 < t < 0: raw sum = 0.903 at # t=-1); earlier evaluation grids simply avoided that window. Enforcing the relation # repairs the branch: validated here against the independent AMFlow reference at # (1,1,4) t=-1 (eps^0 AND eps^1 close to full 60-digit working precision, alongside # the unaffected t=-3,-5 points). z1, z2 stay far from 1/2 on the documented domain. ssum = z1F + z2F + z3F if abs(ssum - 1) > mp.mpf(10) ** (-mp.mp.dps // 2): z3F = 1 - z1F - z2F return dict(e1F=e1F, e2F=e2F, e3F=e3F, kF2=kF2, kFp2=kFp2, tauF=tauF, psi1F=psi1F, z1F=z1F, z2F=z2F, z3F=z3F) def psihat1_of(cv): return abs(cv["psi1F"]) / mp.pi def nome_qC_of(cv): return -mp.e ** (mp.mpc(0, 1) * mp.pi * cv["tauF"]) # --------------------------------------------------------------------------- # 2. AGW tau-iterated-integral engine (port of itint.py: Frame/LSeries/F_word) # --------------------------------------------------------------------------- class Frame: """qbar-series of the kernels g^(2)/g^(3)(z_j, N tau_C), b_2, at truncation Nq (Kronecker-Eisenstein q-series, Broedel-Duhr-Dulat-Tancredi eq. 46 grammar).""" def __init__(self, qC, z1, z2, z3, Nq): self.qC = mp.mpc(qC) self.z = {1: mp.mpc(z1), 2: mp.mpc(z2), 3: mp.mpc(z3)} self.Nq = int(Nq) self._g = {} self._b2 = None self._form = {} self._suffix = {} def _gn_qseries(self, n, z, N): """g^(n)(z, N tau) as qbar-power coefficients (BMSW eq.46 grammar).""" Nq = self.Nq w = mp.exp(_tpi() * z) a = [mp.mpc(0)] * (Nq + 1) pref = -(_tpi()) ** n / mp.factorial(n - 1) a[0] += pref * (-mp.bernoulli(n) / n) sgn = (-1) ** (1 - n) j = 1 while N * j <= Nq: wj = w ** j - sgn * w ** (-j) k = 1 while N * j * k <= Nq: a[N * j * k] += pref * wj * (mp.mpf(k) ** (n - 1)) k += 1 j += 1 return a def gn(self, n, z_idx, N): key = (n, z_idx, N) if key not in self._g: self._g[key] = self._gn_qseries(n, self.z[z_idx], N) return self._g[key] def b2(self): """b_2(tau) = e_2(tau) - 2 e_2(2 tau): coeff of qbar^n = sum of odd divisors.""" if self._b2 is None: Nq = self.Nq c = [mp.mpc(0)] * (Nq + 1) pre = 2 * (_tpi()) ** 2 c[0] = pre * (mp.mpf(1) / 24) for nn in range(1, Nq + 1): s = 0 d = 1 while d <= nn: if nn % d == 0 and d % 2 == 1: s += d d += 1 c[nn] = pre * s self._b2 = c return self._b2 def _qmul(a, b, Nq): """Truncated Cauchy product of two qbar-series. The nonzero entries of b are listed once (mpc comparisons dominate an entry-by-entry loop); the products and their accumulation order are exactly those of the plain double loop.""" out = [mp.mpc(0)] * (Nq + 1) nzb = [(j, bj) for j, bj in enumerate(b) if j <= Nq and bj != 0] for i, ai in enumerate(a): if i > Nq: break if ai == 0: continue for j, bj in nzb: if i + j > Nq: break out[i + j] += ai * bj return out class LSeries: """Polynomial in the formal letter L = log(qbar) with qbar-power-series coefficients.""" def __init__(self, terms, Nq): self.Nq = Nq self.terms = terms @classmethod def one(cls, Nq): q = [mp.mpc(0)] * (Nq + 1) q[0] = mp.mpc(1) return cls({0: q}, Nq) def mul_qseries(self, C): return LSeries({p: _qmul(q, C, self.Nq) for p, q in self.terms.items()}, self.Nq) def regint(self): """Shuffle / tangential-base-point regularised primitive under dqbar/qbar. n>=1 monomials: exact antiderivative (L-tower mixes DOWN via falling factorials); n==0: L raised by one. Validated to 50-52 digits incl. depth-3 shuffle with divergent letters in every slot.""" Nq = self.Nq out = {} def add(pp, qlist): base = out.get(pp) if base is None: out[pp] = list(qlist) else: for i in range(Nq + 1): base[i] += qlist[i] for p, s in self.terms.items(): a = [mp.mpf(1)] for m in range(1, p + 1): a.append(a[-1] * (-(p - m + 1))) for m in range(0, p + 1): pm = p - m q = [mp.mpc(0)] * (Nq + 1) am = a[m] hit = False for n in range(1, Nq + 1): sn = s[n] if sn != 0: q[n] = sn * am / (mp.mpf(n) ** (m + 1)) hit = True if hit: add(pm, q) if s[0] != 0: qr = [mp.mpc(0)] * (Nq + 1) qr[0] = s[0] / (p + 1) add(p + 1, qr) return LSeries(out, Nq) def value_L0(self, qC): s = self.terms.get(0) if s is None: return mp.mpc(0) val = mp.mpc(0) qn = mp.mpc(1) for n in range(0, self.Nq + 1): val += s[n] * qn qn *= qC return val def series_for_form(name, fr): """C(qbar) for an AGW form: C = (dtau-coefficient)/(2 pi i). Cached per frame.""" if name in fr._form: return fr._form[name] Nq = fr.Nq tpi = _tpi() if name == "omega_0_tau": C = [mp.mpc(0)] * (Nq + 1) C[0] = mp.mpc(1) elif name == "eta_2_tau": C = [x / (tpi * tpi) for x in fr.b2()] else: m = re.match(r"omega_(\d)_z(\d)_(2?)tau$", name) if not m: raise ValueError(f"unknown form {name!r}") k = int(m.group(1)) zi = int(m.group(2)) N = 2 if m.group(3) == "2" else 1 if k == 2: C = [N * x / (tpi * tpi) for x in fr.gn(2, zi, N)] elif k == 3: pre = (2 * mp.pi) ** (-1) * 2 * N / (tpi * tpi) C = [pre * x for x in fr.gn(3, zi, N)] else: raise ValueError(f"unsupported omega_{k}") fr._form[name] = C return C def F_word(forms, fr): """F(omega_1,...,omega_k), inner-first recursion S <- regint(C*S), with a per-frame SUFFIX cache (words sharing inner subwords reuse the partial LSeries).""" key = tuple(forms) nf = len(forms) S = None start = 0 for k in range(nf, 0, -1): suf = tuple(forms[nf - k:]) if suf in fr._suffix: S = fr._suffix[suf] start = k break if S is None: S = LSeries.one(fr.Nq) for i in range(nf - start - 1, -1, -1): C = series_for_form(forms[i], fr) S = S.mul_qseries(C).regint() fr._suffix[tuple(forms[i:])] = S return S.value_L0(fr.qC) # --------------------------------------------------------------------------- # 3. Closed-form boundary blocks (wbar_j = exp(2 pi i z_j), on |x|=1 for t<0) # --------------------------------------------------------------------------- def wbars(cv): return [mp.exp(_tpi() * cv[k]) for k in ("z1F", "z2F", "z3F")] def L1_L2_LD(cv): w1, w2, _w3 = wbars(cv) L1 = mp.log(w2) + 2 * mp.log(1 - w1) - 2 * mp.log(1 - w1 * w2) L2 = mp.log(w1) + 2 * mp.log(1 - w2) - 2 * mp.log(1 - w1 * w2) LD = mp.log(-(1 - w1) ** 2 * (1 - w2) ** 2 / (1 - w1 * w2) ** 2) return L1, L2, LD def Li21(x): """Li_{2,1}(x,1) = Nielsen S_{1,2}(x), Koelbig closed reduction -- exact on |x|=1. zeta(3) enters the boundary constant HERE, in closed form.""" x = mp.mpc(x) return (-mp.polylog(3, 1 - x) + mp.log(1 - x) * mp.polylog(2, 1 - x) + mp.mpf(1) / 2 * mp.log(x) * mp.log(1 - x) ** 2 + mp.zeta(3)) def C42_closed(cv): return sum((mp.polylog(2, wj) - mp.polylog(2, 1 / wj)) / (2 * mp.mpc(0, 1)) for wj in wbars(cv)) def C43_closed(cv): """BMSW eq.1645-1652 boundary constant (with the corrected Nielsen Li21).""" s = mp.mpc(0) for wj in wbars(cv): s += (-2 * Li21(wj) - mp.polylog(3, wj) + 2 * Li21(1 / wj) + mp.polylog(3, 1 / wj)) / (2 * mp.mpc(0, 1)) _, _, LD = L1_L2_LD(cv) return s - LD * C42_closed(cv) # --------------------------------------------------------------------------- # 4. Assembly # --------------------------------------------------------------------------- def _qmax(): """QMAX = 0.8995 (exact rational 1799/2000), the same series-domain limit as the eps^0 shared layer sunrise_empl.py (duplicated here, NOT imported: this script stays standalone). Rationale: the q-series tail carries a 1/(1-|q_C|) factor whose supremum on |q_C| < QMAX is 1/(1-0.8995) = 9.9503 < 10 -- strictly one decimal guard digit. Zero-arg function, not a module-level mpf (import-dps footgun).""" return mp.mpf(1799) / 2000 def check_domain(t, qC): """Refuse points outside the documented domain instead of silently truncating a near-divergent q-series.""" if not (mp.im(mp.mpc(t)) == 0 and mp.re(mp.mpc(t)) < 0): raise ValueError(f"t = {t}: only Euclidean t < 0 is supported by this script") if abs(qC) >= _qmax(): raise ValueError(f"t = {t}: |q_C| = {mp.nstr(abs(qC), 6)} >= {mp.nstr(_qmax(), 4)} " "(series-domain limit; t too close to 0^- or too deep)") FRAME_GUARD_DROP = 2 # tol = 10^-(mp.mp.dps - 2); working dps = reported + 10 # => the gate certifies 8 guard digits past the REPORTED dps FRAME_NCAP_MULT = 64 # doubling cap: refuse rather than evaluate beyond 64 * seed BLOCK_TOL_DROP = 6 # block-tripwire tol = 10^-(wdps-6) rel = 10^-(dps+4): # detection floor strictly below a 1e-30 perturbation # for every dps >= 30 (never AT the 1e-30 edge) def _b_ref(msq): """INDEPENDENT re-derivation of the frame constant b = 2*gamma_E + 2*log(m3^2): gamma_E via the digamma route -psi(0,1), the log via 4*log(m3). Used only by _blocks_selfcheck (the block tripwire), never in the value path.""" return -2 * mp.psi(0, mp.mpf(1)) + 4 * mp.log(mp.sqrt(mp.mpf(msq[2]))) def _li2_quad(x): """Li_2(x) = -int_0^1 log(1-x u) du/u (Euler integral; log1p avoids the 0/0 node at u->0). Independent of mp.polylog.""" x = mp.mpc(x) return -mp.quad(lambda u: mp.log1p(-x * u) / u, [0, 1]) def _li3_quad(x): """Li_3(x) = (1/2) int_0^1 x log^2(u)/(1-x u) du. Independent of mp.polylog.""" x = mp.mpc(x) return mp.quad(lambda u: x * mp.log(u) ** 2 / (1 - x * u), [0, 1]) / 2 def _li21_quad(x): """Nielsen S_{1,2}(x) = Li_{2,1}(x,1) = (1/2) int_0^1 log^2(1-x u) du/u. Independent of the Koelbig closed reduction used by Li21().""" x = mp.mpc(x) return mp.quad(lambda u: mp.log1p(-x * u) ** 2 / u, [0, 1]) / 2 def _blocks_selfcheck(t, msq, cv, blocks, b): """Closed-form block tripwire: the USED closed-form constants must match their independent references or the evaluation refuses LOUDLY (nonzero exit) -- a 1e-30 perturbation of b or C_{4,3} must never surface as a silent digit drop in the gate table. Seven compares, all at tol = 10^-(wdps - BLOCK_TOL_DROP) relative: * b vs the digamma re-derivation _b_ref (independent route); * C_4_2, C_4_3 vs tanh-sinh quadrature refs at wdps+10 (independent integral representations of Li_2 / Li_3 / S_{1,2}); * all four used blocks + b vs +10-dps recomputes of the same closed forms (catches assembly-level perturbation; for L1/L2 -- plain logs -- this recompute is the only check wired: DISCLOSED, no independent route). Both routes consume the same curve data cv (treated as exact input), so the compares are independent of cv's own precision. Returns (max rel diff, tol) for the certified print lines.""" wdps = mp.mp.dps tol = mp.mpf(10) ** (-(wdps - BLOCK_TOL_DROP)) with mp.workdps(wdps + 10): bh = _b_ref(msq) L1h, L2h, LDh = L1_L2_LD(cv) C42h = C42_closed(cv) C43h = C43_closed(cv) ws = wbars(cv) ii2 = 2 * mp.mpc(0, 1) C42r = sum((_li2_quad(w) - _li2_quad(1 / w)) / ii2 for w in ws) s = mp.mpc(0) for w in ws: s += (-2 * _li21_quad(w) - _li3_quad(w) + 2 * _li21_quad(1 / w) + _li3_quad(1 / w)) / ii2 C43r = s - LDh * C42r checks = ( ("b used vs digamma ref", b, bh), ("C_4_2 used vs hi recompute", blocks["C_4_2"], C42h), ("C_4_3 used vs hi recompute", blocks["C_4_3"], C43h), ("L1 used vs hi recompute", blocks["L1"], L1h), ("L2 used vs hi recompute", blocks["L2"], L2h), ("C_4_2 closed vs quadrature ref", C42h, C42r), ("C_4_3 closed vs quadrature ref", C43h, C43r), ) worst = mp.mpf(0) bad = [] for name, used, ref in checks: d = abs(mp.mpc(used) - mp.mpc(ref)) / max(mp.mpf(1), abs(mp.mpc(ref))) if d > worst: worst = d if d > tol: bad.append(f"{name}: rel diff {mp.nstr(d, 3)} > tol {mp.nstr(tol, 3)}") if bad: raise RuntimeError( f"row12 closed-form block tripwire FAILED at t={t}, " f"masses^2={tuple(msq)}: " + "; ".join(bad) + " -- frame constant / boundary block does not match its " "independent reference; refusing to print a value") return worst, tol def seed_Nq(qC, safety=0): """STARTING seed only (speed, never answer): the frame-agreement doubling gate in evaluate_point is what certifies the depth at runtime.""" aq = abs(qC) return max(int((mp.mp.dps + safety) / max(-mp.log10(aq), mp.mpf("0.05"))) + 18, 90) def eval_poly(terms, fr, cv, blocks): """Evaluate a parsed word polynomial; blocks = dict of closed-form symbols.""" tot = mp.mpc(0) for coeff_s, sym, forms in terms: c = _rat(coeff_s) base = blocks[sym] if sym is not None else mp.mpc(1) if forms: base = base * F_word(forms, fr) tot += c * base return tot def evaluate_point(t, msq, terms43, terms42): """Full row-12 evaluation at Euclidean t<0 for squared masses msq=(m1^2,m2^2,m3^2). Returns dict with E0, E1, psihat1, J0, J1, Nq + certified frame-agreement data. FRAME-AGREEMENT DOUBLING GATE: seed_Nq is only the starting depth. The full word polynomial is re-evaluated on a second frame at 2*Nq -- an EXACT continuation of the same q-series recurrences (the first Nq coefficients of every kernel/convolution/primitive are bit-identical by construction: identical loops, identical summation order), so the difference IS the truncation-tail estimate (the two-successive-depth agreement). Accept the Nq values only if the E^(0) AND E^(1) agreements beat tol = 10^-(mp.mp.dps - FRAME_GUARD_DROP) relative; otherwise Nq doubles (reusing the 2*Nq evaluation as the new low frame) up to FRAME_NCAP_MULT * seed, then RuntimeError naming t/mass/|q_C|/seed/Nq/bound/tol/cap -- fail closed, never print a value the loop did not certify.""" m1, m2, m3 = (mp.sqrt(mp.mpf(m)) for m in msq) cv = fcurve(t, m1, m2, m3) qC = nome_qC_of(cv) check_domain(t, qC) L1, L2, LD = L1_L2_LD(cv) blocks = {"C_4_2": C42_closed(cv), "C_4_3": C43_closed(cv), "L1": L1, "L2": L2} b = 2 * mp.euler + 2 * mp.log(mp.mpf(msq[2])) # block tripwire: raises (exit 1) BEFORE any value is printed if the # used b / boundary blocks drift from their independent references blk_diff, blk_tol = _blocks_selfcheck(t, msq, cv, blocks, b) def _polys_at(nq): fr = Frame(qC, cv["z1F"], cv["z2F"], cv["z3F"], nq) j42n = eval_poly(terms42, fr, cv, blocks).real j43n = eval_poly(terms43, fr, cv, blocks).real return j42n, j43n n0 = seed_Nq(qC) ncap = FRAME_NCAP_MULT * n0 tol = mp.mpf(10) ** (-(mp.mp.dps - FRAME_GUARD_DROP)) Nq = n0 j42, j43 = _polys_at(Nq) while True: j42h, j43h = _polys_at(2 * Nq) # value-level error via l1-norm of the EXACT prefactors: # E0 = -j42 -> err_E0 = |d42| # E1 = -j43 + b * j42 -> err_E1 = |d43| + |b| |d42| errE0 = abs(j42h - j42) errE1 = abs(j43h - j43) + abs(b) * abs(j42h - j42) sc0 = max(mp.mpf(1), abs(j42)) sc1 = max(mp.mpf(1), abs(-j43 + b * j42)) if errE0 <= tol * sc0 and errE1 <= tol * sc1: break if 2 * Nq >= ncap: raise RuntimeError( f"row12 frame-agreement gate NOT converged at cap: t={t}, " f"masses^2={tuple(msq)}, |q_C|={mp.nstr(abs(qC), 6)}, seed Nq={n0}, " f"pair Nq={Nq} vs {2 * Nq}, cap={ncap}: E0 agreement " f"{mp.nstr(errE0, 3)} (tol {mp.nstr(tol * sc0, 3)}), E1 agreement " f"{mp.nstr(errE1, 3)} (tol {mp.nstr(tol * sc1, 3)}) -- refusing to " "print an uncertified value") Nq *= 2 j42, j43 = j42h, j43h E0 = -j42 E1 = -j43 + b * j42 ph1 = psihat1_of(cv) return {"E0": E0, "E1": E1, "psihat1": ph1, "J0": ph1 * E0, "J1": ph1 * E1, "b": b, "Nq": Nq, "Nq_seed": n0, "Nq_hi": 2 * Nq, "errE0": errE0, "errE1": errE1, "tol": tol, "errJ0": ph1 * errE0, "errJ1": ph1 * errE1, "blk_diff": blk_diff, "blk_tol": blk_tol} # --------------------------------------------------------------------------- # 5. Gates # --------------------------------------------------------------------------- def digits_agree(pred, oracle): rel = abs(pred - oracle) / max(abs(oracle), mp.mpf(1)) if rel == 0: return mp.mp.dps return int(-mp.log10(rel)) def run_gates(data, masses, terms43, terms42, only_t=None): """Gate table over the independent reference points of the given masses (only_t: restrict to these t strings, e.g. ["-5"] for the quick run).""" rows = [] for mass in masses: msq = data["masses_sq"][mass] for g in data["gates"][mass]: if only_t is not None and g["t"] not in only_t: continue t = _rat(g["t"]) t0 = time.time() r = evaluate_point(t, msq, terms43, terms42) wall = time.time() - t0 ora_J0 = mp.mpf(g["J0"]) ora_J1 = mp.mpf(g["J1"]) E0_ora = ora_J0 / r["psihat1"] E1_ora = ora_J1 / r["psihat1"] d0 = digits_agree(r["E0"], E0_ora) d1 = digits_agree(r["E1"], E1_ora) cap = min(g["digits_avail"], mp.mp.dps) rows.append({"mass": mass, "t": g["t"], "d0": d0, "d1": d1, "cap": cap, "E1": r["E1"], "Nq": r["Nq"], "wall": wall}) print(f" ({','.join(str(x) for x in msq)}) t={g['t']:>4}: " f"eps^0 control {d0:>3}d eps^1 GATE {d1:>3}d " f"(cap {cap}, Nq={r['Nq']}, {wall:.1f}s)") print(f" certified: frame agreement Nq={r['Nq']} vs {r['Nq_hi']}: " f"|dE0| = {mp.nstr(r['errE0'], 3)}, |dE1| = {mp.nstr(r['errE1'], 3)} " f"<= tol {mp.nstr(r['tol'], 3)} (rel)") print(f" certified: blocks vs independent refs (digamma b, " f"quadrature C42/C43, hi recomputes): max rel diff " f"{mp.nstr(r['blk_diff'], 3)} <= tol {mp.nstr(r['blk_tol'], 3)}") return rows def main(): ap = argparse.ArgumentParser( description="Row 12: unequal-mass sunrise eps^1, " "E^(1) = -J_4^(3) + (2 gamma_E + 2 log m3^2) J_4^(2)") ap.add_argument("--dps", type=int, default=None, help="reported digits (default 30 for the quick run, 60 with --full or --point; " "+10 guard digits internally)") ap.add_argument("--mass", choices=["112", "123", "114", "all"], default=None, help="squared-mass assignment m1^2 m2^2 m3^2 (default: 112 for the quick run, " "all with --full)") ap.add_argument("--point", type=str, default=None, help="Euclidean t < 0 (rational like -7/2 or decimal; use --point= for " "negatives); mass from --mass (default 112); carries a two-precision " "self-check (the value at dps and at dps+60 printed; FAIL by name below " "dps - 12 agreeing digits)") ap.add_argument("--full", action="store_true", help="all independent reference points (3 masses x 3 points unless --mass) at dps 60; " "with --dps 300 the paper's stated digit counts (paper_claims) are gated, at a " "lower dps they print SKIP by name") ap.add_argument("--check", action="store_true", help="two-precision rule: rerun the table at dps+60 and diff (--point carries this " "check always)") ap.add_argument("--mutate", action="store_true", help="control: perturb the first J_4^(3) coefficient 6 -> 6(1+1e-12); " "the run must exit nonzero") args = ap.parse_args() dps = args.dps if args.dps is not None else (60 if (args.full or args.point is not None) else 30) if dps < 30: ap.error("--dps must be >= 30 (the gate bar is 30 digits)") # argparse exits 2 raw = _pinned_bytes(DATA_FILE, DATA_SHA256) # refused before any computation print(f"[pin] {os.path.basename(DATA_FILE)} sha256 {DATA_SHA256[:16]}... matches") _pinned_bytes(REFIT_FILE, REFIT_SHA256) # the re-fit of the J_4^(3) coefficients: pinned, never read by the value path print(f"[pin] {os.path.basename(REFIT_FILE)} sha256 {REFIT_SHA256[:16]}... matches") mp.mp.dps = dps + 10 # guard digits; reported values at dps t_start = time.time() data = json.loads(raw) terms43 = data["J43_terms"] terms42 = data["J42_terms"] if args.mutate: if terms43[0][0] != "6": print("REFUSED: mutate control expects the first J_4^(3) coefficient to be 6") return EXIT_PIN terms43[0][0] = "6000000000001/1000000000000" # 6 (1 + 1e-12), still an exact rational print("[mutate] CONTROL RUN: first J_4^(3) coefficient 6 -> 6 (1 + 1e-12); " "this run MUST fail the eps^1 gate and exit nonzero") print(f"row 12: unequal sunrise eps^1 | J_4^(3): {len(terms43)} terms " f"(86 depth-3, 13 with eta_2, 14-letter alphabet) | dps={dps}") if args.point is not None: mass = args.mass if args.mass not in (None, "all") else "112" msq = data["masses_sq"][mass] t = _rat(args.point) if not (t < 0): print(f"\npoint: t={args.point}: REFUSED (only Euclidean t < 0 is wired)") return EXIT_USAGE def one(d): mp.mp.dps = d + 10 return evaluate_point(t, msq, terms43, terms42) try: r = one(dps) except ValueError as e: print(f"\npoint: masses^2=({','.join(str(x) for x in msq)}) " f"t={args.point}: REFUSED ({e})") return EXIT_USAGE print(f"\npoint: masses^2=({','.join(str(x) for x in msq)}) t={args.point} (Nq={r['Nq']})") print(f" b = 2*gamma_E + 2*log(m3^2) = {mp.nstr(r['b'], min(dps, 40))}") for k in ("E0", "E1", "psihat1", "J0", "J1"): print(f" {k:>8} = {mp.nstr(r[k], dps)}") print(f" certified: frame agreement Nq={r['Nq']} (seed {r['Nq_seed']}) vs " f"{r['Nq_hi']}: |dE0| = {mp.nstr(r['errE0'], 3)}, " f"|dE1| = {mp.nstr(r['errE1'], 3)}, |dJ1| = {mp.nstr(r['errJ1'], 3)} " f"<= tol {mp.nstr(r['tol'], 3)} (rel)") print(f" certified: blocks vs independent refs (digamma b, quadrature " f"C42/C43, hi recomputes): max rel diff {mp.nstr(r['blk_diff'], 3)} " f"<= tol {mp.nstr(r['blk_tol'], 3)}") # TWO-PRECISION SELF-CHECK (always in --point mode; --check adds nothing here): the # frame-agreement gate and the block tripwire above certify the q-series depth and the # closed-form blocks at ONE working precision, so a loss of working precision inside # the evaluation (rounding in the curve / block arithmetic at that dps) passes both. # The point is re-evaluated at dps + 60, both values are printed, and the run FAILS # by name when the two agree to fewer than dps - 12 digits. Specimen: (1,1,4) at # t = -2 at dps 60 (about 35 digits). fails = [] r2 = one(dps + 60) print(f"\n two-precision self-check: the point re-evaluated at dps {dps + 60} " f"(bar: the two values agree to >= {dps - 12} digits, else FAIL)") for k in ("E0", "E1"): d = min(digits_agree(r[k], r2[k]), dps) ok = d >= dps - 12 print(f" {k} at dps {dps}: {mp.nstr(r[k], dps)}") print(f" {k} at dps {dps + 60}: {mp.nstr(r2[k], dps + 60)}") print(f" {k}: the two agree to {d} digits (bar {dps - 12}) [{'ok' if ok else 'FAIL'}]") if not ok: fails.append(f"two-precision self-check {k}: the values at dps {dps} and dps {dps + 60} " f"agree to {d} digits < {dps - 12} (the printed value is not good to " f"dps {dps}; rerun with a higher --dps)") print(f"\nwall time: {time.time() - t_start:.1f}s (measured)") if fails: print(f"[gate] FAIL (exit {EXIT_FAIL}): " + "; ".join(fails)) return EXIT_FAIL print("[gate] PASS: value certified (frame-agreement doubling + block tripwire) and the " "two-precision self-check stable") return EXIT_PASS full = args.full if args.mass == "all" or (args.mass is None and full): masses = ["112", "123", "114"] else: masses = [args.mass or "112"] only_t = None if full else ["-5"] print("\nGATE TABLE (reference: independent AMFlow values; eps^0 = positive control = the " "row-9/10/11 object)" + ("" if full else " [quick run: t = -5 only]")) rows = run_gates(data, masses, terms43, terms42, only_t) fails = [] # the digits the paper PRINTS (trailing digits truncated): the computed value must # agree to within one unit in the last printed place for lit in data.get("paper_literals", []): hit = [r for r in rows if r["mass"] == lit["mass"] and r["t"] == lit["t"]] if not hit: continue sref = lit["E1"] ref = mp.mpf(sref) nprint = len(sref.replace("-", "").replace(".", "").lstrip("0")) ulp = mp.mpf(10) ** (int(mp.floor(mp.log10(abs(ref)))) - (nprint - 1)) diff = abs(hit[0]["E1"] - ref) ok = diff < ulp if not ok: fails.append(f"paper literal E^(1) ({lit['mass']}) t={lit['t']}: printed {sref}..., " f"computed {mp.nstr(hit[0]['E1'], nprint + 4)}, |diff| = {mp.nstr(diff, 2)} " f"not < {mp.nstr(ulp, 1)}") print(f" paper literal E^(1) ({','.join(str(x) for x in data['masses_sq'][lit['mass']])}) " f"t={lit['t']}: printed {sref}... ({nprint} digits), computed " f"{mp.nstr(hit[0]['E1'], nprint + 4)}...: |diff| = {mp.nstr(diff, 2)} < one unit in " f"the last printed place ({mp.nstr(ulp, 1)}) [{'PASS' if ok else 'FAIL'}] " f"[{lit['where']}]") worst0 = min(min(r["d0"], r["cap"]) for r in rows) worst1 = min(min(r["d1"], r["cap"]) for r in rows) npts = len(rows) print(f"\n MEASURED: eps^1 gate worst {worst1} digits over {npts} point" f"{'s' if npts != 1 else ''} (eps^0 control worst {worst0}d); PASS(>=30): {worst1 >= 30}") if worst1 < 30: fails.append(f"eps^1 gate worst {worst1} digits < 30") if worst0 < 30: fails.append(f"eps^0 positive control worst {worst0} digits < 30") # paper claims (the data file's paper_claims; row-9 form): the digit counts the paper # states for this row, measured here as the worst / best eps^1 agreement (each point # capped by its reference's own digits) per mass and over the nine points. Testable # only at the claim's tier (--full over the three masses at --dps >= the tier's 300, # where every cap is the reference's own count): PASS iff the measured worst >= the # stated worst - 1 (one digit of slack: the reference's own last digit), per mass and # over the nine; at a lower tier SKIP by name with this run's cap. A claim over a # subset this script does not measure is a data error (FAIL), never silently passed. claim_tested = False for cl in data.get("paper_claims", []): if cl.get("measured_over") != "all_nine_reference_points": fails.append(f"paper claim '{cl['text']}': measured_over = {cl.get('measured_over')!r} is " f"not 'all_nine_reference_points' (the only subset this script measures)") continue mt = re.fullmatch(r"--full --dps (\d+)", cl["tier"]) if not mt: fails.append(f"paper claim '{cl['text']}': tier {cl['tier']!r} is not of the form " f"'--full --dps N'") continue tier_dps = int(mt.group(1)) cap_run = max(r["cap"] for r in rows) if not full or dps < tier_dps: why = ((f"quick run, {npts} point{'s' if npts != 1 else ''}" if not full else "the table") + f" at dps {dps}" + ("" if dps >= tier_dps else f" < the tier's {tier_dps}")) print(f" paper claim '{cl['text']}' [{cl['where']}]: not testable in this run ({why}; " f"cap {cap_run} = min(reference digits, dps + 10); the tier is {cl['tier']}) [SKIP]") continue for mass in masses: mrows = [r for r in rows if r["mass"] == mass] w = min(min(r["d1"], r["cap"]) for r in mrows) b = max(min(r["d1"], r["cap"]) for r in mrows) cw, cb = int(cl["per_mass"][mass]["worst"]), int(cl["per_mass"][mass]["best"]) ok = w >= cw - 1 ms = ",".join(str(x) for x in data["masses_sq"][mass]) print(f" paper claim per mass ({ms}): eps^1 worst {w} / best {b} digits over {len(mrows)} " f"points (stated {cw} / {cb}; bar worst >= {cw - 1}) [{'PASS' if ok else 'FAIL'}]") if not ok: fails.append(f"paper claim per mass ({ms}): worst {w}d < {cw - 1}d (stated {cw})") if len(masses) == 3: claim_tested = True nw, nb = int(cl["worst_digits"]), int(cl["best_digits"]) best1 = max(min(r["d1"], r["cap"]) for r in rows) ok = worst1 >= nw - 1 print(f" paper claim '{cl['text']}' [{cl['where']}] at its tier {cl['tier']}: eps^1 worst " f"{worst1} / best {best1} digits over {npts} points (stated {nw} / {nb}; bar worst >= " f"{nw - 1}) [{'PASS' if ok else 'FAIL'}]") if not ok: fails.append(f"paper claim '{cl['text']}': worst {worst1}d over {npts} points < {nw - 1}d " f"(stated {nw})") else: print(f" paper claim '{cl['text']}' [{cl['where']}]: not testable in this run (mass " f"{', '.join(masses)} only; the claim is over the nine points of the three masses; " f"the tier is {cl['tier']}) [SKIP]") if args.check: print(f"\n--check: rerunning the table at dps={dps + 60}") mp.mp.dps = dps + 70 rows2 = run_gates(data, masses, terms43, terms42, only_t) worst_diff = mp.mp.dps for r, r2 in zip(rows, rows2): worst_diff = min(worst_diff, digits_agree(r["E1"], r2["E1"])) stab = min(worst_diff, dps) print(f" --check: E^(1) values stable across dps to >= {stab} digits " f"[{'ok' if stab >= dps - 12 else 'DRIFT'}]") if stab < dps - 12: fails.append(f"--check two-precision drift: {stab} < {dps - 12}") if any(min(r2["d1"], r2["cap"]) < 30 for r2 in rows2): fails.append("--check rerun: an eps^1 gate row fell below 30 digits") print(f"\nwall time: {time.time() - t_start:.1f}s (measured)") if not full: print("[quick run] rerun with --full for the 3 masses x 3 reference points table at dps 60") if fails: print(f"\n[gate] FAIL (exit {EXIT_FAIL}): " + "; ".join(fails)) return EXIT_FAIL print(f"\n[gate] PASS: eps^1 independent-reference gate (worst {worst1} digits, bar 30), eps^0 " f"positive control" + (", the paper's printed digits" if data.get("paper_literals") and any(r["mass"] == "112" and r["t"] == "-5" for r in rows) else "") + (", the paper's stated digit counts" if claim_tested else "") + (" and --check stability" if args.check else "") + " all clear") return EXIT_PASS if __name__ == "__main__": try: sys.exit(main()) except RuntimeError as e: # frame-agreement cap or block tripwire: fail-closed refusal print(f"\n[gate] FAIL (exit {EXIT_FAIL}): {e}") sys.exit(EXIT_FAIL)