#!/usr/bin/env python3 """FRW-triangle (one-loop three-site cosmological correlator, a=2) -- live evaluator for the eps-form campaign's boundary-free deliverables, with live independent oracles. Everything gated here is computed at runtime from EXACT inputs only (exact bivariate-rational connection, exact symbolic gauge chain, exact word lists, closed-form periods); no finite-precision constants enter any computation. What runs at runtime (data: cosmo-data.json.gz, built from this work's epsform_2026-07-04 artifacts; A16 artifact sha256 de035018b116c79a5843aff01275085c840e75bb39ce3050da6c8cffdb1836b7): Point 1 Closed L2 periods. varpi0 = K(m)/sqrt(1-9 lam^2), varpi1 = K(1-m)/sqrt(1-9 lam^2), m = (lam^2-1)/(9 lam^2-1), Wronskian W = pi/(2 lam D4), D4 = (lam^2-1)(9 lam^2-1). Live oracles: (i) direct quadrature of the defining integral K(m) = int_0^{pi/2} dtheta/sqrt(1-m sin^2 theta) (tanh-sinh; never calls ellipk), (ii) L2-annihilation residual by high-order finite differences, (iii) the Wronskian identity (branch-sensitive). Point 2 The symbolic-eps connection A16(lam,eps) (16x16, exact bivariate rationals). (i) EXACT check (Fraction arithmetic on a full bidegree grid -- a rigorous polynomial-identity proof, not numerics): A16[z5,z5] == eps * dlog D4. (ii) The proven sqrt(Q4) twist: at each root of Q4 = 20L^4+39L^3+29L^2+9L+1 (roots found live at working precision) the residue matrix of A16 has rank 1 with sole nonzero eigenvalue eps - 1/2; the residue is computed live from the rational data and gated against eps - 1/2. Point 3 Boundary-free per-eps-order holonomy of the eps-factorized system (the campaign's headline gate). Psi(x;eps) solves the rotated (eps-form) system with Psi(1/2) = I -- the boundary is the IDENTITY, exactly; no boundary constants enter. Prediction route ("this work"): the stored explicit word lists (10/254/3437/47077 words at eps^0..3 over 272 dressed letters + Eichler n_k composites) evaluated as iterated integrals over the closed periods by composite Chebyshev-Lobatto spectral quadrature (acb ball arithmetic). Reference route (live independent oracle): arb Taylor-ODE transport of the ORIGINAL SUB9 system -- the 9x9 block of the stored A16 -- at 2M scalar eps nodes, rotated by the exact gauge chain T9 = blockdiag(U1*Tg, Tzc, 1) (T9 is Laurent in eps, so per-order references are extracted by an eps-Vandermonde solve, exactly as in the campaign gate), with Eichler n_k(7/10) from an independent acb-series Taylor stepper (L2-recurrence route, not the spectral route). Gate: per-order min agreed digits over all nonzero entries, orders eps^0..eps^ORD. Point 4 Correlator boundary values (the lambda=1/2 boundary 9-vector at the campaign point a=2, lam=1/2, eps=-1/20) -- CACHE ONLY, wired 2026-07-06. The values are computable to arbitrary precision from the shipped regeneration chain (oracle_eps/symck/boundary_assemble -- written integrals, standalone python); the deep values printed here are cached from a long high-precision two-leg certified run of the shipped regeneration chain and are clamped to their PER-COMPONENT certified digits (42.92..46.78 d; min-component floor 42.92 d). The cache is sha256-pinned and value-level checked at load: any mutation raises (rc != 0). The cache is NEVER the definition -- see --boundary-recompute, which regenerates any component from the chain at a requested dps (measured costs quoted in --help). Retained literals (gate-only provenance, never used in any computation): campaign_gate table inside cosmo-data.json.gz -- the campaign's measured per-order digits (89.3/87.5/78.9/78.8 at its dps80 leg) for comparison printing only; cosmo-boundary-cache.json -- the sha-pinned deep boundary CACHE (Point 4 display + recompute gates only; never enters Points 1-3). Interface: default run Points 1-4 at --dps (measured wall time printed) --dps D working target digits (default 40); Point-3 routes run at internally scaled precisions (grid D+15, ODE/fit D+90) --point 'lam=3/5,eps=1/7' Point-1 lambda (certified domains (0,1/3) and [1/2,7/10]) and Point-2 eps (any rational) --orders K gate holonomy orders eps^0..eps^K (default 2; K=3 is allowed but slow -- a measured projection is printed before it runs) --check integrity extras (all fail-closed raises): (i) Point-3 rerun at D+40 -- agreement digits must GROW (calibrated +25 floor) and every gated entry must agree with the D run to > D digits (value-level rerun at two genuinely different internal depths); (ii) mutation control (one word scaled by 1+1e-30 -- the gate must collapse to ~30 d); (iii) n_1..n_4 cross-checked LIVE against direct certified mp.quad of the Eichler kernels. --skip-holonomy Points 1-2 + 4 only (no python-flint needed) --boundary-recompute e6,e7,... regenerate boundary components from the chain (the DEFINITION path) at --recompute-dps and gate them against the cache; see the measured-cost table in --help epilog (--chain-dir / COSMO_BOUNDARY_CHAIN locates the chain) Requirements: mpmath, sympy, python-flint (>=0.8, for Point 3). Charter-v2 certification (2026-07-05 axis3 wave2): Every truncation on the value path sits inside a refine-until-bound loop: the legacy depth formulas (Chebyshev-Lobatto degree N, Taylor-ODE depths NT/NT2, quadrature maxdegree, eps-Vandermonde node count M) are STARTING SEEDS only. Each loop accepts only when its certified/measured error bound beats 10^-(dps+guard), escalates by EXACT continuation of the same rule (x1.5 depth, +2 eps-nodes), and raises CertError at its cap naming step/x0/h/bound/tol/N/cap -- a value the loop did not certify is never printed. Tail bounds use the axis-3 pilot construction (WIRING_LOG items 12-13): max |c_n h^n| over the trailing 8 terms x r/(1-r), with r certified at runtime from the live singularity set (all letter-denominator roots, D4 roots, Q4 roots -- found live; Bernstein-ellipse ratio for the Chebyshev panels). Per-step bounds sum to err_total and reach value level through inf-norms of the exact prefactors (T9 rotation norms, |Vm^-1| row sums of the eps-Vandermonde). Agreement gates (K-integral oracle, L2 annihilation, Wronskian, twist residue, holonomy bar, --check legs) are RAISES with measured >=1e6 headroom. Ball radii of every flint->mpmath conversion are part of the certified budgets. Guards are calibrated on measured healthy runs so a healthy run never escalates: values at default dps are bit-identical to the 2026-07-05 build. Honest limitations (see also the page's honest-scope block): * The lambda=1/2 boundary 9-vector (Point 4) is a CACHE, not a runtime derivation: computable to arbitrary precision from the shipped chain; deep values cached from a long high-precision certified run of that chain. Its printed digits are clamped to the PER-COMPONENT two-leg certificates (NOT a blanket floor). The cache currently sits at the 42.92 d certified min-component floor; a deeper (mo7) rerun would refresh it WITHOUT rewiring (same file + sha-pin update). Agreement of any recompute with the cache is capped near ~43-47 d by the measured mo5 inner-fallback systematic (~1e-43, recorded with the cache) until that refresh. * Cached boundary values never enter Points 1-3 (which remain zero-finite-precision-input); Point 4 is display + recompute gates. * The holonomy path is the campaign's frozen gate path [1/2, 7/10] (the word lists are path-specific by convention). * eps^3 (47077 words) is gated by the same machinery but is opt-in (--orders 3) for runtime reasons. Changelog: 2026-07-05 created (from the epsform_2026-07-04 result records; assessment + build for the BootLoops cosmo page). 2026-07-05b (axis3 wave2) charter-v2 hardening: refine-until-bound fail-closed loops on every value-path truncation (legacy formulas demoted to starting seeds => default outputs bit-identical); agreement/Wronskian/L2/residue gates promoted to raises with measured >=1e6 headroom; --check rerun now D/D+40 with value-level rerun agreement + live n_k mp.quad cross-check; ball radii of all flint->mp conversions budgeted; the two flint footgun fixes RETAINED untouched (endpoint- stable M2 without acos; nk_taylor local-precision PI2); TODO consume hook added for the deep boundary values -- lambda=1/2 boundary 9-vector exclusion unchanged, as disclosed. 2026-07-06 deep-boundary fold: correlator-values exclusion LIFTED. Point 4 = cached deep boundary 9-vector (sha-pinned cosmo-boundary-cache.json = byte-copy of this work's step2_boundary_deep.json), printed at per-component certified digits; --boundary-recompute added (chain = definition; measured tier costs in --help). Points 1-3 value paths untouched (default non-correlator output byte-identical). """ import argparse, cmath, gzip, json, os, sys, time from fractions import Fraction as Fr from math import lcm import mpmath as mp import sympy as sp HERE = os.path.dirname(os.path.abspath(__file__)) T0 = time.time() def digits(a, b): """-log10 |a-b|/max(|a|,|b|): measured agreement in decimal digits.""" a, b = mp.mpmathify(a), mp.mpmathify(b) if a == b: return float('inf') den = max(abs(a), abs(b)) return float(-mp.log10(abs(a - b) / den)) def load_bundle(): with gzip.open(os.path.join(HERE, 'cosmo-data.json.gz'), 'rb') as f: return json.loads(f.read()) # =========================================================================== # charter-v2 certification infrastructure (2026-07-05 axis3 wave2) # =========================================================================== class CertError(RuntimeError): """Fail-closed certification failure: a refine-until-bound loop hit its cap, a promoted agreement gate fired, or a ball radius blew its budget.""" class _Escalate(Exception): """internal: prediction-route bound miss -> escalate N and retry.""" # ---- certification guards (calibration measured 2026-07-05; every # guard leaves >=1e6 measured headroom over the healthy bound so the fast # path NEVER escalates on a healthy run) ---- PRED_GUARD = 10 # prediction route: worst-entry certified bound < 10^-(dps+10) # (measured healthy at dps 40: ~1e-59 vs tol 1e-50) REF_GUARD = 12 # reference route: per-eps-order value-level bound < 10^-(dps+12) REF_STEP_SLACK = 5 # ODE/nk per-step tail must beat 10^-(workdps-5) QUAD_GUARD = 6 # point-1 K oracle + --check n_k quadratures: est < 10^-(D+6) _ROOTS_CACHE = {} def _mpf2(me): """(man, exp) -> EXACT dyadic mpf (never float-underflows at any dps).""" m, e = me return mp.mpf((int(m), int(e))) def _absu(z): """certified mpf upper bound on |acb z| (exact dyadic; error terms).""" u = z.abs_upper() return abs(_mpf2(u.mid().man_exp())) + _mpf2(u.rad().man_exp()) def _absuf(z): """fast float upper bound on |acb z| (O(1)-scale propagation maxima).""" try: v = abs(complex(z)) + float(z.real.rad()) + float(z.imag.rad()) except (ValueError, OverflowError, TypeError): return float('inf') if v == 0.0: vv = _absu(z) if vv == 0: return 0.0 fv = float(vv) return fv if fv > 0 else 1e-300 # underflow-safe upper bound return v * (1 + 1e-12) def _mat_maxu(Mt): """certified mpf upper bound on the max |entry| of a 9x9 arb_mat.""" w = mp.mpf(0) for i in range(9): for j in range(9): a = Mt[i, j] v = abs(_mpf2(a.mid().man_exp())) + _mpf2(a.rad().man_exp()) if v > w: w = v return w def _mnorm(Mt): """inf-norm (max row sum) of an mpmath matrix.""" return max(sum(abs(Mt[i, j]) for j in range(Mt.cols)) for i in range(Mt.rows)) def _bernstein_r(lo, hi, sings): """certified Chebyshev-coefficient decay ratio r = 1/rho for the panel [lo,hi]: rho = min over the live singularity set of |w + sqrt(w^2-1)|, w = (2z-lo-hi)/(hi-lo) (Bernstein-ellipse parameter).""" rho = None for z in sings: w = (2*z - lo - hi) / (hi - lo) s = cmath.sqrt(w*w - 1) cand = max(abs(w + s), abs(w - s)) rho = cand if rho is None else min(rho, cand) r = (1.0 / rho) * (1 + 1e-9) if not 0 < r < 0.9: raise CertError(f'panel [{lo:.4f},{hi:.4f}]: certified decay ratio r={r:.3f} ' 'not in (0,0.9) -- a singularity sits on/near the path') return r _DTQUAD = {'probed': False, 'fn': None} def _detransport_quad(): """tools/detransport.quad.quad_refine when the repo tree is present (the standardized E3 helper); on a standalone copy this returns None and the inline maxdegree-ladder below provides the same double-refinement fail-closed semantics (escalation by exact continuation of the nested tanh-sinh rule, mpmath's own successive-level estimate as the gate).""" if not _DTQUAD['probed']: _DTQUAD['probed'] = True try: import importlib.util for up in ('..', '../..', '../../..', '../../../..'): cand = os.path.abspath(os.path.join(HERE, up, 'tools', 'detransport', 'quad.py')) if os.path.exists(cand): spec = importlib.util.spec_from_file_location('_cosmo_dt_quad', cand) mod = importlib.util.module_from_spec(spec) spec.loader.exec_module(mod) _DTQUAD['fn'] = mod.quad_refine break except Exception: _DTQUAD['fn'] = None return _DTQUAD['fn'] def _certified_quad(f, pts, tol, what): """fail-closed mp.quad: the value is accepted ONLY when mpmath's own successive-level error estimate beats tol; refinement = maxdegree ladder (exact continuation of the nested tanh-sinh rule; the default maxdegree is the STARTING seed), then the detransport quad_refine helper when present; CertError at the cap. Returns (value, estimate).""" est_hist = [] for mdeg in (6, 8, 10): # 6 = mpmath default seed v, est = mp.quad(f, pts, error=True, maxdegree=mdeg) est_hist.append(est) if mp.isfinite(est) and est < tol and mp.isfinite(abs(v)): return v, est qr = _detransport_quad() if qr is not None: try: dps_t = max(15, int(mp.floor(-mp.log10(tol))) - 10) cert = qr(f, list(pts), dps_t, guard=10, full_output=True) return cert.value, cert.err_bound except Exception as e: raise CertError(f'{what}: quadrature refine-until-bound FAILED at cap ' f'(maxdegree ladder 6/8/10 estimates ' f'{[mp.nstr(x, 3) for x in est_hist]} >= tol {mp.nstr(tol, 3)}; ' f'detransport.quad_refine escalation also failed: {e})') raise CertError(f'{what}: quadrature refine-until-bound FAILED at cap maxdegree=10 ' f'(estimates {[mp.nstr(x, 3) for x in est_hist]} >= tol ' f'{mp.nstr(tol, 3)}; tools/detransport not present for deeper escalation)') # --------------------------------------------------------------------------- # Deep boundary values: lambda=1/2 boundary 9-vector -- WIRED 2026-07-06. # The former exclusion hook is # lifted: the deep boundary values are consumed as CACHE ONLY. # cache = cosmo-boundary-cache.json, a byte-copy of this work's # step2_boundary_deep.json (sha256-pinned below); per-component # two-leg certificates ride with the values and CLAMP all printing. # defn = the regeneration chain (oracle_eps / symck / boundary_assemble # + the cosmoshard exact shard layer, shipped in chain/) -- exposed # here as --boundary-recompute (measured costs in --help). # Fail-closed: missing cache, sha mismatch (ANY byte mutation, e.g. 1e-30 on # one component), or a value-level integrity miss all raise CertError. # --------------------------------------------------------------------------- BOUNDARY_CACHE_PATH = os.path.join(HERE, 'cosmo-boundary-cache.json') BOUNDARY_CACHE_SHA256 = \ 'fffefe6760d9fc4115bf9d0b3b722aec4f763cef926da55db44ab8780994650a' BOUNDARY_ORDER = ['e6', 'e7', 'e8', 'e9', 'z2', 'z3', 'z4', 'z5', 'z1'] def load_boundary_cache(): """sha-pinned + value-level-checked load of the cached deep boundary. The pin makes the mutation gate structural: a 1e-30 edit of any cached component changes the file bytes -> CertError (rc != 0). On top of the pin, the convergent components' stored two-leg agreement d_2prec is re-MEASURED from the shipped value/value_lo pair, and the divergent certificates must satisfy d_cert == min(d_L, d_2prec).""" if not os.path.exists(BOUNDARY_CACHE_PATH): raise CertError('boundary cache cosmo-boundary-cache.json MISSING next to the ' 'evaluator -- it ships with the page (byte-copy of this work\'s ' 'step2_boundary_deep.json). Restore it or regenerate from the ' 'chain (--boundary-recompute; definition path).') import hashlib raw = open(BOUNDARY_CACHE_PATH, 'rb').read() sha = hashlib.sha256(raw).hexdigest() if sha != BOUNDARY_CACHE_SHA256: raise CertError(f'boundary cache sha256 MISMATCH: {sha[:20]}... != pinned ' f'{BOUNDARY_CACHE_SHA256[:20]}... -- the cache is CACHE ONLY ' '(the definition is the regeneration chain); a mutated/edited ' 'cache is refused, never silently consumed.') cb = json.loads(raw) with mp.workdps(140): for k, ent in cb['convergent'].items(): d = digits(mp.mpf(ent['value']), mp.mpf(ent['value_lo'])) if not abs(d - ent['d_2prec']) < 0.5: raise CertError(f'boundary cache value-level check FAILED: {k} ' f'value/value_lo agreement {d:.3f} d does not reproduce ' f'the stored two-leg cert d_2prec={ent["d_2prec"]}') certs = [ent['d_2prec'] for ent in cb['convergent'].values()] for k, ent in cb['divergent'].items(): if abs(ent['d_cert'] - min(ent['d_L'], ent['d_2prec'])) > 1e-9: raise CertError(f'boundary cache: {k} d_cert != min(d_L, d_2prec)') certs.append(ent['d_cert']) if abs(min(certs) - cb['min_component_cert']) > 1e-9: raise CertError('boundary cache: min_component_cert inconsistent with ' 'per-component certificates') if cb.get('positive_control') != 'PASS': raise CertError('boundary cache: positive control vs the stored 31d vector is ' f'not PASS ({cb.get("positive_control")!r})') return cb def point4(cb): print('== Point 4: correlator boundary values at lambda=1/2 (CACHE, clamped ' 'to per-component certs) ==') print(' computable to arbitrary precision from the shipped chain ' '(oracle_eps/symck/boundary_assemble); deep values cached from a ' 'long high-precision certified run of that chain (two independent ' 'legs, 2026-07-06). The cache is NOT the definition: ' '--boundary-recompute regenerates any component from the ' 'written integrals at requested dps (measured costs: --help).') pt = cb['meta']['point'] print(f" boundary point: a={pt['a']}, lam={pt['lam']}, eps={pt['eps']}, " f"c={pt['c']} (campaign boundary; per-component clamp = its own " f"two-leg certificate, NOT the blanket floor)") for k in BOUNDARY_ORDER: fam = 'convergent' if k in cb['convergent'] else 'divergent' ent = cb[fam][k] cert = ent['d_2prec'] if fam == 'convergent' else ent['d_cert'] nd = int(cert) # printed-digit clamp, per component with mp.workdps(nd + 10): vs = mp.nstr(mp.mpf(ent['value']), nd) if fam == 'convergent': certline = f'two-leg joint cross d_2prec={ent["d_2prec"]}' else: certline = (f'd_cert={ent["d_cert"]} = min(d_L={ent["d_L"]}, ' f'd_2prec={ent["d_2prec"]}); d_tailK={ent["d_tailK"]}') print(f' {k:2s} = {vs}') print(f' [certified {cert:.2f} d ({nd} printed): {certline}; vs ' f'stored-31d control {cb["positive_control_vs_banked"][k]} d]') print(f' min-component certified floor: {cb["min_component_cert"]} d. ' f'The cache currently sits at the 42.92 d certified floor; a deeper ' f'(mo7) rerun refreshes this cache ' f'WITHOUT rewiring (same file + sha-pin update).') print(' Point 4: PASS (sha-pin + value-level integrity checks)') return True # --------------------------------------------------------------------------- # --boundary-recompute: the DEFINITION path. Regenerates boundary components # from the chain's written integrals via the exact shard layer (cosmoshard), # then gates the fresh value against the cache. Sharding is exact # decomposition (independent quadrature segments / disk wedges / tail u-nodes # summed losslessly), so any --recompute-jobs parallelism is bit-equivalent # to a serial run. # --------------------------------------------------------------------------- RECOMPUTE_HELP = """ --boundary-recompute costs (MEASURED, quoted RELATIVE to the probe tier -- one probe-tier conv-pair run = 1 unit; the two lines labeled PROJECTED are rate-fits from measured probe shards, never guesses. Rate-fit ONE shard on your own machine first to price a unit locally; sharding is exact, so any --recompute-jobs parallelism divides the wall without changing the values): conv pair {e6,e7} (one run computes both): probe tier (--recompute-dps <=41: nj48 mm4 mo4, corner/wing pieces mo3, nth128 rdeg5): 1 unit (the baseline; dps36 measured ~0.9 unit) [MEASURED, 2026-07-06 gate; agreement vs cache sits at the measured mm4/mm5 floor ~22.2 d, two-precision stable; fresh-pair internal agreement ~30.5 d] control tier (42-47: nj48 mm5 mo5 nth128 rdeg5): ~1.8 units [agreement vs cache capped ~31.3 d by the measured mm5/mm6 floor (31.24 d)] mid tier (48..65: nj64 mm6 mo5 nth192 rdeg6): PROJECTED ~7-18 units (rate-fit from measured shards; NOT yet run end-to-end) deep tier (>65: legA knobs nj96 mm6 mo5 nth256 rdeg7): ~49 units at dps80; the independent legB config (dps85 nj112 mm6 mo6 nth320 rdeg8) measured ~146 units divergent 7-vector {e8,e9,z2,z3,z4,z5,z1} (one near+tail run computes all 7): control tier: ~1.9 units (near + tail) mid tier: PROJECTED ~x4 control (L15/K90; same rate-fit basis as conv) deep tier: ~54 units (legA); legB measured ~146 units full two-leg certified computation (source of the shipped cache): ~400 units -- a long many-core computation Honesty notes: (i) agreement between ANY recompute and the cache is capped near ~43-47 d by the measured mo5 inner-fallback systematic (~1e-43, recorded with the cache) until a deeper (mo7) refresh; (ii) wall time ~ total work / effective parallel workers -- rate-fit a single shard first on a loaded machine. """ RC_TIERS = { # split='fine' shards the outer segments so every shard fits a short # foreground window (each seg piece pays a FULL outer tanh-sinh ladder -- # measured 295 evals at mo5 / ~150 at mo4 / ~75 at mo3, width-independent # -- so fine splitting multiplies total cpu; the probe tier accepts that # and drops the corner/wing outer ladder one level, mo_cw, with the # per-piece quad estimates recorded and gated in the combine). # split='natural' = the shipped chain's own minimal 5-segment split # (cheapest total cpu; individual shards run tens of minutes -- resumable). 'conv': {'probe': dict(njac=48, mm=4, mo=4, mo_cw=3, nth=128, rdeg=5, bar=20.0, split='fine'), 'control': dict(njac=48, mm=5, mo=5, mo_cw=5, nth=128, rdeg=5, bar=30.0, split='natural'), 'mid': dict(njac=64, mm=6, mo=5, mo_cw=5, nth=192, rdeg=6, bar=38.0, split='natural'), 'deep': dict(njac=96, mm=6, mo=5, mo_cw=5, nth=256, rdeg=7, bar=40.0, split='natural')}, 'near': {'probe': dict(njac=48, mm=4, mo=4, mo_cw=3, nth=96, rdeg=5, L=15, K=75, tdps=100, bar=20.0, split='fine'), 'control': dict(njac=48, mm=5, mo=5, mo_cw=5, nth=96, rdeg=5, L=15, K=75, tdps=100, bar=30.0, split='natural'), 'mid': dict(njac=64, mm=6, mo=5, mo_cw=5, nth=192, rdeg=6, L=15, K=90, tdps=110, bar=38.0, split='natural'), 'deep': dict(njac=96, mm=6, mo=5, mo_cw=5, nth=256, rdeg=7, L=25, K=90, tdps=130, bar=40.0, split='natural')}, } CONV_KEYS = ('e6', 'e7') DIV_KEYS = ('e8', 'e9', 'z2', 'z3', 'z4', 'z5', 'z1') def _rc_tier(dps): if dps <= 41: return 'probe' return 'control' if dps <= 47 else ('mid' if dps <= 65 else 'deep') def _dec_exact(o): """decode a cosmoshard exact-mpf record [sign, man, exp, bc] (or the {'re':..,'im':..} complex form -> real part; imaginary parts of the conv values are quadrature noise at ~1e-100) LOSSLESSLY, independent of the mpmath tuple-constructor signature.""" if isinstance(o, dict): return _dec_exact(o['re']) s, m, e, bc = int(o[0]), int(o[1]), int(o[2]), int(o[3]) with mp.workprec(bc + 16): v = mp.ldexp(mp.mpf(m), e) return -v if s else v def _rc_chain_dir(arg): cands = [arg, os.environ.get('COSMO_BOUNDARY_CHAIN'), os.path.join(HERE, 'chain')] for c in cands: if c and os.path.exists(os.path.join(c, 'cosmoshard.py')): return os.path.abspath(c) raise CertError('regeneration chain not found (need cosmoshard.py + the 8 chain ' 'modules). The chain ships with this page as files/cosmo/chain/; ' 'point --chain-dir or COSMO_BOUNDARY_CHAIN at a copy of that ' 'directory.') def _rc_pieces(lo, hi, cuts): pts = [lo] + list(cuts) + [hi] return [f'{pts[i]}:{pts[i+1]}' for i in range(len(pts) - 1)] def _rc_seg_parts(split, final_hi, ltail_cuts=()): """(part, is_corner_or_wing) pairs. 'natural' = the shipped chain's own 5-segment split (minimal total cpu); 'fine' = window-sized pieces (lam=1/2 geometry: X2mR=0.45, X2=0.5, X2pR=0.55, X3=1).""" parts = [] if split == 'natural': parts += [(p, False) for p in _rc_pieces('0', 'X2mR', [])] parts += [(p, True) for p in _rc_pieces('X2mR', 'X2', [])] parts += [(p, True) for p in _rc_pieces('X2', 'X2pR', [])] parts += [(p, True) for p in _rc_pieces('X2pR', 'X3', [])] parts += [(p, False) for p in _rc_pieces('X3', final_hi, list(ltail_cuts))] else: parts += [(p, False) for p in _rc_pieces('0', 'X2mR', ['1/8', '1/4', '3/8'])] parts += [(p, True) for p in _rc_pieces('X2mR', 'X2', ['91/200', '23/50', '93/200'])] parts += [(p, True) for p in _rc_pieces('X2', 'X2pR', ['101/200', '51/100', '103/200'])] parts += [(p, True) for p in _rc_pieces('X2pR', 'X3', ['5/8', '7/10', '3/4', '4/5', '9/10'])] parts += [(p, False) for p in _rc_pieces('X3', final_hi, list(ltail_cuts) or ['3/2', '3'])] return [('seg:' + p, cw) for p, cw in parts] def _rc_run_jobs(cmds, jobs, quiet): import shutil, subprocess from concurrent.futures import ThreadPoolExecutor pl = shutil.which('prlimit') prefix = [pl, f'--as={20 * 1024**3}'] if pl else [] def run1(ic): i, c = ic out = c[c.index('--out') + 1] if os.path.exists(out): # lossless resume try: json.load(open(out)) return f'[{i+1}/{len(cmds)}] resume-skip {os.path.basename(out)}' except Exception: os.unlink(out) r = subprocess.run(prefix + c, capture_output=True, text=True) if r.returncode != 0: raise CertError(f'recompute shard FAILED rc={r.returncode}: ' f'{" ".join(c[-4:])}\n{r.stderr[-1200:]}') tail = r.stdout.strip().splitlines()[-1] if r.stdout.strip() else '' return f'[{i+1}/{len(cmds)}] {tail}' with ThreadPoolExecutor(max_workers=jobs) as ex: for line in ex.map(run1, enumerate(cmds)): if not quiet: print(' ', line, flush=True) def boundary_recompute(comps, dps, chain_arg, jobs, bar_arg, wdir, cb): import subprocess chain = _rc_chain_dir(chain_arg) cs = os.path.join(chain, 'cosmoshard.py') tier = _rc_tier(dps) comps = [c.strip() for c in comps.split(',') if c.strip()] bad = [c for c in comps if c not in CONV_KEYS + DIV_KEYS] if bad: raise CertError(f'unknown boundary component(s) {bad}; valid: ' f'{",".join(CONV_KEYS + DIV_KEYS)}') want_conv = [c for c in comps if c in CONV_KEYS] want_div = [c for c in comps if c in DIV_KEYS] wdir = os.path.abspath(wdir or f'./cosmo-boundary-recompute-{tier}{dps}') os.makedirs(wdir, exist_ok=True) py = sys.executable t0 = time.time() print(f'== boundary recompute (DEFINITION path): {",".join(comps)} at dps {dps} ' f'(tier {tier}) ==') print(f' chain: {chain}\n workdir: {wdir} jobs: {jobs}') fresh = {} cpu_tot = 0.0 def gjprep(alpha, N, wdps, gj01): tag = f'{alpha.replace("/", "o").replace("-", "m")}_{N}_{wdps}' p = os.path.join(chain, f'gj_cache_{tag}.json') need = True if os.path.exists(p): try: d = json.load(open(p)) need = gj01 and 'gj01' not in d except Exception: need = True if need: cmd = [py, cs, 'gjprep', f'--alpha={alpha}', '--N', str(N), '--wdps', str(wdps)] + (['--gj01'] if gj01 else []) r = subprocess.run(cmd, capture_output=True, text=True) if r.returncode != 0: raise CertError(f'gjprep failed: {r.stderr[-800:]}') print(f' gjprep alpha={alpha} N={N} wdps={wdps} done') if want_conv: kn = RC_TIERS['conv'][tier] gjprep('-11/20', kn['njac'], dps + 20, True) parts = _rc_seg_parts(kn['split'], 'inf', () if kn['split'] == 'natural' else ('3/2', '3')) wed = max(4, kn['nth'] // 16) parts += [(f'disk:{j}:{min(j + wed, kn["nth"])}', False) for j in range(0, kn['nth'], wed)] cmds, outs = [], [] for i, (part, cw) in enumerate(parts): out = os.path.join(wdir, f'conv_{i:03d}.json') outs.append(out) cmds.append([py, cs, 'convshard', '--dps', str(dps), '--njac', str(kn['njac']), '--mm', str(kn['mm']), '--mo', str(kn['mo_cw'] if cw else kn['mo']), '--nth', str(kn['nth']), '--rdeg', str(kn['rdeg']), '--part', part, '--out', out]) print(f' conv plan: {len(cmds)} exact shards (keys e6,e7)') _rc_run_jobs(cmds, jobs, quiet=True) comb = os.path.join(wdir, 'conv_combine.json') r = subprocess.run([py, cs, 'combine', '--out', comb] + outs, capture_output=True, text=True) if r.returncode != 0: raise CertError(f'conv combine failed: {r.stderr[-800:]}') cj = json.load(open(comb)) cpu_tot += cj['cpu_total_s'] with mp.workdps(140): for k in want_conv: fresh[k] = _dec_exact(cj['vals_exact'][k]) print(f' conv combine: cpu_total {cj["cpu_total_s"]} s, quad est ' f'{cj.get("d_quad_est")}') if want_div: kn = RC_TIERS['near'][tier] gjprep('-11/20', kn['njac'], dps + 20, True) L = kn['L'] lt = [c for c in ('3/2', '2', '5/2', '3', '4', '5', '7', '9', '12', '15', '20') if Fr(c) < L] parts = _rc_seg_parts(kn['split'], str(L), lt) wed = max(4, kn['nth'] // 16) parts += [(f'disk:{j}:{min(j + wed, kn["nth"])}', False) for j in range(0, kn['nth'], wed)] cmds, outs = [], [] for i, (part, cw) in enumerate(parts): out = os.path.join(wdir, f'near_{i:03d}.json') outs.append(out) cmds.append([py, cs, 'nearshard', '--dps', str(dps), '--njac', str(kn['njac']), '--mm', str(kn['mm']), '--mo', str(kn['mo_cw'] if cw else kn['mo']), '--L', str(L), '--nth', str(kn['nth']), '--rdeg', str(kn['rdeg']), '--keys', ','.join(DIV_KEYS), '--part', part, '--out', out]) # tail: symck K-sum over u nodes (Nu = K+12), exact u-sharding K, tdps = kn['K'], kn['tdps'] Nu = K + 12 gjprep('-11/20', Nu, tdps + 20, False) gjprep('-1/20', Nu, tdps + 20, False) tout = [] for i, u0 in enumerate(range(0, Nu, 22)): out = os.path.join(wdir, f'tail_{i:02d}.json') tout.append(out) cmds.append([py, cs, 'tailshard', '--K', str(K), '--dps', str(tdps), '--u0', str(u0), '--u1', str(min(u0 + 22, Nu)), '--keys', ','.join(DIV_KEYS), '--out', out]) print(f' near+tail plan: {len(cmds)} exact shards (keys {",".join(DIV_KEYS)}; ' f'L={L}, K={K})') _rc_run_jobs(cmds, jobs, quiet=True) combn = os.path.join(wdir, 'near_combine.json') r = subprocess.run([py, cs, 'combine', '--out', combn] + outs, capture_output=True, text=True) if r.returncode != 0: raise CertError(f'near combine failed: {r.stderr[-800:]}') combt = os.path.join(wdir, 'tail_combine.json') r = subprocess.run([py, cs, 'combine', '--out', combt, '--tails-L', str(L)] + tout, capture_output=True, text=True) if r.returncode != 0: raise CertError(f'tail combine failed: {r.stderr[-800:]}') nj, tj = json.load(open(combn)), json.load(open(combt)) cpu_tot += nj['cpu_total_s'] + tj['cpu_total_s'] with mp.workdps(160): for k in want_div: fresh[k] = (_dec_exact(nj['vals_exact'][k]) + mp.mpf(tj['out'][k]['tails'][str(L)])) print(f' near+tail combines: cpu_total {nj["cpu_total_s"]} + ' f'{tj["cpu_total_s"]} s') # ---- gate vs cache ---- bar = bar_arg if bar_arg is not None else RC_TIERS['conv'][tier]['bar'] okall = True print(f' gate: fresh chain recompute vs sha-pinned cache (bar {bar} d; tier ' f'expectation from MEASURED config crosses: probe ~22.4 d [mm4/mm5], ' f'control ~31.3 d [mm5/mm6 floor], mid/deep ~40-47 d [mo5 systematic ' f'cap, recorded with the cache])') with mp.workdps(160): for k in comps: fam = 'convergent' if k in CONV_KEYS else 'divergent' ent = cb[fam][k] cert = ent['d_2prec'] if fam == 'convergent' else ent['d_cert'] d = digits(fresh[k], mp.mpf(ent['value'])) ok = d >= bar okall = okall and ok print(f' {k:2s}: fresh vs cache {d:7.2f} d (cache cert {cert:.2f} d) ' f'[{"PASS" if ok else "FAIL"}]') print(f' recompute cpu_total {cpu_tot:.0f} s, wall {time.time()-t0:.0f} s') if not okall: raise CertError(f'--boundary-recompute gate FAILED (bar {bar} d) -- the cache ' 'did not reproduce from the definition chain at this tier') print(' boundary recompute: PASS') return 0 # =========================================================================== # closed periods (certified principal branches; periods.py of the campaign) # =========================================================================== def m_mod(l): return (l**2 - 1) / (9*l**2 - 1) def pref(l): return 1 / mp.sqrt(1 - 9*l**2) def D4(l): return (l**2 - 1) * (9*l**2 - 1) def varpi0(l): return mp.ellipk(m_mod(l)) * pref(l) def varpi1(l): return mp.ellipk(1 - m_mod(l)) * pref(l) def varpi0p(l): m = m_mod(l) K = mp.ellipk(m); E = mp.ellipe(m) dK = (E - (1-m)*K) / (2*m*(1-m)) return dK*16*l/(9*l**2-1)**2*pref(l) + K*9*l*pref(l)**3 def point1(lam, dps): print('== Point 1: closed L2 periods at lambda = %s ==' % lam) ok = True for D in (dps, 2*dps): mp.mp.dps = D + 15 l = mp.mpf(lam.numerator) / lam.denominator m = m_mod(l) # live oracle: defining integral of K (tanh-sinh, never calls ellipk). # On (0,1/3) the modulus is m>1: use the reciprocal-modulus decomposition # K(m) = [K(1/m) - i K(1-1/m)]/sqrt(m) (DLMF 19.7.3, principal branch -- # a classical identity, both pieces real quadratures with moduli in (0,1)). # charter v2: refine-until-bound -- accept only when mpmath's own # successive-level estimate beats qtol; escalate maxdegree; raise at cap. qtol = mp.mpf(10) ** (-(D + QUAD_GUARD)) qest = [mp.mpf(0)] def Kquad(mm): v, e = _certified_quad(lambda th: 1/mp.sqrt(1 - mm*mp.sin(th)**2), [0, mp.pi/2], qtol, f'point1 K-integral (dps {D})') qest[0] += e return v if m < 1: Kq = Kquad(m) else: Kq = (Kquad(1/m) - 1j*Kquad(1 - 1/m)) / mp.sqrt(m) d_K = digits(Kq * pref(l), varpi0(l)) # L2 annihilation (high-order FD) h = mp.mpf(10) ** (-(D + 15)//3) def L2(f): d1 = (f(l+h) - f(l-h)) / (2*h) d2 = (f(l+h) - 2*f(l) + f(l-h)) / h**2 c1 = (45*l**4 - 30*l**2 + 1) / (l*D4(l)) c0 = (27*l**2 - 10) / D4(l) return abs(d2 + c1*d1 + c0*f(l)) r0, r1 = L2(varpi0), L2(varpi1) # Wronskian (branch-sensitive) w_num = varpi0(l)*(varpi1(l+h)-varpi1(l-h))/(2*h) - varpi1(l)*(varpi0(l+h)-varpi0(l-h))/(2*h) d_W = digits(w_num, mp.pi/(2*l*D4(l))) fd_floor = (D + 15) // 3 # FD truncation floor dks = f'{d_K:8.1f}' if d_K != float('inf') else f' >={D+15}' print(f' dps={D}: varpi0 vs K-integral oracle: {dks} d ' f'|L2 varpi0|={mp.nstr(r0,3)} |L2 varpi1|={mp.nstr(r1,3)} Wronskian: {d_W:.1f} d (FD floor ~{fd_floor})') if D == dps: print(f' varpi0 = {mp.nstr(varpi0(l), D)}') print(f' varpi1 = {mp.nstr(varpi1(l), D)}') ok = ok and d_K > min(D, 30) and d_W > 28 # promoted raising gates (charter v2; measured headroom >=1e6, see # ROW_REPORT.md calibration: healthy d_K ~ D+15, d_W ~ 2*fd_floor, # |L2| ~ 10^-(D+15-2*fd_floor-1)) l2bar = mp.mpf(10) ** (-(D + 15) + 2*fd_floor + 8) if not d_K >= D + 4: raise CertError(f'point1: K-integral oracle agreement {d_K:.1f} d < raise-bar ' f'{D+4} at dps {D} (quad est {mp.nstr(qest[0],3)}, tol {mp.nstr(qtol,3)})') if not d_W >= 2*fd_floor - 8: raise CertError(f'point1: Wronskian agreement {d_W:.1f} d < raise-bar ' f'{2*fd_floor-8} at dps {D}') if not (r0 < l2bar and r1 < l2bar): raise CertError(f'point1: L2-annihilation residuals |L2 varpi0|={mp.nstr(r0,3)}, ' f'|L2 varpi1|={mp.nstr(r1,3)} above raise-bar {mp.nstr(l2bar,3)} at dps {D}') print(f' [certified] dps={D}: K-quad est {mp.nstr(qest[0],2)} <= tol {mp.nstr(qtol,2)}; ' f'raise-bars armed: d_K>={D+4}, d_W>={2*fd_floor-8}, |L2|<{mp.nstr(l2bar,2)}') print(' Point 1: ' + ('PASS (oracle digits scale with dps)' if ok else 'FAIL')) # NOTE: the promoted RAISES above use the dps-SCALED calibrated bars (they # fire on genuine damage at any dps, never on a healthy run). The legacy # display bar `d_W > 28` does NOT scale down and sits right at the healthy # FD floor for dps <~ 32 (measured: d_W = 2*fd_floor +- 2, e.g. 27.9 at # dps 30, lam=2/7) -- a pre-existing floor of the original evaluator, so # its FAIL keeps the original quiet exit-code semantics (FLAG for owner). return ok # =========================================================================== # Point 2: exact A16 checks + live sqrt(Q4)-twist residue gate # =========================================================================== def _eps_rat(c, ev): n = sum(Fr(v) * ev**j for j, v in enumerate(c['n'])) d = sum(Fr(v) * ev**j for j, v in enumerate(c['d'])) return n / d def point2(A16, epsv, dps): print('== Point 2: symbolic-eps connection A16 (exact rationals) ==') A = A16['A16'] # --- (i) EXACT: A16[z5,z5] == eps * dlog D4 (Fraction grid = rigorous) --- e = A['13,13'] degx = max(len(e['N']), len(e['D'])) + 4 dege = max(max(len(c['n']), len(c['d'])) for c in e['N'] + e['D']) + 2 D4c = [Fr(1), Fr(0), Fr(-10), Fr(0), Fr(9)] # 9x^4-10x^2+1 (ascending) D4pc = [Fr(0), Fr(-20), Fr(0), Fr(36)] # 36x^3-20x def pev(cs, x): return sum(c * x**i for i, c in enumerate(cs)) exact_ok = True for ie in range(dege): ev = Fr(ie + 1, ie + 7) for ix in range(degx): xv = Fr(ix + 2, 2*ix + 5) lhs = pev([_eps_rat(c, ev) for c in e['N']], xv) * pev(D4c, xv) rhs = ev * pev(D4pc, xv) * pev([_eps_rat(c, ev) for c in e['D']], xv) if lhs != rhs: exact_ok = False print(f' exact identity A16[z5,z5] == eps*dlog(D4): ' f'{"PROVED (exact on full bidegree grid)" if exact_ok else "FAIL"}') # --- (ii) live twist-residue gate at the roots of Q4 --- mp.mp.dps = dps + 20 roots = mp.polyroots([20, 39, 29, 9, 1], maxsteps=200, extraprec=dps + 60) target = mp.mpf(epsv.numerator)/epsv.denominator - mp.mpf(1)/2 worst = float('inf'); nrow = set() for rho in roots: # residue matrix: entries with D(rho)=0 (simple zero of Q4 | D) for key, ent in A.items(): i, j = map(int, key.split(',')) Dc = [_eps_rat(c, epsv) for c in ent['D']] Dv = sum(mp.mpf(c.numerator)/c.denominator * rho**k for k, c in enumerate(Dc)) Dn = max(abs(mp.mpf(c.numerator)/c.denominator) for c in Dc) if abs(Dv) < Dn * mp.mpf(10)**(-dps): # pole at rho nrow.add(i) if i == j: Nv = sum(mp.mpf(c.numerator)/c.denominator * rho**k for k, c in enumerate([_eps_rat(c, epsv) for c in ent['N']])) Dp = sum(k * mp.mpf(c.numerator)/c.denominator * rho**(k-1) for k, c in enumerate(Dc) if k) worst = min(worst, digits(Nv/Dp, target)) ok = exact_ok and worst > min(dps - 8, 30) and nrow == {5} print(f' sqrt(Q4) twist: poles only in row 5 (e6 direction): {sorted(nrow)}; ' f'rank-1 residue => spectrum {{eps-1/2, 0^15}}') print(f' residue eigenvalue vs eps-1/2 at eps={epsv}, all 4 live Q4 roots: ' f'min {worst:.1f} d (bar {min(dps-8,30)})') print(' Point 2: ' + ('PASS' if ok else 'FAIL')) if not ok: # promoted to a raise (charter v2) raise CertError(f'point2: FAILED (exact identity {exact_ok}; twist-residue min ' f'{worst:.1f} d vs bar {min(dps-8,30)}; pole rows {sorted(nrow)} vs [5])') return ok # =========================================================================== # Point 3: boundary-free per-eps-order holonomy (prediction vs live oracle) # =========================================================================== def point3(B, dps, orders, mutate=False, quiet=False): from flint import acb, arb, acb_mat, arb_mat, acb_series, arb_series, ctx import math LET = {l['id']: l for l in B['letters']} EICH = B['eichler_kernels'] WORDS = B['words'] IDX = list(range(5, 14)) # SUB9 rows/cols of A16 A = B['A16']['A16'] # ---------- shared exact-parse helpers ---------- xs = sp.Symbol('x') _rc = {} def rcoeffs(lid): """letter dressing r(x) -> (int num coeffs, int den coeffs), exact.""" if lid not in _rc: num, den = sp.fraction(sp.cancel(sp.sympify(LET[lid]['r']))) pn, pd = sp.Poly(num, xs), sp.Poly(den, xs) M = 1 for c in pn.coeffs() + pd.coeffs(): M = lcm(M, sp.Rational(c).q) _rc[lid] = ([int(sp.Rational(c)*M) for c in pn.all_coeffs()][::-1], [int(sp.Rational(c)*M) for c in pd.all_coeffs()][::-1]) return _rc[lid] # ---------- live certified singularity set (letters ∪ D4 ∪ Q4) ---------- # Certified convergence/decay ratios for BOTH routes come from the distance # to the nearest singularity of the integrands: the union of all letter- # denominator roots, the D4 roots (+-1, +-1/3) and the Q4 roots, all found # LIVE (measured 2026-07-05: nearest to the path is x=1/3 at distance 1/6; # the Q4 roots are >= 0.79 away). def _roots_of(desc): key = tuple(desc) if key not in _ROOTS_CACHE: cs = list(desc) while cs and cs[0] == 0: cs = cs[1:] if len(cs) <= 1: _ROOTS_CACHE[key] = [] else: with mp.workdps(60): _ROOTS_CACHE[key] = [complex(r) for r in mp.polyroots(cs, maxsteps=2000, extraprec=400)] return _ROOTS_CACHE[key] SINGS = [] for lid in LET: SINGS += _roots_of(rcoeffs(lid)[1][::-1]) SINGS += _roots_of([9, 0, -10, 0, 1]) # D4 (descending) SINGS += _roots_of([20, 39, 29, 9, 1]) # Q4 (descending) # ======================================================================= # PREDICTION route: word lists as iterated integrals (Chebyshev spectral) # charter v2: N below is the STARTING seed of a refine-until-bound loop; # every cumint carries a certified trailing-window tail bound (last-8 # Chebyshev coefficients x r/(1-r), r = certified Bernstein ratio) which # propagates through the nested iterated integrals to a per-entry bound; # accept only if worst bound < 10^-(dps+PRED_GUARD), else N *= 1.5 (exact # recomputation, same rules) up to 8*N0, then CertError. # ======================================================================= gdps = dps + 15 N0 = max(64, int((gdps + 170) / 3.6) + 1) # calibrated seed: N=64->~65d, N=72->~94d NP = 10 PATH_LEN = mp.mpf(1) / 5 # |7/10 - 1/2|, exact tol_pred = mp.mpf(10) ** (-(dps + PRED_GUARD)) def pred_route(N): ctx.prec = int(gdps * 3.33) + 30 PI = acb.pi() tt = [-acb.cos(PI*k/N) for k in range(N+1)] a, b = acb(1)/2, acb(7)/10 edges = [a + (b-a)*i/NP for i in range(NP+1)] XS = []; PS = [] for p in range(NP): lo, hi = edges[p], edges[p+1]; mid, half = (lo+hi)/2, (hi-lo)/2 st = len(XS); XS.extend([mid + half*tk for tk in tt]); PS.append((st, st+N+1, half)) NX = len(XS) V = acb_mat(N+1, N+1) for j in range(N+1): for k in range(N+1): w = acb.cos(PI*j*(N-k)/N); fac = acb(2)/N if k == 0 or k == N: fac /= 2 if j == 0 or j == N: fac /= 2 V[j, k] = fac*w M1 = acb_mat(N+1, N+1) for j in range(1, N+2): M1[j-1, j-1] = (acb(2) if j == 1 else acb(1))/(2*j) if j+1 <= N: M1[j-1, j+1] = -acb(1)/(2*j) M2 = acb_mat(N+1, N+1) # T_j(t_k)=(-1)^j cos(j pi k/N): endpoint-stable (no acos) for k in range(N+1): for j in range(1, N+2): M2[k, j-1] = ((-1)**j) * (acb.cos(PI*j*k/N) - 1) OP = M2 * (M1 * V) # certified per-panel Chebyshev decay ratios from the live singularity set rpan = [_bernstein_r(float(edges[p].real.mid()), float(edges[p+1].real.mid()), SINGS) for p in range(NP)] rmax = max(rpan) def cumint(fv): out = [None]*NX; carry = acb(0) qerr = mp.mpf(0) for ip, (s, e, half) in enumerate(PS): col = acb_mat([[fv[i]] for i in range(s, e)]) av = OP * col cf = V * col # Chebyshev coeffs (certification only) tw = mp.mpf(0) for jj in range(max(0, N - 7), N + 1): u = _absu(cf[jj, 0]) if u > tw: tw = u rr = rpan[ip] qerr += 2 * _absu(half) * tw * (rr / (1 - rr)) for i in range(e - s): out[s+i] = carry + av[i, 0]*half carry = out[e-1] return out, qerr W0g = [acb.elliptic_k((x**2-1)/(9*x**2-1)) / (1-9*x**2).sqrt() for x in XS] def w0df(x): m = (x**2-1)/(9*x**2-1) K = acb.elliptic_k(m); E = acb.elliptic_e(m) dK = (E - (1-m)*K)/(2*m*(1-m)) pf = 1/(1-9*x**2).sqrt() return dK*16*x/(9*x**2-1)**2*pf + K*9*x*pf**3 W0Dg = [w0df(x) for x in XS] SQg = [(((20*x+39)*x+29)*x*x + 9*x + 1).sqrt() for x in XS] def horner(cs, x): s = acb(0) for c in reversed(cs): s = s*x + c return s NG = {} NGERR = {} _lg = {}; _lgerr = {}; _lgmax = {} def _lbuild(lid, skip=None): """grid values of letter lid; skip=k drops ONE power of n_{k+1} (certification exclusion product only -- never on the value path).""" L = LET[lid]; cn, cd = rcoeffs(lid) out = [] for i, x in enumerate(XS): v = horner(cn, x)/horner(cd, x) if L['w0_pow']: v *= W0g[i]**L['w0_pow'] if L['w0d_pow']: v *= W0Dg[i]**L['w0d_pow'] if L['sqQ4_pow']: v *= SQg[i]**L['sqQ4_pow'] if L['pi_pow']: v *= PI**L['pi_pow'] for k, dk in enumerate(L['n_dress']): dk2 = dk - (1 if skip == k else 0) if dk2: v *= NG[k+1][i]**dk2 out.append(v) return out def lg(lid): if lid not in _lg: _lg[lid] = _lbuild(lid) return _lg[lid] def lgmax(lid): if lid not in _lgmax: _lgmax[lid] = max(_absuf(v) for v in lg(lid)) return _lgmax[lid] def lgerr(lid): if lid not in _lgerr: L = LET[lid] err = mp.mpf(0) for k, dk in enumerate(L['n_dress']): if not dk: continue ex = _lbuild(lid, skip=k) # one power of n_{k+1} removed err += dk * NGERR[k+1] * max(_absuf(v) for v in ex) _lgerr[lid] = err return _lgerr[lid] for k in (1, 2, 3, 4, 5): # eich_5 needs n1: order is significant s = [acb(0)]*NX esum = mp.mpf(0) for lid in EICH[f'eich_{k}']: g = lg(lid) esum += lgerr(lid) s = [u+v for u, v in zip(s, g)] NG[k], qe = cumint(s) NGERR[k] = qe + PATH_LEN * esum _suf = {}; _suferr = {}; _sufmax = {} def word_val(letters): key = tuple(letters) if key in _suf: return _suf[key] if len(letters) == 1: g, qe = cumint(lg(letters[0])) err = qe + PATH_LEN * lgerr(letters[0]) else: ikey = tuple(letters[1:]) inner = word_val(letters[1:]) g, qe = cumint([u*v for u, v in zip(lg(letters[0]), inner)]) err = qe + PATH_LEN * (mp.mpf(lgmax(letters[0])) * _suferr[ikey] + lgerr(letters[0]) * _sufmax[ikey]) _suf[key] = g _suferr[key] = err _sufmax[key] = mp.mpf(max(_absuf(v) for v in g)) return g PRED = {} PERR = {} mut_tag = None predworst = mp.mpf(0); predarg = None for o in range(orders + 1): t1 = time.time() ent = {} for key, v in WORDS[f'eps^{o}'].items(): tot = acb(v['delta']) esum = mp.mpf(0) for iw, wd in enumerate(v['words']): wv = word_val(wd['letters'])[-1] esum += _suferr[tuple(wd['letters'])] if mutate and mut_tag is None: wv = wv * (1 + acb(10)**(-30)) # mutation control: one word only mut_tag = (o, key, iw) tot = tot + wv ent[key] = tot ee = esum + _mpf2(tot.real.rad().man_exp()) + _mpf2(tot.imag.rad().man_exp()) PERR[(o, key)] = ee if ee > predworst: predworst, predarg = ee, (o, key) PRED[o] = ent if not quiet: print(f' prediction eps^{o}: {len(ent)} entries, {len(_suf)} cached integrals, {time.time()-t1:.1f}s') if predworst >= tol_pred: # refine-until-bound gate raise _Escalate(f'eps^{o} entry {predarg[1]} certified bound ' f'{mp.nstr(predworst, 3)} >= tol {mp.nstr(tol_pred, 3)}') return PRED, PERR, predworst, predarg, mut_tag, rmax N = N0 while True: try: PRED, PERR, predworst, predarg, mut_tag, rmax = pred_route(N) break except _Escalate as esc: Nn = int(N * 1.5) + 1 if Nn > 8 * N0: raise CertError(f'point3 prediction route: refine-until-bound FAILED at cap: ' f'{esc} (N={N}, seed {N0}, cap {8*N0}, tol {mp.nstr(tol_pred,3)})') if not quiet: print(f' [escalate] prediction Chebyshev degree N {N} -> {Nn}: {esc}') N = Nn if mutate and not quiet: print(f' [mutation control armed: word {mut_tag} scaled by 1+1e-30]') maxrad = max(max(v.real.rad(), v.imag.rad()) for E in PRED.values() for v in E.values()) if not quiet: print(f' prediction ball radius (arithmetic error, excl. spectral truncation): {float(maxrad):.1e}') print(f' [certified] prediction route: worst entry bound {mp.nstr(predworst, 2)} <= tol ' f'{mp.nstr(tol_pred, 2)} (entry {predarg}, N={N}, seed {N0}, cap {8*N0}, ' f'panel r_max {rmax:.3f}; tail+ball, fail-closed)') # ======================================================================= # REFERENCE route (independent): arb Taylor ODE of the ORIGINAL SUB9 system # ======================================================================= fitdps = dps + 90 # nodes per sign: eps-Vandermonde truncation ~ (M*h)^(2M-k) must beat the # target working floor (measured: M=8 caps eps^2 at ~52 d) # charter v2: M and NT below are STARTING seeds -- per-step Taylor tails # and the value-level eps-Vandermonde bound are certified at runtime, # with escalation (NT x1.5 exact continuation; M += 2 extra eps nodes) # and CertError at the caps. M = orders + 7 + max(0, (dps - 32) // 8) H = Fr(1, 100000) ctx.prec = int(fitdps * 3.33) + 20 RATIO = 0.18 NT = int(fitdps * 2.303 / (-math.log(RATIO))) + 12 ctx.cap = NT + 2 tol_ref = mp.mpf(10) ** (-(dps + REF_GUARD)) # value level, per eps-order tol_step = mp.mpf(10) ** (-(fitdps - REF_STEP_SLACK)) # per ODE-step tail def entry_at_eps(ent, ev): cn = [_eps_rat(c, ev) for c in ent['N']] cd = [_eps_rat(c, ev) for c in ent['D']] Mm = 1 for c in cn + cd: Mm = lcm(Mm, c.denominator) return [int(c*Mm) for c in cn], [int(c*Mm) for c in cd] def shift_series(coeffs, c): ct = arb_series([c, 1]); s = arb_series([0]) for a2 in reversed(coeffs): s = s*ct + a2 return s ONE3 = arb(1)/3 def radius(c): return min(abs(c - ONE3), abs(1 - c), arb('0.35')) def solve_phi(ev): ENT = {} for a2, j in enumerate(IDX): for b2, k in enumerate(IDX): e = A.get(f'{j},{k}') if e is not None: ENT[(a2, b2)] = entry_at_eps(e, ev) # certified step ratios: denominator roots of THIS eps's SUB9 entries evroots = list(SINGS) for (cn, cd) in ENT.values(): evroots += _roots_of(list(cd)[::-1]) Phi = arb_mat(9, 9) for i in range(9): Phi[i, i] = arb(1) errtot = mp.mpf(0) xc = arb(1)/2; target = arb(7)/10 while float(xc) < float(target) - 1e-18: R = radius(xc); h = arb(RATIO)*R if float(xc + h) > float(target): h = target - xc xf = float(xc) dmin = min(abs(complex(xf) - z) for z in evroots) hu = abs(_mpf2(h.mid().man_exp())) + _mpf2(h.rad().man_exp()) rstep = hu / (dmin * (1 - 1e-9)) if not rstep < mp.mpf('0.6'): raise CertError(f'point3 reference ODE: step ratio h/dmin = {mp.nstr(rstep,3)} ' f'not < 0.6 at x0={xf:.6f} (h={mp.nstr(hu,3)}, dmin={dmin:.4f})') NTs = NT while True: # refine-until-bound (charter v2) Am = [arb_mat(9, 9) for _ in range(NTs)] for (i, j), (cn, cd) in ENT.items(): if all(c == 0 for c in cn): continue ser = shift_series(cn, xc) / shift_series(cd, xc) cl = list(ser.coeffs()) for m2 in range(min(NTs, len(cl))): Am[m2][i, j] = cl[m2] phi = [Phi] for m2 in range(NTs): Mx = arb_mat(9, 9) for r in range(m2 + 1): Mx = Mx + Am[r]*phi[m2-r] phi.append(Mx * (arb(1)/(m2+1))) # certified trailing-window tail (pilot construction: last-8 # |phi_m h^m| max-entry terms x r/(1-r), r certified above) tail = mp.mpf(0) hp = hu ** max(0, NTs - 7) for m2 in range(max(0, NTs - 7), NTs + 1): t = _mat_maxu(phi[m2]) * hp if t > tail: tail = t hp = hp * hu tail = tail * (rstep / (1 - rstep)) if tail < tol_step: break NTn = int(NTs * 1.5) + 1 if NTn > 8 * NT: ctx.cap = NT + 2 raise CertError(f'point3 reference ODE step at x0={xf:.6f} (eps={ev}): ' f'certified tail bound {mp.nstr(tail,3)} >= tol ' f'{mp.nstr(tol_step,3)} at N={NTs} (seed {NT}, cap {8*NT}, ' f'h={mp.nstr(hu,3)}, r={mp.nstr(rstep,3)})') if not quiet: print(f' [escalate] reference ODE step x0={xf:.6f}: Taylor depth {NTs} -> ' f'{NTn} (tail {mp.nstr(tail,3)} >= tol {mp.nstr(tol_step,3)})') NTs = NTn ctx.cap = NTs + 2 errtot += tail S = arb_mat(9, 9) for m2 in range(NTs, -1, -1): S = S*h + phi[m2] Phi = S; xc = xc + h if ctx.cap != NT + 2: ctx.cap = NT + 2 radm = mp.mpf(0) for i in range(9): for j in range(9): v = _mpf2(Phi[i, j].rad().man_exp()) if v > radm: radm = v return Phi, errtot, radm # ---- Eichler n_k(7/10) by an INDEPENDENT acb-series Taylor stepper ---- # (w0 series from the L2 recurrence seeded by the closed form per step; # not the spectral route used by the prediction side) def nk_taylor(ndps): ctx.prec = int(ndps * 3.33) + 20 NT2 = int(ndps * 2.303 / (-math.log(RATIO))) + 12 ctx.cap = max(ctx.cap, NT2 + 2) cap0 = ctx.cap tol_nk = mp.mpf(10) ** (-(ndps - REF_STEP_SLACK)) PI2 = acb.pi() # at THIS precision (the outer PI is grid-precision) mp.mp.dps = ndps + 10 eich_data = {k: [(LET[lid], rcoeffs(lid)) for lid in EICH[f'eich_{k}']] for k in (1,2,3,4,5)} nvals = {k: acb(0) for k in (1,2,3,4,5)} errnk = mp.mpf(0) xc = Fr(1, 2); target = Fr(7, 10) def a_from_mp(z): return acb(mp.nstr(z.real, ndps+8), mp.nstr(z.imag, ndps+8)) while xc < target: R = min(abs(float(xc) - 1/3), abs(1 - float(xc)), 0.35) h = Fr(int(RATIO*R*10**6), 10**6) if xc + h > target: h = target - xc xa = acb(xc.numerator)/xc.denominator lm = mp.mpf(xc.numerator)/xc.denominator dmin = min(abs(complex(float(xc)) - z) for z in SINGS) hu = mp.mpf(h.numerator) / h.denominator rstep = hu / (dmin * (1 - 1e-9)) NTs = NT2 while True: # refine-until-bound (charter v2) # w0 series by L2 recurrence: w0'' = -(c1 w0' + c0 w0) num_c1 = shift_acb([1, 0, -30, 0, 45], xa, NTs) # 45x^4-30x^2+1 # exact polys: x*D4 = 9x^5-10x^3+x ; D4 = 9x^4-10x^2+1 xD4 = shift_acb([0, 1, 0, -10, 0, 9], xa, NTs) D4s = shift_acb([1, 0, -10, 0, 9], xa, NTs) c1s = poly_div(num_c1, xD4, NTs) c0s = poly_div(shift_acb([-10, 0, 27], xa, NTs), D4s, NTs) w0s = [a_from_mp(varpi0(lm)), a_from_mp(varpi0p(lm))] for m2 in range(NTs - 2): s1 = sum(c1s[r]*(m2+1-r)*w0s[m2+1-r] for r in range(m2+2) if r < len(c1s)) s0 = sum(c0s[r]*w0s[m2-r] for r in range(m2+1) if r < len(c0s)) w0s.append(-(s1 + s0)/((m2+2)*(m2+1))) w0ds = [ (m2+1)*w0s[m2+1] for m2 in range(NTs-1) ] + [acb(0)] Q4s = shift_acb([1, 9, 29, 39, 20], xa, NTs) SQs = series_sqrt(Q4s, NTs) nser = {} newvals = {} wtail = mp.mpf(0) for k in (1, 2, 3, 4, 5): es = [acb(0)]*NTs for L, (cn, cd) in eich_data[k]: rs = poly_div(shift_acb(cn, xa, NTs), shift_acb(cd, xa, NTs), NTs) v = rs for _ in range(L['w0_pow']): v = series_mul(v, w0s, NTs) for _ in range(L['w0d_pow']): v = series_mul(v, w0ds, NTs) for _ in range(L['sqQ4_pow']): v = series_mul(v, SQs, NTs) if L['pi_pow']: v = [c * PI2**L['pi_pow'] for c in v] for kd, dk in enumerate(L['n_dress']): if dk: assert kd + 1 in nser, f'n{kd+1} dressing before it is built' for _ in range(dk): v = series_mul(v, nser[kd+1], NTs) es = [u+w for u, w in zip(es, v)] nks = [nvals[k]] + [es[m2]/(m2+1) for m2 in range(NTs-1)] # antiderivative nser[k] = nks ha = acb(h.numerator)/h.denominator newvals[k] = series_eval(nks, ha) # certified trailing-window tail of the evaluated series tl = mp.mpf(0) hp = hu ** max(0, len(nks) - 8) for m2 in range(max(0, len(nks) - 8), len(nks)): t = _absu(nks[m2]) * hp if t > tl: tl = t hp = hp * hu tl = tl * (rstep / (1 - rstep)) if tl > wtail: wtail = tl if wtail < tol_nk: break NTn = int(NTs * 1.5) + 1 if NTn > 8 * NT2: ctx.cap = cap0 raise CertError(f'point3 nk_taylor step at x0={float(xc):.6f}: certified tail ' f'{mp.nstr(wtail,3)} >= tol {mp.nstr(tol_nk,3)} at N={NTs} ' f'(seed {NT2}, cap {8*NT2}, h={mp.nstr(hu,3)}, r={mp.nstr(rstep,3)})') if not quiet: print(f' [escalate] nk_taylor step x0={float(xc):.6f}: Taylor depth {NTs} -> ' f'{NTn} (tail {mp.nstr(wtail,3)} >= tol {mp.nstr(tol_nk,3)})') NTs = NTn ctx.cap = max(cap0, NTs + 2) for k in (1, 2, 3, 4, 5): nvals[k] = newvals[k] errnk += 5 * wtail + mp.mpf(10) ** (-(ndps + 6)) # 5 kernels + seed-string slop xc = xc + h if ctx.cap != cap0: ctx.cap = cap0 radnk = mp.mpf(0) for k in (1, 2, 3, 4, 5): v = _mpf2(nvals[k].real.rad().man_exp()) + _mpf2(nvals[k].imag.rad().man_exp()) if v > radnk: radnk = v return [nvals[k] for k in (1, 2, 3, 4, 5)], errnk + radnk def shift_acb(int_coeffs, xa, NT2): """poly with int coeffs (ascending) recentred at xa, as coeff list len<=NT2.""" ser = acb_series([0]) ct = acb_series([xa, 1]) for c in reversed(int_coeffs): ser = ser*ct + c cl = list(ser.coeffs()) cl = cl + [acb(0)]*(NT2 - len(cl)) return cl[:NT2] def series_mul(u, v, NT2): out = [acb(0)]*NT2 for i, ui in enumerate(u): if i >= NT2: break for j2, vj in enumerate(v): if i + j2 >= NT2: break out[i+j2] += ui*vj return out def poly_div(u, v, NT2): inv = [1/v[0]] for m2 in range(1, NT2): s = sum(v[r]*inv[m2-r] for r in range(1, min(m2, len(v)-1)+1)) inv.append(-s/v[0]) return series_mul(u, inv, NT2) def series_sqrt(u, NT2): s0 = u[0].sqrt() out = [s0] for m2 in range(1, NT2): s = sum(out[r]*out[m2-r] for r in range(1, m2)) out.append((u[m2] - s)/(2*s0)) return out def series_eval(cs, h): s = acb(0) for c in reversed(cs): s = s*h + c return s t1 = time.time() nend, nk_err = nk_taylor(fitdps) if not quiet: print(f' reference n_k(7/10) (acb-series Taylor stepper): {time.time()-t1:.1f}s, ' f'max rad {max(float(max(z.real.rad(), z.imag.rad())) for z in nend):.1e}') t1 = time.time() Mcur = M ms = [m2 for m2 in range(-Mcur, Mcur+1) if m2 != 0] phis = {}; oderr = {}; oderad = {} for m2 in ms: phis[m2], oderr[m2], oderad[m2] = solve_phi(m2 * H) if not quiet: print(f' reference ODE (stored A16 SUB9, {len(ms)} eps nodes, {fitdps}d): {time.time()-t1:.1f}s') # ---- T9 rotation + eps-Vandermonde extraction (mpmath at fitdps) ---- t1 = time.time() mp.mp.dps = fitdps SY = sp.symbols('x eps w0 w0d SQ PI n1 n2 n3 n4 n5') def lamb_mat(rows): return [[(None if s in ('0', 0) else sp.lambdify(SY, sp.sympify(s), modules='mpmath')) for s in row] for row in rows] LU1 = lamb_mat(B['gauge']['U1']); LTg = lamb_mat(B['gauge']['Tg']); LTzc = lamb_mat(B['gauge']['Tzc']) def acb2mpc(z, nd): return mp.mpc(mp.mpf(z.real.mid().str(nd, radius=False)), mp.mpf(z.imag.mid().str(nd, radius=False))) nvals_end = [acb2mpc(z, fitdps + 5) for z in nend] def t9_at(l, ev, nv): av = (l, ev, varpi0(l), varpi0p(l), mp.sqrt(((20*l+39)*l+29)*l*l + 9*l + 1), mp.pi) + tuple(nv) def em(L): return mp.matrix([[mp.mpc(0) if L[i][j] is None else L[i][j](*av) for j in range(4)] for i in range(4)]) Tc = em(LU1)*em(LTg); Tz = em(LTzc) T9 = mp.matrix(9, 9) for i in range(4): for j in range(4): T9[i, j] = Tc[i, j]; T9[4+i, 4+j] = Tz[i, j] T9[8, 8] = mp.mpf(1) return T9 l_end, l_half = mp.mpf(7)/10, mp.mpf(1)/2 # measured n_k sensitivity of the endpoint rotation (x10 safety margin; # gauge entries are polynomial in n_k, so an FD probe bounds the slope) _d = mp.mpf(10) ** (-(fitdps // 2)) _ev0 = mp.mpf(H.numerator) / H.denominator _Te0 = t9_at(l_end, _ev0, nvals_end) _Tep = t9_at(l_end, _ev0, [nv + _d for nv in nvals_end]) nk_sens = 10 * _mnorm(_Te0 - _Tep) / _d Hf = mp.mpf(H.numerator) / H.denominator while True: # refine-until-bound on M (charter v2) Psis = [] errPsi = [] for m2 in ms: ev = mp.mpf(m2*H.numerator)/H.denominator Phi = mp.matrix(9, 9) for i in range(9): for j in range(9): s = phis[m2][i, j].mid().str(fitdps + 5, radius=False) Phi[i, j] = mp.mpf(s) Te = t9_at(l_end, ev, nvals_end) Th = t9_at(l_half, ev, [mp.mpf(0)]*5) Thi = Th**-1 PT = Phi * Thi Psis.append(Te * PT) # value-level error budget for this node: ODE step tails + Phi ball # radii through the exact rotation norms, nk error through the # measured sensitivity, + the fitdps+5-digit string conversion. errPsi.append(_mnorm(Te) * _mnorm(Thi) * 9 * (oderr[m2] + oderad[m2]) + nk_sens * nk_err * _mnorm(PT) + mp.mpf(10) ** (-(fitdps - 2))) # Vandermonde solve for eps-Taylor orders 0..len(ms)-1 (need 0..orders) NN = len(ms) Vm = mp.matrix(NN, NN) for r, m2 in enumerate(ms): ev = mp.mpf(m2*H.numerator)/H.denominator for c in range(NN): Vm[r, c] = ev**c REF = {o: mp.matrix(9, 9) for o in range(orders + 1)} chat = mp.mpf(0) for i in range(9): for j in range(9): rhs = mp.matrix([Psis[r][i, j] for r in range(NN)]) co = mp.lu_solve(Vm, rhs) for o in range(orders + 1): REF[o][i, j] = co[o] for p in (NN - 3, NN - 2, NN - 1): if abs(co[p]) > chat: chat = abs(co[p]) Wi = Vm**-1 amp = [sum(abs(Wi[o, r2]) for r2 in range(NN)) for o in range(orders + 1)] maxePsi = max(errPsi) etail = [chat * sum(abs(Wi[o, r2]) * (abs(ms[r2]) * Hf)**NN for r2 in range(NN)) for o in range(orders + 1)] refbound = [amp[o] * maxePsi + etail[o] for o in range(orders + 1)] if max(refbound) < tol_ref: break if Mcur >= M + 8: raise CertError(f'point3 reference: value-level bound {mp.nstr(max(refbound),3)} >= ' f'tol {mp.nstr(tol_ref,3)} at M={Mcur} (seed {M}, cap {M+8}); ' f'components: ODE {mp.nstr(max(oderr.values()),3)}, ball ' f'{mp.nstr(max(oderad.values()),3)}, nk {mp.nstr(nk_err,3)}, ' f'eps-tail {mp.nstr(max(etail),3)}') if not quiet: print(f' [escalate] reference eps-Vandermonde: M {Mcur} -> {Mcur+2} ' f'(bound {mp.nstr(max(refbound),3)} >= tol {mp.nstr(tol_ref,3)})') for mnew in (Mcur+1, Mcur+2, -(Mcur+1), -(Mcur+2)): # exact continuation phis[mnew], oderr[mnew], oderad[mnew] = solve_phi(mnew * H) Mcur += 2 ms = [m2 for m2 in range(-Mcur, Mcur+1) if m2 != 0] if not quiet: print(f' T9 rotation + Vandermonde extraction: {time.time()-t1:.1f}s') print(f' [certified] reference route: value-level bounds ' + '/'.join(mp.nstr(refbound[o], 2) for o in range(orders + 1)) + f' <= tol {mp.nstr(tol_ref, 2)} (ODE tails {mp.nstr(max(oderr.values()), 2)}, ' f'ball {mp.nstr(max(oderad.values()), 2)}, nk {mp.nstr(nk_err, 2)}, ' f'eps-tail {mp.nstr(max(etail), 2)}, amp {mp.nstr(max(amp), 2)}, M={Mcur})') # ---- gate ---- mp.mp.dps = dps + 40 results = {} predmp = {} allpass = True for o in range(orders + 1): worst = float('inf'); arg = None; zbad = 0 seen = set() for key, v in PRED[o].items(): i, j = map(int, key.split(',')) seen.add((i, j)) pv = acb2mpc(v, dps + 35) predmp[(o, key)] = pv rv = REF[o][i, j] d = digits(pv, rv) if (pv != 0 or rv != 0) else mp.inf if d < worst: worst, arg = d, (i, j) for i in range(9): for j in range(9): if (i, j) in seen: continue if abs(REF[o][i, j]) > mp.mpf(10)**(-(dps - 10)): zbad += 1 results[o] = (float(worst), arg, zbad) bar = 30 ok = worst > bar and zbad == 0 allpass = allpass and ok if not quiet: print(f' eps^{o}: min agreed digits {float(worst):8.1f} (argmin entry {arg}, ' f'{len(seen)} gated, zero-mismatches {zbad}) ' f'[{ "PASS" if ok else "FAIL" } bar {bar}]') if not quiet and 1 in PRED and '0,1' in PRED[1]: sv = acb2mpc(PRED[1]['0,1'], dps + 5) print(f' sample word-sum value Psi_1[Js0,Js1](7/10) = {mp.nstr(sv, dps)}') return allpass, results, {'pred': predmp, 'nk': nvals_end, 'pred_bound': predworst, 'ref_bound': refbound} # =========================================================================== def _nk_quad_check(B, dps, nk_ref): """--check (iii): n_1..n_4 cross-checked LIVE against direct certified mp.quad of the Eichler kernels (eich_1..eich_4 are n-undressed => closed integrands over the exact letters; independent of BOTH point-3 routes' machinery: mpmath ellipk/ellipe vs the acb spectral and acb-series routes). Fail-closed: raises CertError below the calibrated bar.""" LETm = {l['id']: l for l in B['letters']} xs = sp.Symbol('x') worst = float('inf') with mp.workdps(dps + 15): qtol = mp.mpf(10) ** (-(dps + QUAD_GUARD)) for k in (1, 2, 3, 4): fs = [] for lid in B['eichler_kernels'][f'eich_{k}']: L = LETm[lid] assert not any(L['n_dress']), f'eich_{k} unexpectedly n-dressed' num, den = sp.fraction(sp.cancel(sp.sympify(L['r']))) fs.append((L, sp.lambdify(xs, num/den, modules='mpmath'))) def integrand(x): tot = mp.mpf(0) for L, fr in fs: v = fr(x) if L['w0_pow']: v *= varpi0(x)**L['w0_pow'] if L['w0d_pow']: v *= varpi0p(x)**L['w0d_pow'] if L['sqQ4_pow']: v *= mp.sqrt(((20*x+39)*x+29)*x*x + 9*x + 1)**L['sqQ4_pow'] if L['pi_pow']: v *= mp.pi**L['pi_pow'] tot += v return tot v, est = _certified_quad(integrand, [mp.mpf(1)/2, mp.mpf(7)/10], qtol, f'--check n_{k} quadrature') d = digits(v, nk_ref[k-1]) worst = min(worst, d) print(f' n_{k}(7/10) vs certified mp.quad: {d:.1f} d (quad est {mp.nstr(est,2)})') bar = dps - 2 if not worst >= bar: raise CertError(f'--check n_k cross-check: min agreement {worst:.1f} d < bar {bar}') print(f' n_k cross-check: PASS (min {worst:.1f} d >= bar {bar})') return worst # =========================================================================== # =========================================================================== # --headline (Point 5): recompute the page's displayed headline numbers from # their producers, vendored under ./headline/ (per-file sha256 pins: # headline/HEADLINE_SHAS.txt). # kappa = -pi/(24 sqrt2) exact kernel/grading facts (kernel_condition, # exact sympy on the shipped c=2 symbolic-eps # connection) + a PHYSICAL measurement of the # 1/eps pole vector from the integrand's UV # tail at lam=1/2 (pole_symck via the shipped # chain's symck) -- gate: >=30 d vs closed form. # e6/e7 at (a,lam)=(3,3/10) live cosmoflow oracle (vendored package) at a # measured setting -- gate: >=6 d vs the # sha-pinned deep-oracle cache (33.7/33.9 d # certified); deep values printed as CACHE. # V at X*=(7/5,7/10,1) live quadrature of the six assembly terms # T1..T6 (vendored engine_asm on the s7poly # engine class), --jobs-parallel subprocesses # at measured per-term settings -- gate: >=10 d # vs the sha-pinned two-rung assembly cache. # The caches are CACHE ONLY (never the definition): both are sha256-pinned # below and refused on any mutation (rc != 0). # =========================================================================== HEADLINE_DIR = os.path.join(HERE, 'headline') ORACLE_CACHE_SHA256 = \ 'd1ac4669abbfbd909fbfc91cdb7f5122c917c3082bd6bd93075f98abfb6f2f02' ASM_CACHE_SHA256 = \ '6006919e7d2e423c3ee4ccc7cad2690462ee3bee74f1e864611e5559ed373a1b' # measured per-term settings (dps, Nphi, maxdeg) and reference walls in # seconds -- every number measured on a 96-core shared machine (load average # 450-650 during measurement, so an idle laptop-class machine should run # faster; receipts in CHANGES.md). Measured digits vs the deep caches at # these settings: T1/T4/T5/T6 12-13 d, T2/T3 12.5 d, e6 10.8 d / e7 11.6 d. # T2/T3 need the higher rung (their middle quadrature converges slowest -- # the deep run saw the same). HEADLINE_V_SETTINGS = {'T1': (12, 16, 3), 'T2': (14, 24, 4), 'T3': (14, 24, 4), 'T4': (12, 16, 3), 'T5': (12, 16, 3), 'T6': (12, 16, 3)} HEADLINE_ORACLE_SETTING = (8, 8, 3) # (dps, Nch, mdeg) for e6/e7 HEADLINE_WALLS = {'kappa': 40, 'e6e7': 502, 'V_terms': 1665, 'V_slowest': 371} def _hl_load_cache(name, pinned_sha): import hashlib p = os.path.join(HEADLINE_DIR, name) if not os.path.exists(p): raise CertError(f'headline cache {name} MISSING under files/cosmo/headline/ ' '-- it ships with the page; restore it (pins: ' 'headline/HEADLINE_SHAS.txt).') raw = open(p, 'rb').read() sha = hashlib.sha256(raw).hexdigest() if sha != pinned_sha: raise CertError(f'headline cache {name} sha256 MISMATCH: {sha[:20]}... != ' f'pinned {pinned_sha[:20]}... -- the cache is CACHE ONLY ' '(the definition is the vendored producer); a mutated/' 'edited cache is refused, never silently consumed.') return json.loads(raw) def headline(jobs): import subprocess, tempfile t0h = time.time() print('== --headline (Point 5): recompute the displayed headline numbers ' 'from their vendored producers ==', flush=True) total = (HEADLINE_WALLS['kappa'] + HEADLINE_WALLS['e6e7'] + HEADLINE_WALLS['V_slowest']) print(f' measured reference walls (96-core shared machine under heavy load): ' f'kappa ~{HEADLINE_WALLS["kappa"]} s + e6/e7 ~{HEADLINE_WALLS["e6e7"]} s ' f'+ V six terms ~{HEADLINE_WALLS["V_slowest"]} s slowest term on ' f'{jobs} workers; V overlaps the other legs -> projected total ' f'~{max(HEADLINE_WALLS["kappa"] + HEADLINE_WALLS["e6e7"], HEADLINE_WALLS["V_slowest"])}' f'-{total} s', flush=True) ocache = _hl_load_cache('oracle-deep-cache.json', ORACLE_CACHE_SHA256) acache = _hl_load_cache('asm-sum-cache.json', ASM_CACHE_SHA256) print(' caches sha-verified (oracle-deep-cache.json, asm-sum-cache.json)', flush=True) py = sys.executable # ---- V leg: start the six term subprocesses FIRST (they overlap the # kappa and e6/e7 legs), collect at the end. print(f'-- V leg started: T1..T6 engine_asm subprocesses, {jobs} workers, ' f'per-term settings {HEADLINE_V_SETTINGS}', flush=True) vdir = tempfile.mkdtemp(prefix='cosmo-headline-') from concurrent.futures import ThreadPoolExecutor def run_term(T): dps_, nphi_, mdeg_ = HEADLINE_V_SETTINGS[T] out = os.path.join(vdir, f'{T}.json') r = subprocess.run([py, os.path.join(HEADLINE_DIR, 'engine_asm.py'), T, str(dps_), str(nphi_), str(mdeg_), out], capture_output=True, text=True) if r.returncode != 0: raise CertError(f'headline V term {T} FAILED rc={r.returncode}:' f'\n{r.stderr[-800:]}') return json.load(open(out)) pool = ThreadPoolExecutor(max_workers=jobs) futs = {T: pool.submit(run_term, T) for T in ('T1', 'T2', 'T3', 'T4', 'T5', 'T6')} # ---- kappa leg ------------------------------------------------------- t0 = time.time() print('-- kappa leg: exact kernel/grading facts + measured UV-tail pole ' 'vector (lam=1/2) --', flush=True) r = subprocess.run([py, os.path.join(HEADLINE_DIR, 'kernel_condition.py')], capture_output=True, text=True) if r.returncode != 0: raise CertError(f'kernel_condition FAILED rc={r.returncode}:\n{r.stderr[-800:]}') for ln in r.stdout.strip().splitlines(): print(' ', ln) with tempfile.NamedTemporaryFile(suffix='.json', delete=False) as tf: kjson = tf.name r = subprocess.run([py, os.path.join(HEADLINE_DIR, 'pole_symck.py'), '--json', kjson], capture_output=True, text=True) if r.returncode != 0: raise CertError(f'pole_symck FAILED rc={r.returncode}:\n{r.stderr[-800:]}') for ln in r.stdout.strip().splitlines(): print(' ', ln) kd = json.load(open(kjson)) os.unlink(kjson) if not kd['agree_closed_d'] >= 30: raise CertError(f'headline kappa gate FAILED: measured pole/closed-form ' f'agreement {kd["agree_closed_d"]} d < bar 30') print(f' kappa = -pi/(24 sqrt2): PASS ({kd["agree_closed_d"]} d measured ' f'vs closed form) [leg wall {time.time()-t0:.0f} s]', flush=True) # ---- e6/e7 leg ------------------------------------------------------- t0 = time.time() dps_o, nch_o, mdeg_o = HEADLINE_ORACLE_SETTING print(f'-- e6/e7 leg: live cosmoflow oracle at (a,lam)=(3,3/10), c=1, ' f'measured setting dps{dps_o}/Nch{nch_o}/mdeg{mdeg_o} --', flush=True) sys.path.insert(0, HEADLINE_DIR) from cosmoflow.oracle import eval_all_triangle rr, _w = eval_all_triangle(a=3, lam=mp.mpf(3) / 10, eps=0, dps=dps_o, Nch=nch_o, mdeg_mid=mdeg_o, mdeg_out=mdeg_o, keys=['e6', 'e7'], cs=(1,)) with mp.workdps(45): for key in ('e6', 'e7'): v, _ = rr[(key, 1)] v = v.real if hasattr(v, 'real') else v ent = ocache[f'{key}_c1'] dref = mp.mpf(ent['val']) dd = digits(v, dref) print(f' {key} live = {mp.nstr(v, 12)} ({dd:.1f} d vs deep cache)') print(f' {key} deep = {ent["val"]} [CACHE, certified ' f'{ent["d_certified"]} d; producer: deep-oracle certified run, ' f'walls {ocache["configs_wall_s"]} s]') if not dd >= 6: raise CertError(f'headline e6/e7 gate FAILED: {key} live vs deep ' f'{dd:.1f} d < bar 6') print(f' e6/e7: PASS (>=6 d live vs 33.7/33.9 d certified cache) ' f'[leg wall {time.time()-t0:.0f} s]', flush=True) # ---- V leg collect --------------------------------------------------- t0 = time.time() with mp.workdps(40): V = mp.mpf(0) for T in ('T1', 'T2', 'T3', 'T4', 'T5', 'T6'): d = futs[T].result() tv = mp.mpf(d['value']) V += tv cent = acache[T] print(f' {T} live = {d["value"]} ({digits(tv, mp.mpf(cent[0])):.1f} d ' f'vs cache; wall {d["wall"]} s)', flush=True) pool.shutdown() vref = mp.mpf(acache['V_final']) dV = digits(V, vref) print(f' V live = {mp.nstr(V, 16)} ({dV:.1f} d vs cache)') print(f' V deep = {acache["V_final"]} [CACHE, internal certificate ' f'{acache["internal_cert_d"]} d; producer: two-rung assembly ' f'quadrature receipt]') if not dV >= 10: raise CertError(f'headline V gate FAILED: live vs cache {dV:.1f} d < bar 10') print(f' V at X*=(7/5,7/10,1): PASS (>={10} d live vs cache) ' f'[V wall {time.time()-t0:.0f} s from collect; terms overlapped ' f'the other legs]', flush=True) print() print(f'HEADLINE: PASS (kappa >=30 d, e6/e7 >=6 d, V >=10 d; total wall ' f'{time.time()-t0h:.0f} s measured)', flush=True) return 0 def main(): ap = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter, epilog=RECOMPUTE_HELP) ap.add_argument('--dps', type=int, default=40) ap.add_argument('--point', type=str, default='lam=3/5,eps=1/7') ap.add_argument('--orders', type=int, default=2, choices=(0, 1, 2, 3)) ap.add_argument('--check', action='store_true') ap.add_argument('--skip-holonomy', action='store_true') ap.add_argument('--headline', action='store_true', help='Point 5: recompute the displayed headline numbers ' '(kappa, e6/e7 at (3,3/10), V at X*) from their vendored ' 'producers under ./headline/ (measured walls printed; ' 'see the header comment above def headline)') ap.add_argument('--jobs', type=int, default=min(6, os.cpu_count() or 4), help='parallel workers for the --headline V-term subprocesses') ap.add_argument('--boundary-recompute', metavar='COMPS', default=None, help='comma list of boundary components (e6,e7,e8,e9,z1..z5) to ' 'REGENERATE from the chain (definition path) and gate vs the ' 'cache; measured costs in the epilog below') ap.add_argument('--recompute-dps', type=int, default=45, help='target dps for --boundary-recompute (tier: <=41 probe, ' '42-47 control, 48-65 mid, >65 deep legA knobs)') ap.add_argument('--recompute-jobs', type=int, default=min(16, os.cpu_count() or 4), help='parallel shard workers (exact decomposition: any value is ' 'bit-equivalent to serial)') ap.add_argument('--recompute-bar', type=float, default=None, help='override the tier PASS bar (digits vs cache)') ap.add_argument('--chain-dir', default=None, help='regeneration chain location (default: $COSMO_BOUNDARY_CHAIN, ' 'then the shipped ./chain next to this script)') ap.add_argument('--recompute-workdir', default=None, help='shard/combine output dir (resumable; default ./cosmo-' 'boundary-recompute-)') args = ap.parse_args() if args.headline: return headline(args.jobs) if args.boundary_recompute: cb = load_boundary_cache() return boundary_recompute(args.boundary_recompute, args.recompute_dps, args.chain_dir, args.recompute_jobs, args.recompute_bar, args.recompute_workdir, cb) pt = dict(kv.split('=') for kv in args.point.split(',')) lam = Fr(pt.get('lam', '3/5')) epsv = Fr(pt.get('eps', '1/7')) if not (Fr(0) < lam < Fr(1, 3) or Fr(1, 2) <= lam < Fr(7, 10)): print(f'lambda={lam} outside the certified domains (0,1/3) u [1/2,7/10); refusing ' '(honest scope: branches are pinned only there).') sys.exit(1) B = load_bundle() CB = load_boundary_cache() # sha-pinned deep boundary CACHE (Point 4); # fail-closed on any mutation -- rc != 0. print(f'FRW-triangle eps-form evaluator (dps {args.dps}; point lam={lam}, eps={epsv})') print(f'data: cosmo-data.json.gz [A16 artifact sha256 {B["meta"]["A16_sha256"][:16]}...]') print() ok1 = point1(lam, args.dps) print() ok2 = point2(B['A16'], epsv, args.dps) print() ok3 = True if args.skip_holonomy: print('== Point 3: skipped (--skip-holonomy) ==') else: try: import flint # noqa: F401 have_flint = True except ImportError: print('== Point 3: SKIPPED -- python-flint not installed (pip install python-flint) ==') have_flint = False if have_flint: if args.orders >= 3: print('== Point 3 note: --orders 3 evaluates all 47077 eps^3 words ' '(measured: 259 s total wall at --dps 40 on the reference box) ==') print(f'== Point 3: boundary-free per-eps-order holonomy, path 1/2 -> 7/10, ' f'orders eps^0..eps^{args.orders} ==') print(' (prediction = stored exact word lists over closed periods; oracle = live ' 'arb ODE of the stored A16 SUB9 block + exact gauge rotation)') ok3, res, ex = point3(B, args.dps, args.orders) # explicit key into the sha-pinned cache block (order- # independent; the key name is data inside the pinned file) cg = B.get('campaign_gate', {}).get('dps80_lane', {}) if cg: camp = '/'.join(f"{cg[f'eps^{o}']['min_agreed_digits']:.1f}" for o in range(min(args.orders, 3) + 1) if f'eps^{o}' in cg) print(f' (campaign gate, stored provenance only: {camp} d at its dps80 leg)') if not ok3: # promoted to a raise (charter v2) raise CertError('point3 holonomy gate FAILED (bar 30 / zero-mismatch; ' 'see the table above)') if args.check: print() print('== --check: dps-rerun D/D+40, mutation control, n_k cross-check ==') ok3b, res2, ex2 = point3(B, args.dps + 40, args.orders, quiet=True) grew = all(res2[o][0] > res[o][0] + 25 for o in range(args.orders + 1)) for o in range(args.orders + 1): print(f' eps^{o}: {res[o][0]:.1f} d @dps{args.dps} -> {res2[o][0]:.1f} d @dps{args.dps+40}') print(f' dps-scaling: {"PASS (digits grow)" if grew else "FAIL"}') if not grew: raise CertError('--check: D/D+40 rerun digits did not grow by the ' 'calibrated +25 floor (healthy growth ~ +40)') # value-level rerun: the same constants at two genuinely different # internal depths (charter v2 point 3) must agree beyond D digits. with mp.workdps(args.dps + 60): cw, carg = float('inf'), None for kk, pv in ex['pred'].items(): pv2 = ex2['pred'].get(kk) if pv2 is None: continue d = digits(pv, pv2) if d < cw: cw, carg = d, kk print(f' rerun value agreement (all {len(ex["pred"])} gated entries, ' f'D vs D+40): min {cw:.1f} d [{"PASS" if cw > args.dps else "FAIL"} ' f'bar {args.dps}]') if not cw > args.dps: raise CertError(f'--check rerun: entry {carg} agrees only {cw:.1f} d ' f'< dps bar {args.dps}') okm, resm, _exm = point3(B, args.dps, 0, mutate=True, quiet=True) print(f' mutation control (one eps^0 word x(1+1e-30)): eps^0 gate ' f'{resm[0][0]:.1f} d -- {"PASS (collapses to ~30)" if resm[0][0] < 33 else "FAIL"}') if not resm[0][0] < 33: raise CertError('--check mutation control FAILED to collapse: the gate ' 'is not consuming the live word values') _nk_quad_check(B, args.dps, ex['nk']) ok3 = ok3 and ok3b and grew and resm[0][0] < 33 print() ok4 = point4(CB) print() verdict = ok1 and ok2 and ok3 and ok4 print(f'OVERALL: {"PASS" if verdict else "FAIL"} (wall {time.time()-T0:.0f}s)') return 0 if verdict else 1 if __name__ == '__main__': sys.exit(main())