#!/usr/bin/env python3 """hexabox_vop.py — the sec255 block-analytic VoP layer-recursion engine, packaged with the bundle hexabox-vop-data.json.gz (no external imports). What it does at runtime, from the exact rational graded path DE: * derives every block's homogeneous fundamental EXACTLY (rational solution vectors / sqrt-letter radical columns / triangular+conjugate reductions), certifying each strategy symbolically (zero-residual identities over fmpq); * Tier 1: builds the fully-expanded symbolic closed forms n[i,k](t) = sum C*sqrt(rad)*R(t)*G_word(t) for layers -4..0 and re-proves the zero-remainder exact-DE certificate for the top rows; * Tier 2: evaluates the same VoP construction on a Chebyshev spectral contour (layers to eps^0 of the top sector), i.e. the layered iterated integrals over closed-form kernels — precision grows with NC and mp.mp.dps. All kernels are exact fmpq rationals (bundle["exactA"]); bundle["hybrid_kernels"] is empty since the (97,0,2..4) succession and the numeric-kernel path below is kept only to read superseded bundles, whose output it labels. Public deps: mpmath, sympy, python-flint. CHANGELOG 2026-07-05: (i) boundary-seed PSLQ naming (try_rational / try_logring) pinned at NAMING_DPS=80 inside mp.workdps — ambient dps now tracks the caller's --dps, naming stays byte-identical across caller precisions (seed strings are 65-70 d; naming is re-certified downstream by the exact-DE certificates + the endpoint checks); (ii) setup_engine resets NUM/hybrid_used so the refine-until-bound loop can re-run the SAME recursion at a refined NC; (iii) top_table_bound: a posteriori trailing-window Chebyshev tail bounds of the top-sector solution tables (certified-bound machinery). REACH CHECK 2026-09-07 (Vop255.reach_check; every value path is untouched where it does not fire). The recursion V_i^(k) = sum_j sum_{m in SUPP(i,j)} A_ij^(m) V_j^(k-m) reads only the layers m the bundle carries: a layer m ABSENT from SUPP(i,j) that lies within the reach of the requested layer k -- k - k_min(j) >= m, where k_min(j) is the lowest live layer of the source row j (its first nonzero fit-point seed order plus L_j; -1 for a sunrise row, whose closed containers start there) -- used to be dropped silently. Now a coupling (i,j) whose max m present is below k - k_min(j) is SHORT and the recursion refuses by name (RuntimeError 'VoP reach: ...') BEFORE the sum: (a) for every coupling at every layer above KMAX_RECORD = 3 (the layer through which this bundle's supports and seeds were computed and checked against the independent oracle, eps^0 of the top sector); (b) for the top-sector couplings at every layer (on this bundle each of the 294 ends exactly at the layer-3 reach, so a missing top layer is a drop); (c) for a row computed above KMAX_RECORD with no fit-point seed at that layer (the seeds stop at eps^2 = layer L_i + 2; the sunrise closed containers at layer 5). NOT a short coupling: an absent m BELOW max m present (a zero coefficient of the exact expansion; no coupling of this bundle has such a gap), and a non-top coupling whose support ends below the layer-3 reach (350 of them; the default run's agreement with the independent oracle establishes those absent layers as zero coefficients -- the bundle carries no per-entry expansion order, so such a complete expansion cannot be told from a dropped layer below KMAX_RECORD). Once per run the per-row k_min table and the reach summary are printed ('[vop/reach]' lines); with --kmax 5 on the bundle served before 2026-09-07 every top-sector coupling was short at layer 4 and the run refused before any layer was computed. SUCCESSION 2026-09-07: KMAX_RECORD is a per-bundle declaration -- the class default 3 is the served bundle's; a driver that has gated a bundle deeper declares that layer on the instance before the check (hexabox-evaluate.py declares 5 for the succeeded bundle hexabox-vop-data.json.gz by sha, the layer through which the gate of record checked it against both oracle records), so that above the record layer the same three rules apply and the non-top couplings whose supports end below the reach at layers 4, 5 (the (61,61) class, complete expansions by the oracle reading) are not refused; --kmax 6 on that bundle refuses by name at layer 6. """ import gzip, json, os, sys, time from fractions import Fraction import mpmath as mp import flint import sympy as sp NAMING_DPS = 80 # naming runs at 80 digits, the precision at which the seeds were named # flipped: lives in the blog dir natively import hexabox_vop_lib as VL from hexabox_vop_lib import (FQ, FP, fq, fp, Rat, R_ZERO, R_ONE, ZERO_P, ONE_P, LETTERS, reg_letter, reg_const, CONSTS, rad_dlog_half, cf_zero, cf_const, cf_add, cf_scale, cf_mulrat, cf_mulalg, cf_diff, cf_int, _neg_cf, hermite, Engine, solve_rat_solutions, _mpc2acb, _acb2mpc) VARS = ['s12', 's23', 's34', 's45', 's15', 'mm'] def _neg(A): return _neg_cf(A) class Vop255: """One full live closure of the sec255 graded path DE from the bundle.""" def __init__(self, bundle_path, log=print): self.log = log t0 = time.time() raw = gzip.open(bundle_path, 'rb').read() B = json.loads(raw) self.B = B spec = B["D_spec"] self.N = spec["N"]; self.TOP = spec["TOP"]; self.SUN = spec["SUN"] self.SUNset = set(self.SUN); self.L = spec["L"] self.SECTORS = spec["sectors"] self.MASTERS = [tuple(m) for m in spec["masters"]] # exact rational DE entries (i,j,m) -> Rat (ALL exact fmpq) self.RatA = {} for key, (Ps, Qs) in B["exactA"].items(): i, j, m = (int(x) for x in key.split(',')) self.RatA[(i, j, m)] = Rat(fp([Fraction(c) for c in Ps]), fp([Fraction(c) for c in Qs])) self.COLS = {}; self.SUPP = {} for (i, j, m) in self.RatA: self.COLS.setdefault(i, set()).add(j) self.SUPP.setdefault((i, j), []).append(m) # the 3 disclosed certified-numeric kernels (row 97 <- row 0, m=2,3,4) self.HYBRID = {} for key, (Ps, Qs) in B["hybrid_kernels"].items(): i, j, m = (int(x) for x in key.split(',')) self.HYBRID[(i, j, m)] = ([mp.mpf(c) for c in Ps], [mp.mpf(c) for c in Qs]) self.COLS.setdefault(i, set()).add(j) self.hybrid_used = [] # fit-point boundary seeds: derived values (see bundle _fold); the # retired interim strings are kept as a compare-on-load cross-check self.V0 = {} for si, row in B["boundary_Pfit"].items(): i = int(si) for so, (re_, im_) in row.items(): self.V0[(i, int(so) + self.L[i])] = mp.mpc(mp.mpf(re_), mp.mpf(im_)) if "_fold" in B: cache = B.get("boundary_Pfit_interim_cache", {}) worst = mp.mpf('inf'); nchk = 0 with mp.workdps(100): for si, row in cache.items(): for so, (re_, im_) in row.items(): new = B["boundary_Pfit"][si][so] nv = mp.mpc(mp.mpf(new[0]), mp.mpf(new[1])) ov = mp.mpc(mp.mpf(re_), mp.mpf(im_)) sc = max(abs(nv), abs(ov)); nchk += 1 if sc == 0: continue d = abs(nv - ov) ag = mp.mpf(100) if d == 0 else -mp.log10(d/sc) worst = min(worst, ag) if ag < 60: raise RuntimeError( "boundary cache mismatch: boundary_Pfit row %s " "order %s, derived value vs retired interim " "string agree to only %s digits (below 60; at " "least 69 expected) — the data bundle has been " "altered; refusing to run" % (si, so, mp.nstr(ag, 4))) self.log("[vop] boundary seeds: %d derived values cross-checked on " "load against the retired interim strings, worst agreement " "%s digits (threshold 60)." % (nchk, mp.nstr(worst, 4))) # boundary-constant log ring self.LOG_TAGS = ['1', 'log2', 'log3', 'log5', 'log7', 'log11', 'log13', 'log17', 'ipi'] self.LOG_BASIS = [mp.mpf(1), mp.log(2), mp.log(3), mp.log(5), mp.log(7), mp.log(11), mp.log(13), mp.log(17), mp.pi*1j] for t_, v_ in zip(self.LOG_TAGS, self.LOG_BASIS): reg_const(t_, mp.mpc(v_)) # blocks + topological order self.BLOCKS = {} for i in range(self.N): if i in self.SUNset: continue self.BLOCKS.setdefault(self.SECTORS[i], []).append(i) self.BLOCKS = {s: sorted(r) for s, r in self.BLOCKS.items()} row_block = {i: s for s, rows in self.BLOCKS.items() for i in rows} edges = {s: set() for s in self.BLOCKS} for (i, j) in self.SUPP: if i in self.SUNset or j in self.SUNset: continue si, sj = row_block[i], row_block[j] if si != sj: edges[si].add(sj) self.ORDER = []; seen = set() def visit(s): if s in seen: return seen.add(s) for p in sorted(edges[s]): visit(p) self.ORDER.append(s) for s in sorted(self.BLOCKS): visit(s) # same-layer (m=0) ancestor set of TOP: rows whose top layer matters self.C21 = set(); stack = list(self.TOP) while stack: i = stack.pop() if i in self.C21: continue self.C21.add(i) for j in self.COLS.get(i, []): if (i, j, 0) in self.RatA and j not in self.C21 and j not in self.SUNset: stack.append(j) self._build_sunrise() self.BAS_CACHE = {} self.NCF = {} # symbolic containers (i,k) self.NUM = {} # numeric node tables (i,k) self.boundary_record = {} self.log(f"[vop] bundle loaded: {len(self.RatA)} exact entries + " f"{len(self.HYBRID)} certified-numeric, {len(self.BLOCKS)} blocks " f"({time.time()-t0:.1f}s)") # ── analytic sunrise closed containers (sources of the recursion) ── def _build_sunrise(self): SY = {v: sp.Symbol(v) for v in VARS} tt = sp.Symbol('t') Pa, Pb = self.B["path"]["Pfit"], self.B["path"]["gate"] _epsC = lambda e: -mp.gamma(1-e)**3*mp.gamma(1+2*e)/(2*(2*e-1)*mp.gamma(3-3*e)) _cT = mp.taylor(_epsC, 0, 16) ccoef = {k: _cT[k+1] for k in range(-1, 16)} self.SUN_CF = {} for j in self.SUN: e = sp.sympify(self.B["sun_p2"][str(self.SECTORS[j])], locals=SY) ex = e.subs({SY[v]: sp.Rational(Pa[v]) + tt*(Pb[v]-Pa[v]) for v in VARS}) cs = [Fraction(int(sp.Rational(c).p), int(sp.Rational(c).q)) for c in reversed(sp.Poly(ex, tt).all_coeffs())] P2R = Rat(fp(cs)) P2_0 = P2R.ev_fq(0) P2_0mp = mp.mpf(int(P2_0.p))/int(P2_0.q) L0 = mp.log(-P2_0mp) if P2_0mp < 0 else mp.log(P2_0mp) - mp.pi*1j li = reg_letter(P2R.n) if P2R.n.degree() >= 1 else None for k in range(-1, 6): C = cf_zero() for q in range(0, k+2): val = mp.mpc(0) for r in range(-1, k+1): p = k - r if q > p: continue val += ccoef[r] * (-2)**p * L0**(p-q) / mp.factorial(p-q) if abs(val) < mp.mpf(10)**-70: continue if li is None and q > 0: continue tag = reg_const(f"SUN{j}k{k}q{q}", mp.mpc(val)) word = (('L', li),)*q if li is not None else () C = cf_add(C, {(tag, (), word): -P2R}) self.SUN_CF[(j, k)] = C maxe = mp.mpf(0) for j in self.SUN: for k in range(-1, 3): v = sum(CONSTS[t2]*r.ev_mpc(mp.mpc(0)) for (t2, rd, w), r in self.SUN_CF[(j, k)].items() if w == ()) maxe = max(maxe, abs(v - self.V0.get((j, k), mp.mpc(0)))) assert maxe < mp.mpf(10)**-60, f"sunrise boundary check FAIL {maxe}" self.sunrise_t0_err = maxe # ── boundary-seed PSLQ naming (fit-side; identical to the original derivation) ── # 2026-07-05: naming runs pinned at NAMING_DPS = 80, the precision at which the # seeds were named, inside mp.workdps, so the CLASSIFICATION of the 65-70 digit # seed strings is byte-identical no matter what ambient dps the caller set # (ambient now tracks --dps). Naming precision is bounded by the seed STRING # length, not by --dps; every naming decision is re-certified downstream by the # zero-remainder exact-DE certificates and the independent endpoint checks. def try_rational(self, x, maxc=10**15, tol=None): with mp.workdps(NAMING_DPS): tol = tol or mp.mpf(10)**-52 if abs(mp.im(x)) > tol: return None if mp.re(x) == 0: return Fraction(0) # PSLQ rejects zero vectors r = mp.pslq([mp.re(x), mp.mpf(1)], maxcoeff=maxc, maxsteps=10**6) if not r or r[0] == 0: return None v = Fraction(-r[1], r[0]) if abs(mp.re(x) - mp.mpf(v.numerator)/v.denominator) < tol: return v return None def try_logring(self, x, tol=None): with mp.workdps(NAMING_DPS): tol = tol or mp.mpf(10)**-48 vec = [mp.re(x)] + [mp.re(b) for b in self.LOG_BASIS[1:8]] + [mp.mpf(1)] r = mp.pslq(vec, maxcoeff=10**4, maxsteps=10**6) out = None if r and r[0] != 0: den = -r[0] co = {self.LOG_TAGS[a+1]: Fraction(r[a+1], den) for a in range(7)} co['1'] = Fraction(r[8], den) chk = sum(mp.mpf(f.numerator)/f.denominator * self.LOG_BASIS[self.LOG_TAGS.index(t2)] for t2, f in co.items()) if abs(mp.re(x)-chk) < tol: out = co if out is None: return None qpi = self.try_rational(mp.im(x)/mp.pi, maxc=10**6) if qpi is None and abs(mp.im(x)) > tol: return None if qpi: out['ipi'] = qpi return {t2: f for t2, f in out.items() if f != 0} # ── exact block-strategy machinery (verbatim port) ── def eig_generic(self, Amat): def ef(f): rts = mp.polyroots([mp.mpf(int(c.p))/int(c.q) for c in reversed(f.coeffs())], maxsteps=300, extraprec=80) r0 = rts[0]; nb = len(Amat) Mm = mp.matrix(nb, nb); K = 8; eps = mp.mpf('1e-18') for a in range(nb): for b in range(nb): if Amat[a][b].is_zero(): continue s = mp.mpc(0) for kk in range(K): z = r0 + eps*mp.exp(2j*mp.pi*kk/K) s += (z-r0)*Amat[a][b].ev_mpc(z) Mm[a, b] = s/K try: return mp.eig(Mm, left=False, right=False) except Exception: return [] return ef @staticmethod def certify_generic(Amat, vec, rad): nb = len(Amat) half = rad_dlog_half(rad) if rad else R_ZERO for a in range(nb): resid = vec[a].deriv() + vec[a]*half for b in range(nb): if Amat[a][b].is_zero(): continue resid = resid - Amat[a][b]*vec[b] if not resid.is_zero(): return False return True def scalar_strategy(self, b): if b.is_zero(): return ('FULL', [[R_ONE]], [()]) H, logs, rems = hermite(b) if not H.is_zero() or rems: raise RuntimeError("scalar mu non-elementary") num = ONE_P; den = ONE_P; rad = [] for li, cl in logs.items(): c2 = Fraction(int(cl.p), int(cl.q)) if c2.denominator == 1: if c2 >= 0: num = num * LETTERS[li]**c2.numerator else: den = den * LETTERS[li]**(-c2.numerator) elif c2.denominator == 2: fl = (c2.numerator - 1)//2 if fl >= 0: num = num * LETTERS[li]**fl else: den = den * LETTERS[li]**(-fl) rad.append(li) else: raise RuntimeError(f"scalar non-half-integer exponent {c2}") col = Rat(num, den); rad = tuple(sorted(rad)) assert self.certify_generic([[b]], [col], rad), "scalar cert fail" return ('FULL', [[col]], [rad]) def find_solutions(self, Amat): nb = len(Amat); facs = {} for a in range(nb): for b in range(nb): if Amat[a][b].is_zero(): continue c, fl = Amat[a][b].d.factor() for f, e in fl: facs[tuple(str(x) for x in f.coeffs())] = f ef = self.eig_generic(Amat) halfset = [] for k2, f in facs.items(): for x in ef(f): fr = float(mp.re(x))*2 if abs(fr - round(fr)) < 1e-5 and round(fr) % 2 != 0: halfset.append(f); break cols = []; rads = [] for v in solve_rat_solutions(Amat, None, eig_fn=ef): if len(cols) < nb and self.certify_generic(Amat, v, ()): cols.append(v); rads.append(()) if len(cols) < nb and halfset: radli = tuple(sorted(reg_letter(f) for f in halfset)) for v in solve_rat_solutions(Amat, radli, eig_fn=ef): if len(cols) < nb and self.certify_generic(Amat, v, radli): cols.append(v); rads.append(radli) return cols, rads @staticmethod def matmul_rat(X, Y): n1 = len(X); n2 = len(Y[0]) out = [[R_ZERO for _ in range(n2)] for _ in range(n1)] for a in range(n1): for b in range(n2): acc = R_ZERO for c in range(len(Y)): if X[a][c].is_zero() or Y[c][b].is_zero(): continue acc = acc + X[a][c]*Y[c][b] out[a][b] = acc return out @staticmethod def rat_matinv(Rm): nb = len(Rm) if nb == 1: return [[Rm[0][0].inv()]], Rm[0][0] if nb == 2: det = Rm[0][0]*Rm[1][1] - Rm[0][1]*Rm[1][0] di = det.inv() return [[Rm[1][1]*di, -(Rm[0][1]*di)], [-(Rm[1][0]*di), Rm[0][0]*di]], det det = (Rm[0][0]*(Rm[1][1]*Rm[2][2]-Rm[1][2]*Rm[2][1]) - Rm[0][1]*(Rm[1][0]*Rm[2][2]-Rm[1][2]*Rm[2][0]) + Rm[0][2]*(Rm[1][0]*Rm[2][1]-Rm[1][1]*Rm[2][0])) di = det.inv() def cof(r2, c2): rs = [x for x in range(3) if x != r2]; cs2 = [x for x in range(3) if x != c2] v = Rm[rs[0]][cs2[0]]*Rm[rs[1]][cs2[1]] - Rm[rs[0]][cs2[1]]*Rm[rs[1]][cs2[0]] return v if (r2+c2) % 2 == 0 else -v return [[cof(c2, r2)*di for c2 in range(3)] for r2 in range(3)], det def strategy_for(self, Amat, tag=''): nb = len(Amat) if all(Amat[a][b].is_zero() for a in range(nb) for b in range(nb)): return ('FULL', [[R_ONE if a == c else R_ZERO for a in range(nb)] for c in range(nb)], [()]*nb) if nb == 1: return self.scalar_strategy(Amat[0][0]) cols, rads = self.find_solutions(Amat) if len(cols) == nb: return ('FULL', cols, rads) ratcols = [cols[c] for c in range(len(cols)) if rads[c] == ()] if not ratcols: radlis = [rads[c] for c in range(len(cols)) if rads[c] != ()] if radlis: radli = radlis[0] g = ONE_P for li in radli: g = g * LETTERS[li] half = Rat(g.derivative(), g*FP([2])) Ag = [[(Amat[a][b] - half if a == b else Amat[a][b]) for b in range(nb)] for a in range(nb)] sub = self.strategy_for(Ag, tag+"/conj") return ('CONJ', radli, sub) raise RuntimeError(f"{tag}: no rational homogeneous column") f = len(ratcols) from itertools import combinations T = Tinv = None for unit_idx in combinations(range(nb), nb-f): Tc = [[(ratcols[c][a] if c < f else (R_ONE if a == unit_idx[c-f] else R_ZERO)) for c in range(nb)] for a in range(nb)] try: Ti, det = self.rat_matinv(Tc) if det.is_zero(): continue if det.ev_fq(0) == 0: continue T, Tinv = Tc, Ti break except ZeroDivisionError: continue if T is None: raise RuntimeError(f"{tag}: no invertible T padding found") Tp = [[T[a][c].deriv() for c in range(nb)] for a in range(nb)] AT = self.matmul_rat(Amat, T) ATmTp = [[AT[a][c] - Tp[a][c] for c in range(nb)] for a in range(nb)] Bm = self.matmul_rat(Tinv, ATmTp) for c in range(f): for a in range(nb): assert Bm[a][c].is_zero(), f"{tag}: B col {c} not zero" corner = [[Bm[f+a][f+b] for b in range(nb-f)] for a in range(nb-f)] sub = self.strategy_for(corner, tag+f"/corner{nb-f}") return ('TRI', T, Tinv, Bm, f, sub) def block_strategy(self, s): if s in self.BAS_CACHE: return self.BAS_CACHE[s] rows = self.BLOCKS[s]; nb = len(rows) Amat = [[self.RatA.get((rows[a], rows[b], 0), R_ZERO) for b in range(nb)] for a in range(nb)] st = self.strategy_for(Amat, tag=f"sec{s}{rows}") self.BAS_CACHE[s] = st return st # ── Tier-1 symbolic closure (containers) ── def solve_lin(self, strat, Svec, w0cf): if strat[0] == 'CONJ': _, radli, sub = strat invfac = R_ONE for li in radli: fpoly = LETTERS[li] invfac = invfac * Rat(FP([fpoly(fq(0))]), fpoly) Smod = [cf_mulalg(S, invfac, radli) if S else cf_zero() for S in Svec] Z = self.solve_lin(sub, Smod, w0cf) return [cf_mulalg(z, R_ONE, radli) if z else cf_zero() for z in Z] if strat[0] == 'FULL': _, cols, rads = strat nb = len(cols) Rm = [[cols[c][a] for c in range(nb)] for a in range(nb)] Rinv, det = self.rat_matinv(Rm) Ys = [cf_zero() for _ in range(nb)] for c in range(nb): Tc = cf_zero() for b in range(nb): if Rinv[c][b].is_zero() or not Svec[b]: continue Tc = cf_add(Tc, cf_mulrat(Svec[b], Rinv[c][b])) if rads[c] and Tc: invfac = R_ONE for li in rads[c]: fpoly = LETTERS[li] invfac = invfac * Rat(FP([fpoly(fq(0))]), fpoly) Tc = cf_mulalg(Tc, invfac, rads[c]) Ic = cf_int(Tc) if Tc else cf_zero() C0 = cf_zero() for b in range(nb): if Rinv[c][b].is_zero() or not w0cf[b]: continue rv = Rinv[c][b].ev_fq(0) if rv != 0: C0 = cf_add(C0, cf_scale(w0cf[b], rv)) inner = cf_add(C0, Ic) if not inner: continue for a in range(nb): if Rm[a][c].is_zero(): continue Ys[a] = cf_add(Ys[a], cf_mulalg(inner, Rm[a][c], rads[c])) return Ys _, T, Tinv, Bm, f, sub = strat nb = len(T) TS = []; W0 = [] for c in range(nb): acc = cf_zero() for b in range(nb): if Tinv[c][b].is_zero() or not Svec[b]: continue acc = cf_add(acc, cf_mulrat(Svec[b], Tinv[c][b])) TS.append(acc) acc0 = cf_zero() for b in range(nb): if Tinv[c][b].is_zero() or not w0cf[b]: continue rv = Tinv[c][b].ev_fq(0) if rv != 0: acc0 = cf_add(acc0, cf_scale(w0cf[b], rv)) W0.append(acc0) wsub = self.solve_lin(sub, TS[f:], W0[f:]) w = [None]*nb for c in range(nb-f): w[f+c] = wsub[c] for a in range(f): integ = TS[a] for b in range(f, nb): if Bm[a][b].is_zero() or not w[b]: continue integ = cf_add(integ, cf_mulrat(w[b], Bm[a][b])) w[a] = cf_add(W0[a], cf_int(integ) if integ else cf_zero()) Ys = [cf_zero() for _ in range(nb)] for a in range(nb): for c in range(nb): if T[a][c].is_zero() or not w[c]: continue Ys[a] = cf_add(Ys[a], cf_mulrat(w[c], T[a][c])) return Ys def getn(self, j, k): if j in self.SUNset: return self.SUN_CF.get((j, k), cf_zero()) if k >= -1 else cf_zero() return self.NCF.get((j, k), cf_zero()) # ── the reach check (2026-09-07): an absent in-reach layer refuses by name instead of dropping ── KMAX_RECORD = 3 # the layer through which the bundle's layer supports and fit-point seeds were # computed and checked against the independent oracle (the default run: layers # -4..3 = eps^-4..eps^0 of the top sector); every layer above it is governed # by reach_check for every coupling, the top-sector couplings at every layer def reach_check(self, KMAX, k=None, rows=None, report=False): """The reach check (module docstring). Built once per run: KMIN[j] = the lowest live layer of source row j (its first fit-point seed order with |seed| >= 1e-65, plus L_j; -1 for a sunrise row, whose closed containers start there; None when no seed is live), KSEED[j] = the last seeded layer (seeds stop at eps^2 = layer L_j + 2; sunrise containers at layer 5), MAXM[(i,j)] = the highest layer m the bundle carries for the coupling (exact entries and the disclosed certified-numeric kernels alike). Coupling (i,j) is SHORT at layer k when MAXM < k - KMIN[j]: the sum would read the live source V_j^(k-m) for an m the bundle lacks. Checked for every coupling at k > KMAX_RECORD and for the top-sector couplings at every k; a non-top coupling at k <= KMAX_RECORD is not checked (see the docstring: complete expansions the bundle cannot tell from a drop). A row computed at k > KMAX_RECORD with no fit-point seed at layer k is a SHORT SEED; a sunrise source read above its containers' last layer is SHORT too. k is None: project every layer -4..KMAX over the rows each layer computes (the top layer restricted to the same-layer ancestors of TOP exactly as run_numeric does) and refuse before any layer is computed; k given: check `rows` at that layer (the call before each block's sum). report=True prints the per-row table and the reach summary once per run. Raises RuntimeError('VoP reach: ...') naming every short coupling and seed.""" if not hasattr(self, 'KMIN'): thr = mp.mpf(10)**-65 self.KMIN = {}; self.KSEED = {}; self.MAXM = {} for j in range(self.N): if j in self.SUNset: self.KMIN[j] = -1 self.KSEED[j] = max(kk for (jj, kk) in self.SUN_CF if jj == j) continue ks = [kk for (jj, kk), x in self.V0.items() if jj == j and abs(x) >= thr] self.KMIN[j] = min(ks) if ks else None self.KSEED[j] = max([kk for (jj, kk) in self.V0 if jj == j] or [-5]) for (i, j), ms in self.SUPP.items(): self.MAXM[(i, j)] = max(ms) for (i, j, m) in self.HYBRID: self.MAXM[(i, j)] = max(self.MAXM.get((i, j), m), m) self._reach_reported = False KR = self.KMAX_RECORD; TOPset = set(self.TOP) def layer_rows(kk): out = [] for s in self.ORDER: rr = self.BLOCKS[s] if KMAX >= 3 and kk == KMAX and not any(i in self.C21 for i in rr): continue out.extend(rr) return out plan = [(k, list(rows))] if k is not None else [(kk, layer_rows(kk)) for kk in range(-4, KMAX+1)] short = []; seeds = []; checked = set() for kk, rr in plan: for i in rr: if kk > KR and i not in self.SUNset and self.KSEED[i] < kk: seeds.append((i, kk, self.KSEED[i])) if kk <= KR and i not in TOPset: continue for j in self.COLS.get(i, []): checked.add((i, j)) M = self.MAXM[(i, j)]; kj = self.KMIN[j] if kj is not None and M < kk - kj: short.append(('coupling', i, j, kk, M, kk - kj, kj)) if j in self.SUNset and (i, j) in self.SUPP and kk - min(self.SUPP[(i, j)]) > self.KSEED[j]: short.append(('sunrise', i, j, kk, M, kk - min(self.SUPP[(i, j)]), self.KSEED[j])) if report and not self._reach_reported: self._reach_reported = True cons = {} for (i, j), ms in self.SUPP.items(): c = cons.setdefault(j, [0, 99, -1]) c[0] += 1; c[1] = min(c[1], min(ms)); c[2] = max(c[2], self.MAXM[(i, j)]) gaps = sum(max(ms) - min(ms) + 1 - len(ms) for ms in self.SUPP.values()) self.log(f"[vop/reach] per-row k_min (the lowest live layer of every source row) and the reach at " f"KMAX {KMAX} (record layer {KR}; {len(self.SUPP)} couplings, {gaps} absent layers below a " f"coupling's max m present):") for j in range(self.N): kj = self.KMIN[j] if j in self.SUNset: src = "sunrise closed containers" elif kj is None: src = "no live seed" else: src = f"seed eps^{kj - self.L[j]:+d}, L {self.L[j]}" c = cons.get(j) self.log(f"[vop/reach] row {j}: k_min {kj} ({src}), seeded to layer {self.KSEED[j]}, " + (f"consumed by {c[0]} couplings (m present {c[1]}..{c[2]})" if c else "consumed by no coupling") + f", max m needed at KMAX {KMAX}: {(KMAX - kj) if kj is not None else 'n/a'}") if short or seeds: by_row = {} for t in short: by_row.setdefault(t[1], []).append(t) for i in sorted(by_row): items = sorted(by_row[i], key=lambda t: (t[2], t[3])) first = {} for t in items: first.setdefault((t[0], t[2]), t) self.log(f"[vop/reach] short of the reach, row {i} ({len(first)} coupling(s)): " + "; ".join( (f"({t[1]},{t[2]}) m present {t[4]}, needed {t[5]} at layer {t[3]} (k_min({t[2]}) = {t[6]})" if t[0] == 'coupling' else f"({t[1]},{t[2]}) sunrise source layer {t[5]} at layer {t[3]} above the containers' last layer {t[6]}") for t in first.values())) if seeds: self.log("[vop/reach] fit-point seed absent: " + "; ".join( f"row {i} at layer {kk} (seeded to layer {ks})" for (i, kk, ks) in seeds)) if report: n_c = len({(t[1], t[2]) for t in short}); n_s = len({i for (i, kk, ks) in seeds}) self.log(f"[vop/reach] reach complete for KMAX {KMAX}: {'yes' if not (short or seeds) else 'NO'} " f"({n_c} short coupling(s), {n_s} row(s) without a seed at a computed layer; " f"{len(checked)} couplings checked)") if short or seeds: lead = sorted(short, key=lambda t: (t[1] not in TOPset, t[6] if t[6] is not None else 99, t[3], t[1], t[2])) lead = [t for n, t in enumerate(lead) if all(t[1:3] != u[1:3] for u in lead[:n])][:8] names = "; ".join( (f"coupling ({t[1]},{t[2]}) layer {t[4] + 1} absent within the reach of layer {t[3]} " f"(k_min({t[2]}) = {t[6]}, max m present {t[4]})" if t[0] == 'coupling' else f"coupling ({t[1]},{t[2]}) reads sunrise layer {t[5]} at layer {t[3]} above the containers' last layer {t[6]}") for t in lead) if seeds: names += ("; " if names else "") + f"row {seeds[0][0]} has no fit-point seed at layer {seeds[0][1]}" raise RuntimeError( f"VoP reach: {len({(t[1], t[2]) for t in short})} coupling(s) short of the reach at KMAX {KMAX} and " f"{len({i for (i, kk, ks) in seeds})} row(s) without a seed at a computed layer -- {names}" + (f" [+{len({(t[1], t[2]) for t in short}) - len(lead)} more coupling(s), every one named in the [vop/reach] lines above]" if len({(t[1], t[2]) for t in short}) > len(lead) else "") + ": the bundle's expansion stops below the requested layer; enter the layers or lower --kmax") return short, seeds def close_symbolic(self, KSYM=0): """Layers -4..KSYM as fully-expanded exact containers (Tier 1).""" self.reach_check(KSYM) for k in range(-4, KSYM+1): tL = time.time() for s in self.ORDER: rows = self.BLOCKS[s] self.reach_check(KSYM, k=k, rows=rows) Svec = [] for i in rows: S = cf_zero() for j in self.COLS.get(i, []): for m in self.SUPP.get((i, j), []): if m == 0 and j in rows: continue src = self.getn(j, k-m) if not src: continue S = cf_add(S, cf_mulrat(src, self.RatA[(i, j, m)])) for (fi, fj, fm) in self.HYBRID: if fi == i and self.getn(fj, k-fm): raise RuntimeError( f"certified-numeric kernel ({fi},{fj},{fm}) needed at " f"layer {k}: Tier-1 symbolic scope must stay below it") Svec.append(S) v0 = [] for i in rows: x = self.V0.get((i, k)) v0.append(None if (x is None or abs(x) < mp.mpf(10)**-65) else x) if all(v is None for v in v0) and all(not S for S in Svec): continue strat = self.block_strategy(s) v0cf = [] for a, i in enumerate(rows): if v0[a] is None: v0cf.append(cf_zero()); continue q = self.try_rational(v0[a]) if q is not None: v0cf.append({('1', (), ()): Rat(fp([q]))}) self.boundary_record[f"{i},{k}"] = f"exact {q}" else: lr = self.try_logring(v0[a]) if lr is not None: C = cf_zero() for t2, f in lr.items(): C = cf_add(C, {(t2, (), ()): Rat(fp([f]))}) v0cf.append(C) self.boundary_record[f"{i},{k}"] = "logring " + str(lr) else: tag = reg_const(f"B{i}k{k}", v0[a]) v0cf.append(cf_const(tag)) self.boundary_record[f"{i},{k}"] = f"tagged {mp.nstr(v0[a],35)}" Ys = self.solve_lin(strat, Svec, v0cf) for a, i in enumerate(rows): if Ys[a]: self.NCF[(i, k)] = Ys[a] nnz = sum(1 for (i2, k2) in self.NCF if k2 == k) self.log(f"[vop/T1] layer {k}: {nnz} expanded rows ({time.time()-tL:.1f}s)") return self.NCF def de_certificate(self, rows, kmax=0): """Zero-remainder exact-DE check of the Tier-1 expanded containers.""" out = {} self.reach_check(kmax) for i in rows: for k in range(-1, kmax+1): Y = self.NCF.get((i, k)) if Y is None: continue D = cf_diff(Y) self.reach_check(kmax, k=k, rows=[i]) for j in self.COLS.get(i, []): for m in self.SUPP.get((i, j), []): src = self.getn(j, k-m) if not src: continue D = cf_add(D, _neg(cf_mulrat(src, self.RatA[(i, j, m)]))) bad = [kk for kk, r in D.items() if not r.is_zero()] out[(i, k)] = (len(Y), not bad) return out # ── Tier-2 numeric route: same construction on the spectral contour ── def setup_engine(self, NC): polys = [] def walk(st): if st[0] == 'FULL': cols = st[1]; nb = len(cols) Rm = [[cols[c][a] for c in range(nb)] for a in range(nb)] _, det = self.rat_matinv(Rm) if det.n.degree() >= 1: polys.append(det.n) elif st[0] == 'CONJ': walk(st[2]) else: _, det = self.rat_matinv(st[1]) if det.n.degree() >= 1: polys.append(det.n) walk(st[5]) for s in self.ORDER: walk(self.block_strategy(s)) roots = set() for p in polys: cs = [mp.mpf(int(c.p))/int(c.q) for c in reversed(p.coeffs())] try: rr = mp.polyroots(cs, maxsteps=300, extraprec=80) except Exception: continue for r in rr: if abs(mp.im(r)) < 1e-9 and 1e-6 < mp.re(r) < 1-1e-6: roots.add(round(float(mp.re(r)), 8)) self.det_roots = sorted(roots) self.eng = Engine(sorted(set(self.B["real_poles"]) | roots), NC=NC) self.NC = NC; self.NSEG = len(self.eng.segs) self._RV = {}; self._RT = {}; self._ST = {}; self._HT = {} # 2026-07-05: reset the numeric state so a refine-until-bound re-run at a # refined NC is the SAME recursion on fresh tables (no stale-NC leftovers) self.NUM = {}; self.hybrid_used = [] return self.eng def top_table_bound(self, window=8): """A posteriori certified spectral tail bound of the Tier-2 TOP-sector solution tables (2026-07-05): per (row, layer) node table and per contour segment, the trailing-`window` Chebyshev coefficient max, geometric-envelope bound tail*r/(1-r) with r measured from the last two windows and clipped to [1/4, 3/4] (hexabox_vop_lib.tail_bound_geom). Returns (worst, (i,k,seg)).""" worst = mp.mpf(0); wkey = None for (i, k), tab in sorted(self.NUM.items()): if i not in self.TOP: continue for s in range(self.NSEG): tail, prev = self.eng.coeff_tail(tab[s], window) b = VL.tail_bound_geom(tail, prev, window) if b > worst: worst, wkey = b, (i, k, s) return worst, wkey def _rat_tab(self, R): key = R.key() if key in self._RV: return self._RV[key] ncf = [_mpc2acb(mp.mpf(int(c.p))/int(c.q)) for c in reversed(R.n.coeffs())] dcf = [_mpc2acb(mp.mpf(int(c.p))/int(c.q)) for c in reversed(R.d.coeffs())] Z0 = flint.acb(0); out = [] for s in range(self.NSEG): row = [] for k in range(self.NC+1): z = self.eng.zga[s][k] num = Z0 for c in ncf: num = num*z + c den = Z0 for c in dcf: den = den*z + c row.append(num/den) out.append(row) self._RV[key] = out return out def _rad_tab(self, rad): if rad in self._RT: return self._RT[rad] rv = self.eng.rad_vals(rad) out = [[_mpc2acb(x) for x in row] for row in rv] self._RT[rad] = out return out def _tab_add(self, A, B): if A is None: return B if B is None: return A return [[A[s][k]+B[s][k] for k in range(self.NC+1)] for s in range(self.NSEG)] def _tab_mulrat(self, A, R): if A is None or R.is_zero(): return None rv = self._rat_tab(R) return [[A[s][k]*rv[s][k] for k in range(self.NC+1)] for s in range(self.NSEG)] def _tab_mulalg(self, A, R, rad): if A is None or R.is_zero(): return None rv = self._rat_tab(R) if rad: rt = self._rad_tab(rad) return [[A[s][k]*rv[s][k]*rt[s][k] for k in range(self.NC+1)] for s in range(self.NSEG)] return [[A[s][k]*rv[s][k] for k in range(self.NC+1)] for s in range(self.NSEG)] def _tab_int(self, A): if A is None: return None out = []; cum = flint.acb(0) for s in range(self.NSEG): g = [A[s][k]*self.eng.dza[s] for k in range(self.NC+1)] segcum = self.eng._cumcheb_acb(g) out.append([cum+segcum[k] for k in range(self.NC+1)]) cum += segcum[-1] return out def _tab_const(self, val): v = _mpc2acb(mp.mpc(val)) return [[v]*(self.NC+1) for _ in range(self.NSEG)] def _sun_tab(self, j, k): if (j, k) in self._ST: return self._ST[(j, k)] C = self.SUN_CF.get((j, k)) if not C: self._ST[(j, k)] = None; return None out = None for (tag, rad, word), r in C.items(): wv = self.eng.word_vals_acb(word) rv = self._rat_tab(r) cv = _mpc2acb(CONSTS[tag]) t2 = [[cv*rv[s][kk]*wv[s][kk] for kk in range(self.NC+1)] for s in range(self.NSEG)] out = self._tab_add(out, t2) self._ST[(j, k)] = out return out def _hyb_tab(self, key): """Node table of a CERTIFIED-NUMERIC kernel (disclosed; not exact).""" if key in self._HT: return self._HT[key] Pc, Qc = self.HYBRID[key] Pa = [_mpc2acb(c) for c in reversed(Pc)] Qa = [_mpc2acb(c) for c in reversed(Qc)] Z0 = flint.acb(0); out = [] for s in range(self.NSEG): row = [] for kk in range(self.NC+1): z = self.eng.zga[s][kk] num = Z0 for c in Pa: num = num*z + c den = Z0 for c in Qa: den = den*z + c row.append(num/den) out.append(row) self._HT[key] = out self.hybrid_used.append(key) self.log(f"[vop/T2] CERTIFIED-NUMERIC kernel {key} entered the recursion " f"(degQ={len(Qc)-1} stored fit — precision through it is capped)") return out def _getn_tab(self, j, k): if j in self.SUNset: return self._sun_tab(j, k) if k >= -1 else None return self.NUM.get((j, k)) def solve_lin_num(self, strat, Stabs, w0vals): if strat[0] == 'CONJ': _, radli, sub = strat invfac = R_ONE for li in radli: fpoly = LETTERS[li] invfac = invfac * Rat(FP([fpoly(fq(0))]), fpoly) Smod = [self._tab_mulalg(S, invfac, radli) if S is not None else None for S in Stabs] Z = self.solve_lin_num(sub, Smod, w0vals) rt = self._rad_tab(radli) return [([[z[s][k]*rt[s][k] for k in range(self.NC+1)] for s in range(self.NSEG)] if z is not None else None) for z in Z] if strat[0] == 'FULL': _, cols, rads = strat nb = len(cols) Rm = [[cols[c][a] for c in range(nb)] for a in range(nb)] Rinv, det = self.rat_matinv(Rm) Ys = [None]*nb for c in range(nb): Tc = None for b in range(nb): if Rinv[c][b].is_zero() or Stabs[b] is None: continue Tc = self._tab_add(Tc, self._tab_mulrat(Stabs[b], Rinv[c][b])) if rads[c] and Tc is not None: invfac = R_ONE for li in rads[c]: fpoly = LETTERS[li] invfac = invfac * Rat(FP([fpoly(fq(0))]), fpoly) Tc = self._tab_mulalg(Tc, invfac, rads[c]) Ic = self._tab_int(Tc) c0 = mp.mpc(0) for b in range(nb): if Rinv[c][b].is_zero(): continue rv0 = Rinv[c][b].ev_fq(0) c0 += mp.mpf(int(rv0.p))/int(rv0.q) * w0vals[b] inner = Ic if abs(c0) > mp.mpf(10)**-72: inner = self._tab_add(inner, self._tab_const(c0)) if inner is None: continue for a in range(nb): if Rm[a][c].is_zero(): continue Ys[a] = self._tab_add(Ys[a], self._tab_mulalg(inner, Rm[a][c], rads[c])) return Ys _, T, Tinv, Bm, f, sub = strat nb = len(T) TS = []; W0 = [] for c in range(nb): acc = None for b in range(nb): if Tinv[c][b].is_zero() or Stabs[b] is None: continue acc = self._tab_add(acc, self._tab_mulrat(Stabs[b], Tinv[c][b])) TS.append(acc) c0 = mp.mpc(0) for b in range(nb): if Tinv[c][b].is_zero(): continue rv0 = Tinv[c][b].ev_fq(0) c0 += mp.mpf(int(rv0.p))/int(rv0.q) * w0vals[b] W0.append(c0) wsub = self.solve_lin_num(sub, TS[f:], W0[f:]) w = [None]*nb for c in range(nb-f): w[f+c] = wsub[c] for a in range(f): integ = TS[a] for b in range(f, nb): if Bm[a][b].is_zero() or w[b] is None: continue integ = self._tab_add(integ, self._tab_mulrat(w[b], Bm[a][b])) wa = self._tab_int(integ) if abs(W0[a]) > mp.mpf(10)**-72: wa = self._tab_add(wa, self._tab_const(W0[a])) w[a] = wa Ys = [None]*nb for a in range(nb): for c in range(nb): if T[a][c].is_zero() or w[c] is None: continue Ys[a] = self._tab_add(Ys[a], self._tab_mulrat(w[c], T[a][c])) return Ys def run_numeric(self, KMAX=3): """Layers -4..KMAX of every block on the contour (Tier 2).""" self.reach_check(KMAX) for k in range(-4, KMAX+1): tL = time.time() for s in self.ORDER: rows = self.BLOCKS[s] if KMAX >= 3 and k == KMAX and not any(i in self.C21 for i in rows): continue self.reach_check(KMAX, k=k, rows=rows) Stabs = [] for i in rows: S = None for j in self.COLS.get(i, []): for m in self.SUPP.get((i, j), []): if m == 0 and j in rows: continue src = self._getn_tab(j, k-m) if src is None: continue S = self._tab_add(S, self._tab_mulrat(src, self.RatA[(i, j, m)])) for key in self.HYBRID: if key[0] != i: continue src = self._getn_tab(key[1], k-key[2]) if src is None: continue pt = self._hyb_tab(key) S = self._tab_add(S, [[src[s2][k2]*pt[s2][k2] for k2 in range(self.NC+1)] for s2 in range(self.NSEG)]) Stabs.append(S) w0 = [] for i in rows: x = self.V0.get((i, k)) w0.append(mp.mpc(0) if (x is None or abs(x) < mp.mpf(10)**-65) else x) if all(abs(x) < mp.mpf(10)**-72 for x in w0) and all(S is None for S in Stabs): continue strat = self.block_strategy(s) Ys = self.solve_lin_num(strat, Stabs, w0) for a, i in enumerate(rows): if Ys[a] is not None: self.NUM[(i, k)] = Ys[a] self.log(f"[vop/T2] layer {k}: " f"{sum(1 for (i2,k2) in self.NUM if k2==k)} rows ({time.time()-tL:.1f}s)") return self.NUM def endpoint(self, i, k): tab = self.NUM.get((i, k)) if tab is None: return None return _acb2mpc(tab[self.NSEG-1][self.NC]) def node_value(self, i, k, s, kk): """(z, n[i,k](z)) at contour node (segment s, node kk) — a live function value of the layered iterated integral at an interior path point.""" tab = self.NUM.get((i, k)) if tab is None: return None, None return self.eng.zg[s][kk], _acb2mpc(tab[s][kk]) def value_at(self, i, k, z): """n[i,k](z) at ANY point z on the spectral contour (real t outside the detour windows, or on the detours): barycentric Chebyshev interpolation of the node table — spectrally accurate for the analytic integrand, so precision grows with NC (up to the seed/kernel caps disclosed above). Raises ValueError if z is off the contour (domain limit).""" tab = self.NUM.get((i, k)) if tab is None: return None s, u = self.eng.locate(mp.mpc(z)) return self.eng.interp([_acb2mpc(x) for x in tab[s]], u) def real_domain(self): """Real path intervals covered by the contour (the detour windows in between are excluded — the domain limit for real-t evaluation).""" out = [] for (za, zb) in self.eng.segs: if abs(mp.im(za)) < 1e-12 and abs(mp.im(zb)) < 1e-12: out.append((float(mp.re(za)), float(mp.re(zb)))) return out