#!/usr/bin/env python3 """Equal-mass two-loop sunrise diagram J111(2-2eps, t), m^2 = mu^2 = 1 — evaluate the eps^0 closed form AT RUNTIME from the kinematic input t and compare against a held-out certified oracle. Standalone: python3 + mpmath. Closed form (this work): J^(0)(t) = -(psi1/pi) * ( (3/2)*sqrt(3)*L(chi_{-3},2) + I(1,f3; q_C) ) Modular ingredients on Gamma_1(6), conventions of arXiv:1704.08895: e1(q) = 1/6 + sum_{m>=1} (sum_{d|m} chi_{-3}(d)) q^m, e2(q) = e1(q^2) psi1/pi = 2*sqrt(3)*(e1 + e2) f3 = 36*sqrt(3)*(e1^3 - e1^2 e2 - 4 e1 e2^2 + 4 e2^3) = sum_{n>=1} a_n q^n t(q) = 9 q prod_n (1-q^{6n})^8 (1-q^n)^4 (1-q^{2n})^{-8} (1-q^{3n})^{-4} I(1,f3;q)= sum_{n>=1} a_n / n^2 q^n Everything above is built from scratch here: the q-series are generated at runtime, the hauptmodul t(q) is inverted by Newton iteration to get q_C(t), and the Eichler sum is evaluated term by term. Kinematic branches: * Euclidean t < 0: q_C is real negative (cusp-connected branch), Newton seeded at q = t/9. * Physical t > 9 with +i0: q_C is complex. It is obtained by ANALYTIC CONTINUATION at runtime — Newton-tracking q along a t-path from the Euclidean anchor t = -3 through the upper half t-plane (arc height 3, which stays clear of the singular points t = 0, 1, 9), ending at real t on the +i0 side. This mirrors the certified oracle's own two-path branch tracking (winding k = 0). Oracle (held out, embedded as literals below so the script is standalone): first 128 digits of the certified 1024-bit ball-arithmetic evaluation of the all-orders Gamma_1(6) literature representation (arXiv:1704.08895), fetched only after the blind bootstrap result was committed. Repo artifact: stage1/reveal/out/reference_sunrise_merged.tsv (d0=2, eps_power=0 rows) independently confirmed by held-out AMFlow runs (e.g. stage1/reveal/heldout/out/sunrise_d2_tm3_g150_o11b.json, ~160 certified digits at t=-3). The t = -9 entry is the eps^0 midpoint of a further AMFlow run (goal 150, eps_order 11, 2026-09-06; ~160 certified digits, its goal-100 twin agreeing to ~110) made only after this script's value at t = -9 had been filed: a sixth never-fit point beside the withheld t = -3 run. The agreement printed below is RECOMPUTED at runtime as -log10(|computed - oracle| / |oracle|); it is capped by mp.dps and the 128-digit oracle excerpts, not by the method. Archived (offline, NOT computed by this script): the full certified reveal gate compared blind-vs-oracle in 1024-bit ball arithmetic and found 304-305 digit agreement across all 50 ball-overlap tests (10 points x 5 eps-orders). Interface (evaluate at YOUR point and precision, no code edits needed): python3 sunrise-evaluate.py # gate demo + dps-doubling demo python3 sunrise-evaluate.py --point -7.25 --dps 200 # arbitrary point/precision python3 sunrise-evaluate.py --point 12 --point -5 --dps 220 or from Python (file name has a hyphen, so load via importlib): f(t, dps) returns J^(0)(t) as mpf (Euclidean) / mpc (physical). Domain: t = p^2/m^2 with t < 0 (Euclidean; for t < -3 the nome is tracked along the real t-axis from the anchor t = -3, which crosses no singular points) or t > 9 (physical, +i0 via the complex-arc continuation above). 0 <= t <= 9 (singular points t = 0, 1, 9 and the between-thresholds strip) is NOT covered by this script. Precision is limited only by dps: the q-series length is auto-sized from a MEASURED truncation-tail bound (top built coefficient x |q_C|^N; the eta-quotient coefficients grow like 10^(~1.5 sqrt(n)), so this is checked, not assumed) — more digits just means more series terms (no stored numeric table caps it). Points with |q_C| >= 0.92 (t -> 0^-, -inf, or 9^+) are refused with an explicit message rather than silently truncated. Dependency: pip mpmath. """ import mpmath as mp mp.mp.dps = 110 NTERMS = 700 # default q-series length; f() auto-extends it per (t, dps) QMAX = 0.92 # refuse |q_C| beyond this (series-domain limit, see docstring) # ---- axis-3 hardening knobs (2026-07-06; all seeds/margins, no caps) -------- TAIL_GUARD = 8 # certified truncation tolerance = 10^-(dps+TAIL_GUARD) NUM_GUARD = 12 # raising gate on Newton/roundoff nome floor = 10^-(dps-NUM_GUARD) CERT_CAP = 8 # max x1.5 series extensions in the certified loop (fail-closed) GATE_MARGIN = 8 # gate table RAISES below min(dps, ORACLE_DIGITS) - GATE_MARGIN CRANK_STEP = 40 # in-run D -> D+CRANK_STEP two-precision crank gate (raising) CRANK_MARGIN = 6 # crank gate RAISES below dps - CRANK_MARGIN (calibrated 2026-07-06) # ---- held-out certified oracle, eps^0, first 128 significant digits -------- # Provenance: /Physics/Bootstrap/stage1/reveal/out/reference_sunrise_merged.tsv # (re_mid / im_mid columns, d0=2, eps_power=0; certified literature representation, # never used in the bootstrap fit). t=-3 also confirmed by an independent AMFlow # held-out run, stage1/reveal/heldout/out/sunrise_d2_tm3_g150_o11b.json. # t=-9: the AMFlow run of 2026-09-06 (goal 150, eps_order 11; output sha256 # 04901ddcac467744..., config 84feb6f28583e0a8..., goal-100 twin 4b75d285e8308d74...), # made after this script's value at t=-9 was filed; first 128 digits of its # eps^0 midpoint. REF = { -3: "-2.0589766979254918724182799691488449049958852132591171652425844676811576183594946409331691857618941814230305619245144390988781967", -5: "-1.9161259714466100706039079919221565571280358799665609058226481378073315747706936274128691068126757888869334454365408686839810674", -9: "-1.6979464611954261219419261948835053057039346015060125104535582050775611024877819046960812360785531969330402744676323938496644114", } REF12 = ("-2.3758041831615816443545900918692305823990797520896306758079022409394144293423681886048566410679938493335318929795821108366616888", "-4.5915809030764085059917091047911926064725476772484902124680632162033592129573825126580474150924580552355219596040181594985234149") 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: for j in range(0, N + 1 - i): if b[j]: r[i + j] += ai * b[j] return r def euler_sparse(N): """Euler function prod(1-x^n) by the pentagonal number theorem: [(exp, sign)].""" 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): # sp[0] == (0, 1) 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 series(N=NTERMS): """Coefficient arrays (length N+1): psi1/pi, a_n of f3, t(q).""" 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) 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))] 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]] # eta-quotient weight w=1 shift, x9 return psi1, f3, t _SERIES_CACHE = {} def get_series(N=NTERMS): """series(N), cached per (working precision, N).""" key = (mp.mp.dps, N) if key not in _SERIES_CACHE: _SERIES_CACHE[key] = series(N) return _SERIES_CACHE[key] def horner(c, q): acc = mp.mpf(0) for ck in reversed(c): acc = acc * q + ck return acc def _newton(q, ttarget, tser, dser, iters=60): tol = mp.mpf(10) ** (-mp.mp.dps + 5) for _ in range(iters): dq = (horner(tser, q) - ttarget) / horner(dser, q) q = q - dq if abs(dq) < tol: break return q def nome_from_t(tser, tval): """Euclidean branch (t < 0): real negative nome, Newton seeded at q = t/9.""" dser = [k * tser[k] for k in range(1, len(tser))] return _newton(mp.mpf(tval) / 9, mp.mpf(tval), tser, dser) def nome_physical(tser, tval, t_anchor=-3, height=3, steps=400): """Physical branch (t > 9, +i0): analytic continuation of the nome computed at runtime. Newton-track q along the t-path t(u) = t_anchor + (tval - t_anchor)*u + i*height*sin(pi*u), u: 0 -> 1, an arc through the upper half t-plane (+i0 prescription) clear of the singular points t = 0, 1, 9; endpoint polished at real t = tval.""" dser = [k * tser[k] for k in range(1, len(tser))] q = _newton(mp.mpf(t_anchor) / 9, mp.mpf(t_anchor), tser, dser) # real anchor for k in range(1, steps + 1): u = mp.mpf(k) / steps re = t_anchor + (tval - t_anchor) * u # height=0 -> real path (deep-Euclidean tracking, stays real mpf) tt = mp.mpc(re, height * mp.sin(mp.pi * u)) if height else mp.mpf(re) q = _newton(q, tt, tser, dser, iters=12) return _newton(q, mp.mpf(tval), tser, dser) # polish on the real axis def J0_from_q(q, psi1, f3): """eps^0 closed form evaluated at nome q (real or complex, |q| < 1).""" P = horner(psi1, q) eich = mp.fsum(f3[n] / mp.mpf(n) ** 2 * q ** n for n in range(1, len(f3))) L2 = (mp.zeta(2, mp.mpf(1) / 3) - mp.zeta(2, mp.mpf(2) / 3)) / 9 # L(chi_{-3},2) return -P * (mp.mpf(3) / 2 * mp.sqrt(3) * L2 + eich) def nome_for_t(tser, tval): """q_C(t) on the demo branches. For t < -3 the direct seed q = t/9 leaves the unit disk and Newton can land on a SPURIOUS root of the truncated series, so the nome is tracked along the real t-axis from the t = -3 anchor instead (no singular points on t < 0).""" if -3 <= tval < 0: return nome_from_t(tser, tval) if tval < -3: return nome_physical(tser, tval, height=0) # real-axis tracking if tval > 9: return nome_physical(tser, tval) raise ValueError("domain: Euclidean t < 0 or physical t > 9 (+i0); " "0 <= t <= 9 (singular points t = 0, 1, 9) is not covered") def J0(tval, psi1, f3, tser): """J^(0)(t) for Euclidean t < 0 (real branch) or physical t > 9 (+i0).""" return J0_from_q(nome_for_t(tser, tval), psi1, f3) # -------- certified truncation bounds (axis-3 hardening, 2026-07-06) -------- # Construction (2026-07-06): # (certified geometric tail bounds, refine-until-bound, fail-closed raises), # adapted from Taylor transports to the Gamma_1(6) q-series. Every quantity # below is a BOUND, not an estimate: # (i) |psi1_n| <= 3*sqrt(3)*n for n >= 1: |e1_n| <= sigma_0(n) <= n and # e2_n = e1_{n/2}, so |e1_n + e2_n| <= (3/2) n. # (ii) |f3_n| <= 1080*sqrt(3)*n^5: f3 = 36 sqrt3 (e1^3 - e1^2 e2 # - 4 e1 e2^2 + 4 e2^3); each triple convolution has at most # (n+1)(n+2)/2 <= 3 n^2 terms (n >= 1), every factor <= max(1/6, idx) # <= n, and the absolute coefficient sum is 1+1+4+4 = 10. # (iii) |t_n| <= 9 G(rho) rho / rho^n for any qa < rho < 1: the eta-quotient # coefficients are dominated termwise by the all-inverted majorant # G(x) = prod_k (1-x^k)^-4 (1-x^2k)^-8 (1-x^3k)^-4 (1-x^6k)^-8 # (|coeffs of (1-x^k)^m| = C(m,j) <= C(m+j-1,j) = coeffs of # (1-x^k)^-m), G has nonnegative coefficients, so [x^m] G <= # G(rho)/rho^m (Cauchy); t = 9 x * (shifted quotient). log G(rho) is # a finite sum plus the termwise remainder -log(1-x) <= x/(1-x). # (iv) nome error, Newton-Kantorovich: with m = T_t + r_res (T_t = t-series # tail at qa, r_res = measured Newton residual |t_trunc(q_C) - t|), # a' = |t'_trunc(q_C)| - T_t' (lower bound on |t'|), b = D2 (bound on # |t''| on |z| <= qb = qa + r0: triangle inequality on the computed # coefficients plus the majorant tail): if h = m*b/a'^2 <= 1/4 the true nome # lies within 2m/a' of q_C. The truncation part dq_trunc = 2 T_t / a' # feeds the certified budget; the Newton/roundoff part dq_num = # 2 r_res / a' is floored by the working precision (= dps, unchanged # for byte-identity of the recorded values) and is gated separately at # 10^-(dps-NUM_GUARD), plus end-to-end by the D/D+40 crank gate. # Total certified truncation budget (c0 = (3/2) sqrt3 L(chi-3,2); PA/EA/DJ = # absolute-series majorants: partial sums at qb plus their closed-form tails): # |J0 err| <= T_psi*(c0 + EA) + PA*T_eich + DJ*dq_trunc # f() RAISES (fail-closed) if this cannot be driven under 10^-(dps+TAIL_GUARD) # within CERT_CAP x1.5 series extensions. def _logG(rho): """Rigorous upper bound on log G(rho) (majorant product of (iii)).""" s, k, rk = mp.mpf(0), 1, rho while rk > mp.mpf(10) ** -30: s -= (4 * mp.log(1 - rk) + 8 * mp.log(1 - rk ** 2) + 4 * mp.log(1 - rk ** 3) + 8 * mp.log(1 - rk ** 6)) k += 1 rk = rho ** k # k > K remainder, termwise -log(1-x) <= x/(1-x) <= x/(1-rho^{K+1}): s += (4 * rk / (1 - rho) + 8 * rk ** 2 / (1 - rho ** 2) + 4 * rk ** 3 / (1 - rho ** 3) + 8 * rk ** 6 / (1 - rho ** 6)) / (1 - rk) return s def _t_scale_pick(qa, M, r0): """(T_t, rho, S_t) minimizing the certified t-series tail over a small rho grid; every grid member gives a valid bound, so the min is one.""" best = None for rho in (mp.sqrt(qa), qa ** (mp.mpf(1) / 3), qa ** mp.mpf("0.25"), mp.mpf("0.96")): if not (qa + r0 + mp.mpf("0.005") < rho < mp.mpf("0.995")): continue with mp.workdps(50): St = 9 * mp.exp(_logG(rho)) * rho # |t_n| <= St / rho^n y = qa / rho Tt = St * y ** (M + 1) / (1 - y) # sum_{n>M} |t_n| qa^n if best is None or Tt < best[0]: best = (Tt, rho, St) return best def _certify(psi1, f3, tser, q, tval): """Certified error budget at the evaluation point (see block comment). Returns a dict of BOUNDS; never raises (caller gates fail-closed).""" qa = abs(q) M = len(tser) - 1 s3 = mp.sqrt(3) r0 = mp.mpf("0.01") * (1 - qa) qb = qa + r0 # derivative bounds hold on |z| <= qb def tail_n1(C, x): # C * sum_{n>M} n x^n (exact identity) return C * x ** (M + 1) * ((M + 1) - M * x) / (1 - x) ** 2 def tail_np(C, x, p): # C * sum_{n>M} n^p x^n <= C (M+1)^p x^{M+1} p!/(1-x)^{p+1} fac = {2: 2, 3: 6, 4: 24}[p] # (1+j)^p <= p! C(j+p, p) return C * mp.mpf(M + 1) ** p * x ** (M + 1) * fac / (1 - x) ** (p + 1) T_psi = tail_n1(3 * s3, qa) # psi1 value tail at qa T_eich = tail_np(1080 * s3, qa, 3) # Eichler value tail at qa # absolute-series partial sums at qb (majorants for values + derivatives; # triangle inequality on the ACTUAL computed coefficients — rigorous) PAp = abs(psi1[0]) PAd = EAp = EAd = TA2 = mp.mpf(0) xp = mp.mpf(1) # qb^(n-1) running xpp = mp.mpf(1) # qb^(n-2) running (from n = 2) for n in range(1, M + 1): a1, a3 = abs(psi1[n]), abs(f3[n]) PAd += n * a1 * xp EAd += a3 / n * xp if n >= 2: TA2 += n * (n - 1) * abs(tser[n]) * xpp # |t''| partial at qb xpp *= qb xp *= qb PAp += a1 * xp EAp += a3 / mp.mpf(n) ** 2 * xp PA = PAp + tail_n1(3 * s3, qb) PAd = PAd + tail_np(3 * s3, qb, 2) / qb EA = EAp + tail_np(1080 * s3, qb, 3) EAd = EAd + tail_np(1080 * s3, qb, 4) / qb c0 = (mp.mpf(3) / 2) * s3 * (mp.zeta(2, mp.mpf(1) / 3) - mp.zeta(2, mp.mpf(2) / 3)) / 9 DJ = PAd * (c0 + EA) + PA * EAd # |dJ0/dq| on |z| <= qb T_t, rho, St = _t_scale_pick(qa, M, r0) y = qa / rho T_td = (St / qa) * y ** (M + 1) * ((M + 1) - M * y) / (1 - y) ** 2 x2 = qb / rho # |t''| on |z| <= qb: computed partial sum (TA2) + majorant tail (n > M, # n(n-1) <= n^2, |t_n| <= St/rho^n) D2 = TA2 + (St / qb ** 2) * mp.mpf(M + 1) ** 2 * x2 ** (M + 1) \ * 2 / (1 - x2) ** 3 dser = [k * tser[k] for k in range(1, len(tser))] r_res = abs(horner(tser, q) - tval) # measured Newton residual a_lo = abs(horner(dser, q)) - T_td # lower bound on |t'(q_C)| bad = mp.mpf(10) ** (mp.mp.dps + 100) # sentinel: forces escalation/raise if a_lo > 0: m_tot = T_t + r_res h = m_tot * D2 / a_lo ** 2 dq_trunc = 2 * T_t / a_lo dq_num = 2 * r_res / a_lo if not (h <= mp.mpf("0.25") and dq_trunc + dq_num <= r0): dq_trunc = dq_num = bad # Kantorovich bracket not certified else: h, dq_trunc, dq_num = bad, bad, bad err_trunc = T_psi * (c0 + EA) + PA * T_eich + DJ * dq_trunc err_num = DJ * dq_num return {"T_psi": T_psi, "T_eich": T_eich, "T_t": T_t, "rho": rho, "r_res": r_res, "h": h, "dq_trunc": dq_trunc, "dq_num": dq_num, "err_trunc": err_trunc, "err_num": err_num, "N": M} def f(t, dps=110, nterms=None, cert_out=None): """Public entry point: J^(0)(t) at kinematic point t = p^2/m^2, dps digits. t : number or numeric string; t < 0 (Euclidean -> mpf) or t > 9 (physical +i0 -> mpc). 0 <= t <= 9 raises ValueError (see header). dps : working precision in decimal digits. The q-series length is auto-sized from a measured truncation-tail bound (see loop below), so precision is capped by dps only, never a stored table. nterms : optional q-series length SEED override (skips the measured auto-sizing; the certified refine-until-bound loop below still extends it, fail-closed, if the certified budget misses). cert_out : optional dict, filled with the certified error budget (BOUNDS, see _certify). Certification always runs and RAISES RuntimeError (fail-closed) if the certified truncation budget cannot clear 10^-(dps+TAIL_GUARD), or if the Newton/roundoff nome floor reaches 10^-(dps-NUM_GUARD). Refuses |q_C| >= QMAX=0.92 (t -> 0^-, -inf, or 9^+) explicitly.""" old = mp.mp.dps mp.mp.dps = dps try: tval = mp.mpf(t) N = nterms if nterms else max(NTERMS, 2 * dps) psi1, f3, tser = get_series(N) q = nome_for_t(tser, tval) qa = abs(q) if qa >= QMAX: raise ValueError( f"|q_C(t)| = {float(qa):.3f} >= QMAX = {QMAX}: too close to a cusp " "(t -> 0^-, -inf, or 9^+); outside this script's series domain") # Auto-size the q-series by a MEASURED truncation-tail bound. The # hauptmodul's eta-quotient coefficients grow like 10^(~1.5 sqrt(n)) # (partition-type), so a bare dps/-log10|q| estimate undersizes the # series; bound the tail with the actually-computed top coefficients # and extend until it clears the requested dps. for _ in range(8): M = len(tser) - 1 cmax = max(abs(tser[M]), abs(f3[M]), abs(psi1[M]), mp.mpf(1)) tail = cmax * qa ** M / (1 - qa) if nterms is not None or tail < mp.mpf(10) ** (-(dps + 8)): break slope = -mp.log10(qa) - mp.log10(cmax) / M # net decay, digits/term if slope < mp.mpf("0.02"): raise ValueError( f"|q_C(t)| = {float(qa):.3f}: series would need >{50 * (dps + 8)} " "terms at this dps; outside this script's practical domain") M = int(mp.ceil((dps + 18) / slope)) + 50 psi1, f3, tser = get_series(M) dser = [k * tser[k] for k in range(1, len(tser))] q = _newton(q, tval, tser, dser) # re-polish nome on longer series qa = abs(q) # ---- certified truncation budget: refine until BOUND, fail-closed -- tol = mp.mpf(10) ** (-(dps + TAIL_GUARD)) esc = 0 cert = _certify(psi1, f3, tser, q, tval) while not (cert["err_trunc"] < tol): esc += 1 if esc > CERT_CAP: raise RuntimeError( "FAIL-CLOSED: certified truncation bound " f"{mp.nstr(cert['err_trunc'], 3)} >= tol {mp.nstr(tol, 3)}" f" at t={t}, dps={dps}, N={len(tser) - 1} after " f"{CERT_CAP} x1.5 series extensions") M = int(1.5 * (len(tser) - 1)) + 50 psi1, f3, tser = get_series(M) dser = [k * tser[k] for k in range(1, len(tser))] q = _newton(q, tval, tser, dser) cert = _certify(psi1, f3, tser, q, tval) num_gate = mp.mpf(10) ** (-(dps - NUM_GUARD)) if not (cert["err_num"] < num_gate): raise RuntimeError( "FAIL-CLOSED: nome Newton/roundoff floor " f"{mp.nstr(cert['err_num'], 3)} >= gate " f"{mp.nstr(num_gate, 3)} at t={t}, dps={dps} (residual " f"{mp.nstr(cert['r_res'], 3)})") cert["esc"], cert["dps"] = esc, dps if cert_out is not None: cert_out.update(cert) return J0_from_q(q, psi1, f3) finally: mp.mp.dps = old def agreed(v, oracle): """Recomputed agreement: -log10(|v - oracle| / |oracle|).""" d = abs(v - oracle) return mp.mp.dps if d == 0 else int(mp.floor(-mp.log10(d / abs(oracle)))) ORACLE_DIGITS = sum(c.isdigit() for c in REF[-3]) # 128 stored oracle digits (counted, not assumed) def _oracle_for(tval, dps): """Stored-oracle value for a gate point (parsed at full excerpt precision), or None. Never fabricates: only t = -3, -5, -9, 12 have stored oracle literals.""" with mp.workdps(max(dps, ORACLE_DIGITS) + 10): if tval in (-3, -5, -9): return mp.mpf(REF[int(tval)]) if tval == 12: return mp.mpc(mp.mpf(REF12[0]), mp.mpf(REF12[1])) return None def _print_cert(cert, indent=" "): """[BOUND]/[cert] lines for a completed evaluation (all values BOUNDS).""" d = cert["dps"] print(f"{indent}[BOUND] certified q-series truncation tails: psi1 <= " f"{mp.nstr(cert['T_psi'], 3)}, Eichler <= " f"{mp.nstr(cert['T_eich'], 3)}, t-series <= " f"{mp.nstr(cert['T_t'], 3)} (rho = {mp.nstr(cert['rho'], 3)}); " f"nome dq_trunc <= {mp.nstr(cert['dq_trunc'], 3)} " f"(Newton-Kantorovich h = {mp.nstr(cert['h'], 3)} <= 1/4)") print(f"{indent}[cert] |J0 truncation error| <= " f"{mp.nstr(cert['err_trunc'], 3)} < 1e-{d + TAIL_GUARD} -- BOUND, " f"not estimate (N = {cert['N']}, escalations = {cert['esc']}, " f"fail-closed at cap {CERT_CAP}); Newton/roundoff nome floor <= " f"{mp.nstr(cert['err_num'], 3)} (raising gate 1e-{d - NUM_GUARD}); " f"end-to-end roundoff gated by the D/D+{CRANK_STEP} crank") def _crank_gate(v, t, dps, fails): """In-run two-precision crank gate D -> D+CRANK_STEP (raising).""" v40 = f(t, dps=dps + CRANK_STEP) with mp.workdps(dps + CRANK_STEP + 20): dcr = agreed(v, v40) bar = dps - CRANK_MARGIN ok = dcr >= bar print(f" [gate] crank D/D+{CRANK_STEP} (dps {dps} -> " f"{dps + CRANK_STEP}) at t={t}: agreement = {dcr} digits " f"(raising bar {bar}) {'PASS' if ok else 'FAIL'}") if not ok: fails.append(f"crank D/D+{CRANK_STEP} at t={t}: {dcr} < {bar}") if __name__ == "__main__": import argparse import time ap = argparse.ArgumentParser( description="eps^0 equal-mass two-loop sunrise J^(0)(t), evaluated at " "runtime from the Gamma_1(6) closed form. Domain: t < 0 or " "t > 9 (+i0). No arguments -> gate demo + dps-doubling demo.") ap.add_argument("--point", action="append", metavar="T", help="kinematic point t = p^2/m^2 (repeatable); t<0 or t>9") ap.add_argument("--dps", type=int, default=110, help="decimal digits of working precision (default 110)") args = ap.parse_args() if args.point: fails = [] cranks = [] for ts in args.point: tv = float(ts) t0 = time.perf_counter() cert = {} try: v = f(ts, dps=args.dps, cert_out=cert) except ValueError as e: print(f" t={ts}: DOMAIN ERROR -- {e}") continue dt = time.perf_counter() - t0 print(f" t={ts}, dps={args.dps} [{dt:.1f} s wall, measured]") if isinstance(v, mp.mpc): print(f" Re = {mp.nstr(v.real, args.dps)}") print(f" Im = {mp.nstr(v.imag, args.dps)}") else: print(f" J0 = {mp.nstr(v, args.dps)}") _print_cert(cert) o = _oracle_for(tv, args.dps) if o is not None: d = agreed(v, o) print(f" agreement vs stored {ORACLE_DIGITS}-digit oracle " f"(recomputed) = {d} digits " f"[cap = min(dps, {ORACLE_DIGITS})]") bar = min(args.dps, ORACLE_DIGITS) - GATE_MARGIN if d < bar: fails.append(f"oracle gate at t={ts}: {d} < {bar}") else: print(" (no stored oracle literal at this point -- value " "printed without an agreement claim)") cranks.append((ts, v)) # crank evals (higher dps) run AFTER all default-dps evaluations so # the latter see pristine constant caches (byte-identity of values) for ts, v in cranks: _crank_gate(v, ts, args.dps, fails) if fails: print("[gate] FAIL (fail-closed raising gates):") for m in fails: print(" - " + m) raise SystemExit(1) print(f"[gate] PASS: all raising gates clear (crank D/D+{CRANK_STEP} " f">= dps-{CRANK_MARGIN}; oracle table, where stored, >= " f"min(dps,{ORACLE_DIGITS})-{GATE_MARGIN})") raise SystemExit(0) # ---- default: gate demo (quick), then the dps-doubling honesty demo ---- # All gates below RAISE (nonzero exit) on miss — fail-closed. A 1e-30 # mutation of any REF/REF12 oracle string trips the table gate (rc=1). fails = [] gate_bar = min(args.dps, ORACLE_DIGITS) - GATE_MARGIN t_start = time.perf_counter() print(f"eps^0 closed form evaluated at runtime (dps={args.dps}, q-series auto-sized)") print(f"vs held-out certified oracle ({ORACLE_DIGITS}-digit excerpts, provenance in header):\n") cranks = [] # (t, v): crank evals run AFTER the demos so the for t in (-3, -5): # default-dps values see pristine constant caches cert = {} v = f(t, dps=args.dps, cert_out=cert) o = _oracle_for(t, args.dps) d = agreed(v, o) print(f" t={t} (Euclidean):") print(f" computed = {mp.nstr(v, 50)}") print(f" oracle = {REF[t][:52]}...") print(f" agreement (recomputed) = {d} digits\n") _print_cert(cert) if d < gate_bar: fails.append(f"oracle gate at t={t}: {d} < {gate_bar}") cranks.append((t, v)) cert12 = {} v12 = f(12, dps=args.dps, cert_out=cert12) o12 = _oracle_for(12, args.dps) d12 = agreed(v12, o12) print(" t=12 (physical, +i0; nome analytically continued at runtime):") print(f" computed Re = {mp.nstr(v12.real, 50)}") print(f" computed Im = {mp.nstr(v12.imag, 50)}") print(f" oracle Re = {REF12[0][:52]}...") print(f" oracle Im = {REF12[1][:52]}...") print(f" agreement (recomputed) = {d12} digits\n") _print_cert(cert12) if d12 < gate_bar: fails.append(f"oracle gate at t=12: {d12} < {gate_bar}") cranks.append((12, v12)) gate_s = time.perf_counter() - t_start print(f"gate-demo wall time: {gate_s:.1f} s (measured)\n") # dps-doubling demo: same point, twice the precision -> agreement digits # must GROW until they hit the stored-oracle-excerpt cap. t0 = time.perf_counter() d1 = agreed(f(-3, dps=args.dps), _oracle_for(-3, args.dps)) v2 = f(-3, dps=2 * args.dps) d2 = agreed(v2, _oracle_for(-3, 2 * args.dps)) dbl_s = time.perf_counter() - t0 print(f"dps-doubling check at t=-3 (dps {args.dps} -> {2 * args.dps}):") print(f" agreement at dps={args.dps:<4d}= {d1} digits") print(f" agreement at dps={2 * args.dps:<4d}= {d2} digits") print(f" stored oracle excerpt = {ORACLE_DIGITS} digits, so the doubled-dps") print(f" agreement is CAPPED at ~{ORACLE_DIGITS} by the oracle string, not the") print(" method (more digits = more q-series terms, auto-sized).") print(f"doubling-demo wall time: {dbl_s:.1f} s (measured)\n") dbl_bar = min(2 * args.dps, ORACLE_DIGITS) - 10 if d2 < dbl_bar: fails.append(f"doubling gate: {d2} < {dbl_bar}") if d2 < d1 - 3: # -3: both capped by the oracle string fails.append(f"doubling gate: agreement shrank ({d1} -> {d2})") # in-run D/D+40 two-precision crank gates (raising; run last, see above) for t, v in cranks: _crank_gate(v, t, args.dps, fails) print() print(f"Note: agreement above is limited by mp.dps and the {ORACLE_DIGITS}-digit oracle") print("excerpts, not the method. The full offline certified reveal gate") print("(1024-bit ball arithmetic, 6000 q-terms) agreed at 304-305 digits on") print("all 50 ball-overlap tests -- an ARCHIVED statistic, not computed here.") if fails: print("\n[gate] FAIL (fail-closed raising gates):") for m in fails: print(" - " + m) raise SystemExit(1) print(f"\n[gate] PASS: oracle table >= min(dps,{ORACLE_DIGITS})-" f"{GATE_MARGIN} at 3/3 points; crank D/D+{CRANK_STEP} >= " f"dps-{CRANK_MARGIN} at 3/3; doubling >= min(2*dps," f"{ORACLE_DIGITS})-10 and non-shrinking. A 1e-30 mutation of any " "stored oracle string exits nonzero (mutation rc-gate).")