#!/usr/bin/env python3 r"""direct_linear_extract.py -- DIRECT LINEAR extraction of the threshold connection coefficients c_alpha at >= 100 working digits (python-flint acb balls at >= 1400 bits), printed beside the spectral projector's digits. STAMP (date -u at emission): 2026-09-10T13:58:46Z extends: the threshold-banana / banana-evaluate code set (rowscripts): the (1,1,1,9) K3 row-31 and the (1,1,1,1,16) CY3 row-32 operators OF RECORD (PF_119_EXACT.json / pf_cy3_data.json, theta_z form, read by object with sha256 pins in the fixtures), the multinomial-squared period series generated exactly here and checked against the operator by EXACT annihilation before anything is transported; the projector's digits are read from the record objects and printed beside (ROW31_ARB_VALUE_OF_RECORD.json, c32_record.json, GATE.json, LEGII_CY3.json). Nothing here overwrites or imports the served evaluators. reference: an independent reference derivation of 2026-09-10, on file with the authors (a 6x6 direct solve at 1400 bits reaching 414 digits) -- the claim tested here is the one the served record motivates the projector with: "the unipotent tower is ~1e4 larger than the fractional admixture, so a direct solve fishes the answer out of subtraction noise" (a working-precision cap of the dps-30 first pass, not a structural one). WHAT IT DOES (--family K3 | CY3, --direct-linear BITS): 1. c_n exactly (ints) from the multinomial-squared definition; the record operator L = sum_j P_j(z) theta_z^j (z = -1/s) applied to sum_n (-1)^n c_n z^n must give residual EXACTLY 0 over --ncheck coefficients (refused otherwise, exit 3); the threshold indicial polynomial and its exponents from the same operator. 2. the s-chart D-form p_k(s) (exact ints): theta_z = -s d/ds, s^degz P_j(-1/s), Stirling numbers of the second kind for theta_s^j -> s^i d^i/ds^i. 3. the seed jet [f, f', ..., f^(r-1)] at s_b (real, above the threshold (sum sqrt M)^2) from the exact series, as balls whose radius carries the RIGOROUS tail bound |c_n| <= thr^n. 4. Taylor-series transport of the jet along the record's path [s_b, s_b + i lift, s_dec + i lift, s_dec] (Im s > 0: the physical side) in acb at BITS, step |h| <= --frac x (distance to the nearest root of p_r), the recursion in MIDPOINT arithmetic (a ball recurrence of ~10^3 terms blows its radii up exponentially); the truncation tail of every step is estimated by the trailing-window geometric envelope (the served evaluators' form) and printed; the same transport is repeated at BITS/2 + 100 and the two-precision agreement is the error estimate. The CERTIFICATE of the transport is the separate ore_algebra oracle, not this file. 5. the local Frobenius basis at s = 0 EXACTLY over Q: every fractional branch Phi_alpha = s^alpha (1 + ...) and every unipotent block (exponent rho_0 of multiplicity mu: the mu solutions [eps^j] s^{rho_0+eps} sum_n a_n(eps) s^n, a_n(eps) truncated power series in eps with the resonant divisions done exactly -- a resonance whose numerator does not vanish to the denominator's order is refused by name); their jets at s_dec (principal branch s^alpha > 0, log s real). 6. ONE r x r linear solve (acb_mat.solve) of jet(varpi_0) = sum_k x_k jet(basis_k): the direct linear extraction; C_alpha = x_alpha; c_alpha = C_alpha e^{i pi alpha} (the t = e^{-i pi} s convention of record). 7. printed beside: the ball radius (digits), the two-precision agreement, the agreement with the record strings and with the closed form (calpha_rings structures = the analytic result of threshold_hankel_tail.py), the size of the unipotent coefficients relative to c_alpha (the "1e4 wall"), and the PROJECTOR'S digits of record read by object (Route A / Route B / the IBP-operator leg). Exit 0 iff every record-string agreement >= min(50, its cap) and the two-precision agreement >= 100 d; a named FAIL exits 1; a pin mismatch or a failed annihilation exits 3; usage 2. --selftest: the clean run at 1400 bits on both families, then the planted controls (a) the record operator's coefficient of z^1 in the theta^2 term shifted by +1 (the served evaluator's planted control "2,1,1") -> the exact annihilation must FAIL by name (exit 3); (b) one digit of a record string flipped -> that comparison must FAIL by name. USAGE python3 direct_linear_extract.py --family K3 --direct-linear 1400 [--sdec 1/4] [--sb 80] [--lift 10] python3 direct_linear_extract.py --family CY3 --direct-linear 1400 --sdec 1/4 python3 direct_linear_extract.py --selftest Self-contained: python3 + python-flint (+ mpmath for the record-string digit counts) + the fixtures file. """ import argparse import hashlib import json import math import os import subprocess import sys import time from fractions import Fraction as Fr from flint import acb, arb, acb_poly, acb_mat, fmpq, fmpz, ctx STAMP = "2026-09-10T13:58:46Z" EXIT_FAIL, EXIT_USAGE, EXIT_PIN = 1, 2, 3 FIXTURES_DEFAULT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "fixtures", "direct_linear_fixtures.json") FIXTURES_SHA256 = "c3bb23092d1367f9eb94f0b250c3a93752593b2ea4f391d11e9f04aa28908e22" FAMILIES = {"K3": (1, 1, 1, 9), "CY3": (1, 1, 1, 1, 16)} def sha256_of(path): h = hashlib.sha256() with open(path, "rb") as f: for ch in iter(lambda: f.read(1 << 20), b""): h.update(ch) return h.hexdigest() def stamp_now(): return subprocess.run(["date", "-u", "+%Y-%m-%dT%H:%M:%SZ"], capture_output=True, text=True, check=True).stdout.strip() # --------------------------------------------------------------------------- exact series + operator def multinomial_squared_series(msq, N): prod = [Fr(1)] for M in msq: ser = [Fr(M) ** k / Fr(math.factorial(k)) ** 2 for k in range(N + 1)] new = [Fr(0)] * (N + 1) for i, a in enumerate(prod): if a == 0: continue for j in range(N + 1 - i): new[i + j] += a * ser[j] prod = new out = [] for n in range(N + 1): c = prod[n] * Fr(math.factorial(n)) ** 2 assert c.denominator == 1 out.append(int(c)) return out def extend_by_recurrence(Pj, c_head, N): """c_n for n > len(c_head)-1 from the operator's recurrence on a_n = (-1)^n c_n (exact integer division, asserted), the leading coefficient sum_j P_{j,0} n^j never zero for n >= 1 beyond the MUM exponents.""" degz = max(len(co) - 1 for co in Pj.values()) a = [(-1) ** n * c_head[n] for n in range(len(c_head))] for n in range(len(a), N + 1): tot = 0 for j, co in Pj.items(): for k, p in enumerate(co): if p and k >= 1 and 0 <= n - k < n: tot += p * (n - k) ** j * a[n - k] lead = sum(co[0] * n ** j for j, co in Pj.items() if co and co[0]) assert lead != 0 and tot % lead == 0, f"recurrence at n={n}: non-integer step" a.append(-tot // lead) return [(-1) ** n * a[n] for n in range(len(a))] def annihilation_residuals(Pj, a, ncheck): """L = sum_j P_j(z) theta_z^j on sum_n a_n z^n: the z^n coefficient sum_j sum_k P_{jk} (n-k)^j a_{n-k}.""" degz = max(len(co) - 1 for co in Pj.values()) res = [] for n in range(degz, ncheck): tot = 0 for j, co in Pj.items(): for k, p in enumerate(co): if p and 0 <= n - k < len(a): tot += p * (n - k) ** j * a[n - k] res.append(tot) return res def s_chart_dform(Pj, order): """p_i(s) (exact ints, low -> high) with L = sum_i p_i(s) d^i/ds^i, from the theta_z form with z = -1/s.""" degz = max(len(co) - 1 for co in Pj.values()) # Q_j(s) = s^degz P_j(-1/s) = sum_k p_jk (-1)^k s^{degz-k} Q = {} for j, co in Pj.items(): q = [0] * (degz + 1) for k, p in enumerate(co): q[degz - k] += p * (-1) ** k Q[j] = q # Stirling numbers of the second kind: theta^j = sum_i S(j,i) s^i D^i ; theta_z^j = (-theta_s)^j S = [[0] * (order + 1) for _ in range(order + 1)] S[0][0] = 1 for n in range(1, order + 1): for k in range(1, n + 1): S[n][k] = k * S[n - 1][k] + S[n - 1][k - 1] p = [[0] * (degz + order + 1) for _ in range(order + 1)] for j, q in Q.items(): for i in range(j + 1): if S[j][i] == 0: continue for e, c in enumerate(q): if c: p[i][e + i] += c * (-1) ** j * S[j][i] # strip common power of s while all(pi[0] == 0 for pi in p): p = [pi[1:] for pi in p] return p def indicial_threshold(Pj, order): """The s-chart local form at the threshold s = 0: L[s^e] = sum_i R_i(e) s^{e+i} with R_i(x) = (-1)^{degz-i} sum_j P_{j,degz-i} (-x)^j (z = -1/s, theta_z = -theta_s); R_0 is the indicial polynomial. (The record's branch series a32_series / abranch are in the t = -s chart: a_n^(s) = (-1)^n a_n^(t).)""" degz = max(len(co) - 1 for co in Pj.values()) R = [] for i in range(degz + 1): coeffs = [Fr(0)] * (order + 1) for j, co in Pj.items(): if degz - i < len(co): c = co[degz - i] if c: # z^k = (-1)^k s^{-k} with k = degz - i (the s-chart, z = -1/s); theta_z = -theta_s coeffs[j] += Fr(c) * (-1) ** (degz - i) * (-1) ** j R.append(coeffs) # R_i(rho) = sum_j coeffs[j] rho^j : L[s^e] = sum_i R_i(e) s^{e+i} return R def poly_eval(coeffs, x): tot = Fr(0) for c in reversed(coeffs): tot = tot * x + c return tot def exponents_from_indicial(R0): """Rational roots with multiplicities of R_0(rho) (the ladder's exponents are rational); found by exact rational root search on the primitive integer polynomial, refusing if a non-rational factor remains.""" from math import gcd den = 1 for c in R0: den = den * c.denominator // gcd(den, c.denominator) P = [int(c * den) for c in R0] while P and P[-1] == 0: P.pop() roots = {} # candidates p/q with p | P[0]..., use the standard rational root theorem on the deflated polynomial def divisors(n): n = abs(n) out = set() for d in range(1, int(math.isqrt(n)) + 1): if n % d == 0: out.add(d); out.add(n // d) return out poly = list(P) changed = True while changed and len(poly) > 1: changed = False # strip rho = 0 roots if poly[0] == 0: roots[Fr(0)] = roots.get(Fr(0), 0) + 1 poly = poly[1:] changed = True continue for q in sorted(divisors(poly[-1])): for p in sorted(divisors(poly[0])): for sgn in (1, -1): r = Fr(sgn * p, q) if poly_eval([Fr(c) for c in poly], r) == 0: # deflate new = [] acc = Fr(0) for c in reversed(poly): acc = acc * r + c new.append(acc) new = [c for c in reversed(new[:-1])] # quotient coefficients (low -> high) # normalise to integers d2 = 1 for c in new: d2 = d2 * c.denominator // gcd(d2, c.denominator) poly = [int(c * d2) for c in new] roots[r] = roots.get(r, 0) + 1 changed = True break if changed: break if changed: break if len(poly) > 1: raise RuntimeError(f"indicial polynomial has a non-rational factor {poly}") return roots # --------------------------------------------------------------------------- eps-series helpers (Fraction) def ser_mul(a, b, K): out = [Fr(0)] * K for i, x in enumerate(a): if x == 0 or i >= K: continue for j, y in enumerate(b): if i + j >= K: break out[i + j] += x * y return out def ser_div(num, den, K, where=""): """num/den as eps-series mod eps^K, allowing a common eps-valuation (the log-free resonance).""" v = 0 while v < len(den) and den[v] == 0: v += 1 if v == len(den): raise RuntimeError("division by the zero series" + where) if any(num[i] != 0 for i in range(min(v, len(num)))): raise RuntimeError(f"resonance{where}: the numerator does not vanish to order eps^{v} of the denominator " "(a logarithmic resonance -> the eps^v lift is needed; not implemented)") n2 = list(num[v:]) + [Fr(0)] * v d2 = list(den[v:]) + [Fr(0)] * v out = [Fr(0)] * K inv0 = 1 / d2[0] for n in range(K): acc = n2[n] if n < len(n2) else Fr(0) for k in range(1, n + 1): if k < len(d2): acc -= d2[k] * out[n - k] out[n] = acc * inv0 return out def poly_at_shift(coeffs, x0, K): """R(x0 + eps) as an eps-series mod eps^K (exact).""" out = [Fr(0)] * K # Taylor: sum_j c_j (x0+eps)^j for j, c in enumerate(coeffs): if c == 0: continue # (x0+eps)^j = sum_i C(j,i) x0^{j-i} eps^i for i in range(min(j, K - 1) + 1): out[i] += c * math.comb(j, i) * x0 ** (j - i) return out def frobenius_block(R, rho0, mult, N, extra): """a_n(eps) mod eps^K (K = mult + extra) for s^{rho0+eps} sum a_n s^n, a_0 = 1: the exact recursion sum_i R_i(rho + n - i) a_{n-i} = 0 with the resonant divisions exact.""" K = mult + extra a = [[Fr(1)] + [Fr(0)] * (K - 1)] for n in range(1, N + 1): acc = [Fr(0)] * K for i in range(1, len(R)): if n - i < 0: break t = ser_mul(poly_at_shift(R[i], rho0 + n - i, K), a[n - i], K) acc = [x + y for x, y in zip(acc, t)] den = poly_at_shift(R[0], rho0 + n, K) a.append([-x for x in ser_div(acc, den, K, where=f" at n={n}, rho0={rho0}")]) return a def jets_block(a, rho0, mult, s_dec, order, prec): """jets (acb lists, k = 0..order-1) of y_j = [eps^j] s^{rho0+eps} sum a_n(eps) s^n, j < mult, at real s_dec.""" ctx.prec = prec sq = fmpq(s_dec.numerator, s_dec.denominator) sa = arb(sq) Ls = sa.log() K = len(a[0]) ex = [Ls ** i / arb(math.factorial(i)) for i in range(K)] res = [[None] * order for _ in range(mult)] for k in range(order): G = [arb(0)] * K for n, an in enumerate(a): # falling factorial (rho0+n+eps)^{(k)} as eps-series ffp = [Fr(1)] + [Fr(0)] * (K - 1) for i in range(k): base = rho0 + n - i new = [Fr(0)] * K for d in range(K): new[d] += ffp[d] * base if d + 1 < K: new[d + 1] += ffp[d] ffp = new pr = ser_mul(an, ffp, K) e = rho0 + n - k sp_ = sa ** arb(fmpq(e.numerator, e.denominator)) for d in range(K): if pr[d] != 0: G[d] += arb(fmpq(pr[d].numerator, pr[d].denominator)) * sp_ for j in range(mult): tot = arb(0) for i in range(j + 1): tot += ex[i] * G[j - i] res[j][k] = acb(tot) return res def frobenius_fractional(R, alpha, N): a = [Fr(1)] for n in range(1, N + 1): acc = Fr(0) for i in range(1, len(R)): if n - i < 0: break acc += poly_eval(R[i], alpha + n - i) * a[n - i] den = poly_eval(R[0], alpha + n) assert den != 0 a.append(-acc / den) return a def jets_fractional(a, alpha, s_dec, order, prec): ctx.prec = prec sa = arb(fmpq(s_dec.numerator, s_dec.denominator)) out = [] for k in range(order): tot = arb(0) for n, q in enumerate(a): e = alpha + n ff = Fr(1) for i in range(k): ff *= (e - i) cq = q * ff e2 = e - k tot += arb(fmpq(cq.numerator, cq.denominator)) * sa ** arb(fmpq(e2.numerator, e2.denominator)) out.append(acb(tot)) return out # --------------------------------------------------------------------------- transport in acb def _mid(x): return acb(x.real.mid(), x.imag.mid()) class Transport: def __init__(self, p_int, prec, sing=None): ctx.prec = prec self.r = len(p_int) - 1 self.P = [acb_poly([acb(fmpz(c)) for c in co]) for co in p_int] if sing is None: lead = acb_poly([acb(fmpz(c)) for c in p_int[-1]]) sing = [complex(z) for z in lead.roots()] self.sing = list(sing) def dist(self, s0): s0 = complex(s0) return min(abs(s0 - x) for x in self.sing) def step(self, s0, jet, s1, envelope_window=8): r = self.r h = acb(s1) - acb(s0) shift = acb_poly([acb(s0), acb(1)]) Q = [p(shift) for p in self.P] Qc = [[q[j] for j in range(q.degree() + 1)] for q in Q] rho = abs(complex(s1) - complex(s0)) / self.dist(s0) N = int(ctx.prec * math.log(2) / -math.log(rho)) + 40 u = [None] * (N + r + 1) fact = acb(1) for k in range(r): u[k] = _mid(jet[k]) / fact fact = fact * (k + 1) lead = Qc[r][0] for m in range(0, N + 1): acc = acb(0) for k in range(r + 1): qk = Qc[k] for j in range(min(m, len(qk) - 1) + 1): if k == r and j == 0: continue idx = m - j + k ff = 1 for i in range(k): ff *= (idx - i) acc += qk[j] * u[idx] * ff ff = 1 for i in range(r): ff *= (m + r - i) u[m + r] = _mid(-acc / (lead * ff)) # MIDPOINT recursion: a ball recurrence of ~10^3 terms blows its radii # up exponentially (the wrapping effect); the transport error is # ESTIMATED by the two-precision agreement, not certified -- the # certificate is the ore_algebra oracle (ore_transport_oracle.py) out = [] hp_abs = abs(complex(s1) - complex(s0)) for k in range(r): tot = acb(0) hp = acb(1) for n in range(0, N + r + 1 - k): ff = 1 for i in range(k): ff *= (n + k - i) tot += u[n + k] * ff * hp hp = hp * h # trailing-window geometric envelope of the truncated tail (the served evaluators' form), in arb tailmax = arb(0) habs = arb(hp_abs) for n in range(N + r + 1 - k - envelope_window, N + r + 1 - k): ff = 1 for i in range(k): ff *= (n + k - i) term = abs(u[n + k]) * arb(ff) * habs ** n tailmax = tailmax.max(term) if hasattr(tailmax, "max") else (term if term > tailmax else tailmax) env = tailmax * arb(rho) / arb(1 - rho) self.envelope_max = max(getattr(self, "envelope_max", 0.0), float(env)) out.append(_mid(tot)) return out def path(self, waypoints, jet, frac=0.4): cur = complex(waypoints[0]) curx = waypoints[0] nsteps = 0 for wp in waypoints[1:]: target = complex(wp) while True: d = self.dist(cur) rem = target - cur if abs(rem) <= frac * d: jet = self.step(curx, jet, wp) cur, curx = target, wp nsteps += 1 break nxt = cur + rem / abs(rem) * frac * d nxt = complex(round(nxt.real * 1024) / 1024, round(nxt.imag * 1024) / 1024) jet = self.step(curx, jet, nxt) cur, curx = nxt, nxt nsteps += 1 return jet, nsteps def seed_jet(c, thr, sb, order, prec): """exact partial sums as balls + the rigorous tail bound (|c_n| <= thr^n).""" ctx.prec = prec N = len(c) - 1 sbq = Fr(sb) out, tails = [], [] for k in range(order): tot = Fr(0) for n in range(N + 1): ff = 1 for i in range(k): ff *= (-n - i) tot += Fr(c[n] * ff) / sbq ** (n + k) rho = math.exp(k / (N + 1)) * float(thr / sbq) assert rho < 1 T = float(sbq) ** (-k) * (N + 1 + k) ** k * float(thr / sbq) ** (N + 1) / (1 - rho) b = arb(fmpq(tot.numerator, tot.denominator)) b = b.add_error(arb(T)) if hasattr(b, "add_error") else arb(b.mid(), T) out.append(acb(b)) tails.append(T) return out, tails def digits_between(a_str, b_str, cap): import mpmath as mp with mp.workdps(int(cap) + 20): a, b = mp.mpf(a_str), mp.mpf(b_str) if b == 0: return float(cap) if a == 0 else min(float(cap), float(-mp.log10(abs(a)))) d = abs(a - b) / abs(b) return min(float(cap), float(-mp.log10(d))) if d > 0 else float(cap) def run_family(name, bits, fixtures, sb=None, sdec=None, lift=None, frac=0.4, ncheck=200, plant_operator=None, string_flip=None, quiet=False): fam = fixtures["families"][name] msq = tuple(fam["msq"]) opsrc = fam["operator"]["source"] src_path = resolve_source(opsrc["path"]) if src_path: got = sha256_of(src_path) if got != opsrc["sha256"]: sys.stderr.write(f"REFUSED (exit 3): operator source {opsrc['path']} sha256 {got[:16]} != pin {opsrc['sha256'][:16]}\n") sys.exit(EXIT_PIN) Pj = {int(j): [int(x) for x in co] for j, co in fam["operator"]["theta_form_Pj"].items()} if plant_operator: j, k, dv = plant_operator Pj[j][k] += dv order = max(Pj) degz = max(len(co) - 1 for co in Pj.values()) m = [math.isqrt(M) for M in msq] thr = Fr(sum(m) ** 2) sb = Fr(sb if sb is not None else fam["path_of_record"]["sb"]) lift = Fr(lift if lift is not None else fam["path_of_record"]["lift"]) sdec = Fr(sdec if sdec is not None else fam["path_of_record"]["sdec"][0]) rep = {"family": name, "msq": list(msq), "bits": bits, "sb": str(sb), "lift": str(lift), "sdec": str(sdec), "frac": frac} t0 = time.time() # 1. exact annihilation Nser = int(bits * math.log(2) / math.log(float(sb / thr))) + 60 c_def = multinomial_squared_series(msq, ncheck + degz + 5) # the definition, exact, for the certificate a_alt = [(-1) ** n * c_def[n] for n in range(len(c_def))] res = annihilation_residuals(Pj, a_alt, ncheck) nz = sum(1 for x in res if x != 0) c = None if nz == 0: c = extend_by_recurrence(Pj, c_def, max(Nser, ncheck) + 5) # the CERTIFIED operator's own recurrence beyond assert c[:len(c_def)] == c_def rep["annihilation"] = {"ncheck": ncheck, "nonzero_residuals": nz, "verdict": "EXACT 0" if nz == 0 else "FAIL"} if not quiet: print(f"\n== {name} msq={list(msq)}: record operator order {order} degz {degz} ({opsrc['path']} {opsrc['sha256'][:16]}{'' if src_path else ' -- source not present here, the pinned copy in the fixtures used'}); exact annihilation over {ncheck} coefficients: {'residual EXACTLY 0' if nz == 0 else f'{nz} NONZERO residuals'}") if nz: rep["fails"] = [f"{name}: the operator does not annihilate the period series ({nz} nonzero residuals of {len(res)})"] return rep, EXIT_PIN R = indicial_threshold(Pj, order) exps = exponents_from_indicial(R[0]) rep["threshold_exponents"] = {str(k): v for k, v in sorted(exps.items())} if not quiet: print(f" threshold exponents (indicial R_0): {rep['threshold_exponents']}") p_int = s_chart_dform(Pj, order) # singular points in s: s = 0 (the threshold) and s = -1/z_k over the roots z_k of the (squarefree) leading # theta-coefficient P_order(z); double precision suffices for the step-size control import numpy as _np zroots = _np.roots(list(reversed([float(c) for c in Pj[order]]))) sing_s = [0j] + [complex(-1.0 / z) for z in zroots if abs(z) > 1e-300] rep["singular_points_s"] = [f"{z.real:.6g}{z.imag:+.6g}j" for z in sing_s] # 2-4. seed + transport at two precisions results = {} for prec in (bits, bits // 2 + 100): ctx.prec = prec T = Transport(p_int, prec, sing=sing_s) jet0, tails = seed_jet(c[:Nser + 1], thr, sb, order, prec) wps = [complex(float(sb), 0), complex(float(sb), float(lift)), complex(float(sdec), float(lift)), complex(float(sdec), 0)] # exact waypoints as acb from rationals wps_acb = [acb(arb(fmpq(sb.numerator, sb.denominator)), arb(0)), acb(arb(fmpq(sb.numerator, sb.denominator)), arb(fmpq(lift.numerator, lift.denominator))), acb(arb(fmpq(sdec.numerator, sdec.denominator)), arb(fmpq(lift.numerator, lift.denominator))), acb(arb(fmpq(sdec.numerator, sdec.denominator)), arb(0))] t1 = time.time() jet, nsteps = T.path(wps_acb, jet0, frac=frac) t_tr = time.time() - t1 # 5. local basis at s = 0 Nthr = int(prec * math.log(2) / -math.log(float(sdec) / 4.0)) + 40 # the nearest true singularity of the ladder at s = 4 cols, labels = [], [] for rho, mult in sorted(exps.items()): if rho.denominator == 1: a = frobenius_block(R, rho, mult, Nthr, extra=2) J = jets_block(a, rho, mult, sdec, order, prec) for j in range(mult): cols.append(J[j]); labels.append(f"unip rho={rho} j={j}") else: a = frobenius_fractional(R, rho, Nthr) cols.append(jets_fractional(a, rho, sdec, order, prec)); labels.append(f"frac alpha={rho}") assert len(cols) == order, (len(cols), order) A = acb_mat(order, order) for col in range(order): for row in range(order): A[row, col] = cols[col][row] b = acb_mat(order, 1) for row in range(order): b[row, 0] = jet[row] x = A.solve(b) out = {} for i, lab in enumerate(labels): xi = x[i, 0] if lab.startswith("frac"): al = Fr(lab.split("=")[1]) ph = (acb(0, 1) * acb.pi() * acb(fmpq(al.numerator, al.denominator))).exp() cval = xi * ph out[str(al)] = {"C_principal": xi, "c": cval} else: out[lab] = {"x": xi} results[prec] = {"out": out, "nsteps": nsteps, "transport_s": round(t_tr, 2), "Nser": Nser, "Nthr": Nthr, "seed_tails": [f"{t:.2e}" for t in tails], "truncation_envelope_max": f"{getattr(T, 'envelope_max', 0.0):.2e}"} if not quiet: print(f" prec {prec} bits: transport {nsteps} steps {t_tr:.1f} s (Nser {Nser}, Nthr {Nthr}); seed tail bounds {results[prec]['seed_tails']}; largest per-step truncation envelope {results[prec]['truncation_envelope_max']} (midpoint arithmetic: an estimate, not a certificate)") hi, lo = bits, bits // 2 + 100 fails = [] rep["coefficients"] = {} ndig = int(hi * 0.30103) - 5 for key in results[hi]["out"]: if "C_principal" not in results[hi]["out"][key]: continue ch, cl = results[hi]["out"][key]["c"], results[lo]["out"][key]["c"] re_h = ch.real.str(ndig, radius=False) im_h = ch.imag.str(20, radius=False) rad = float(ch.rad()) mag = abs(complex(ch)) ball_d = (-math.log10(rad / mag) if (rad > 0 and mag > 0) else float("inf")) two_prec = digits_between(re_h, cl.real.str(int(lo * 0.30103) - 5, radius=False), lo * 0.30103) entry = {"c_re": re_h, "c_im": im_h, "solve_rounding_radius": f"{rad:.3e}", "solve_rounding_digits": round(ball_d, 1) if ball_d != float("inf") else "inf", "two_precision_d": round(two_prec, 1), "note": "the transport is midpoint arithmetic; the digits are the two-precision agreement and the agreement with the records / the closed form"} # unipotent size ratio (the "1e4 wall") unip_max = max(abs(complex(v["x"])) for k2, v in results[hi]["out"].items() if "x" in v) entry["unipotent_over_c_ratio"] = f"{unip_max / abs(complex(ch)):.3e}" if not quiet: print(f" c_{{{key}}} (direct linear, {hi} bits) = {re_h[:80]}... ; Im {im_h}; solve rounding radius {rad:.2e}; two-precision ({hi} vs {lo} bits) {two_prec:.1f} d; |unipotent|/|c| = {entry['unipotent_over_c_ratio']}") if two_prec < 100: fails.append(f"{name} c_{{{key}}}: two-precision agreement {two_prec:.1f} d < 100") entry["vs_records"] = [] for rs in fam.get("record_strings", {}).get(key, []): s = rs["value"] if string_flip and string_flip == (name, key, rs["label"]): s = flip_digit(s, 30) cap = min(int(rs["certified_digits"]), ndig) d = digits_between(re_h, s, cap) bar = min(50, cap - 2) ok = d >= bar entry["vs_records"].append({"label": rs["label"], "certified_digits": rs["certified_digits"], "agree_d": round(d, 1), "cap": cap, "bar": bar, "PASS": ok, "projector_note": rs.get("note")}) if not quiet: print(f" vs {rs['label']} ({rs['certified_digits']} d; {rs['source']['sha256'][:16]}): {d:.1f} d (cap {cap}) {'PASS' if ok else 'FAIL'} -- {rs.get('note','')}") if not ok: fails.append(f"{name} c_{{{key}}} vs {rs['label']}: {d:.1f} d < bar {bar}") cf = fam.get("closed_forms", {}).get(key) if cf: import mpmath as mp with mp.workdps(ndig + 20): val = eval(cf["mpmath_expr"], {"mp": mp, "pi": mp.pi, "sqrt": mp.sqrt, "gamma": mp.gamma, "mpf": mp.mpf}) d = digits_between(re_h, mp.nstr(val, ndig + 10), ndig) entry["vs_closed_form"] = {"expr": cf["display"], "agree_d": round(d, 1), "cap": ndig} if not quiet: print(f" vs the closed form {cf['display']} (threshold_hankel_tail.py, analytic): {d:.1f} d (cap {ndig} = the working precision)") rep["coefficients"][key] = entry rep["projector_of_record"] = fam.get("projector_of_record") if not quiet and fam.get("projector_of_record"): print(" the spectral projector's digits of record, by object:") for line in fam["projector_of_record"]: print(" -", line) rep["wall_s"] = round(time.time() - t0, 1) rep["fails"] = fails return rep, (EXIT_FAIL if fails else 0) def resolve_source(path): """A fixture source path is relative to the project root (BOOTSTRAP_ROOT, else the tree above this file's tools/ or blog/ dir); returns the existing absolute path or None (the pinned string is then used as is).""" if os.path.isabs(path): return path if os.path.exists(path) else None roots = [os.environ.get("BOOTSTRAP_ROOT", "")] here = os.path.dirname(os.path.abspath(__file__)) for _ in range(10): roots.append(here) here = os.path.dirname(here) for r in roots: if r and os.path.exists(os.path.join(r, path)): return os.path.join(r, path) return None def flip_digit(s, k): out = list(s) seen = 0 for i, ch in enumerate(out): if ch.isdigit() and (seen or ch != "0"): seen += 1 if seen == k: out[i] = str((int(ch) + 1) % 10) return "".join(out) raise ValueError("string too short") def load_fixtures(path): if not os.path.exists(path): sys.stderr.write(f"fixtures {path} missing (exit 3)\n") sys.exit(EXIT_PIN) got = sha256_of(path) if got != FIXTURES_SHA256: # 2026-09-11: unconditional (the self-arming template clause removed with its siblings') sys.stderr.write(f"REFUSED (exit 3): fixtures sha256 {got[:16]} != pin {FIXTURES_SHA256[:16]}\n") sys.exit(EXIT_PIN) return json.load(open(path)), got def main(argv=None): ap = argparse.ArgumentParser(description="direct linear extraction of c_alpha at >= 1400 bits (see the docstring)") ap.add_argument("--family", default="K3") ap.add_argument("--direct-linear", type=int, default=None, metavar="BITS", help="working precision in bits (>= 1400)") ap.add_argument("--sb", default=None) ap.add_argument("--sdec", default=None) ap.add_argument("--lift", default=None) ap.add_argument("--frac", type=float, default=0.4) ap.add_argument("--ncheck", type=int, default=200) ap.add_argument("--fixtures", default=FIXTURES_DEFAULT) ap.add_argument("--selftest", action="store_true") ap.add_argument("--json", default=None) a = ap.parse_args(argv) print(f"direct_linear_extract.py STAMP {STAMP}; python-flint {__import__('flint').__version__}") fixtures, fsha = load_fixtures(a.fixtures) print(f"fixtures {a.fixtures} sha256 {fsha[:16]}...") report = {"stamp": STAMP, "fixtures_sha256": fsha} rc = 0 if a.selftest: bits = a.direct_linear or 1400 print(f"\n== SELFTEST at {bits} bits") clean = {} for name in ("K3", "CY3"): rep, r = run_family(name, bits, fixtures, ncheck=a.ncheck, frac=a.frac, quiet=True) clean[name] = rep print(f" clean {name}: {'PASS' if r == 0 else 'FAIL ' + str(rep.get('fails'))}; " + "; ".join(f"c_{{{k}}} two-prec {v['two_precision_d']} d, vs records {[x['agree_d'] for x in v['vs_records']]}, vs closed {v.get('vs_closed_form',{}).get('agree_d')}" for k, v in rep["coefficients"].items())) if r: rc = EXIT_FAIL # (a) planted operator coefficient (the served control 2,1,1): exact annihilation must FAIL rep_a, r_a = run_family("K3", 600, fixtures, ncheck=a.ncheck, plant_operator=(2, 1, 1), quiet=True) ok_a = (r_a == EXIT_PIN and rep_a["annihilation"]["verdict"] == "FAIL") print(f" planted operator (theta^2, z^1, +1): {'REFUSED by name (annihilation FAIL, rc 3)' if ok_a else 'NOT as expected'}") # (b) planted record-string digit rep_b, r_b = run_family("K3", 600, fixtures, ncheck=a.ncheck, string_flip=("K3", "3/2", "ROW31_ARB_VALUE_OF_RECORD.value_string"), quiet=True) ok_b = (r_b == EXIT_FAIL and any("ROW31_ARB_VALUE_OF_RECORD.value_string" in f for f in rep_b["fails"]) and len(rep_b["fails"]) == 1) print(f" planted record-string digit (K3, digit 30): {'FAILED by name on that string only' if ok_b else 'NOT as expected: ' + str(rep_b.get('fails'))}") report["selftest"] = {"clean": clean, "planted_operator": {"rc": r_a, "annihilation": rep_a["annihilation"], "PASS": ok_a}, "planted_string": {"rc": r_b, "fails": rep_b.get("fails"), "PASS": ok_b}} ok = (rc == 0) and ok_a and ok_b report["selftest"]["PASS"] = ok print("SELFTEST:", "PASS" if ok else "FAIL") rc = 0 if ok else EXIT_FAIL else: if a.direct_linear is None: ap.error("--direct-linear BITS (>= 1400) or --selftest") if a.direct_linear < 1400: ap.error("--direct-linear must be >= 1400 bits (the >= 100 working digits this extraction targets, with margin)") if a.family not in FAMILIES: ap.error("--family K3 | CY3") rep, rc = run_family(a.family, a.direct_linear, fixtures, sb=a.sb, sdec=a.sdec, lift=a.lift, frac=a.frac, ncheck=a.ncheck) report["run"] = rep print("\nVERDICT:", "PASS" if rc == 0 else "FAIL") for f in rep.get("fails", []): print(" FAIL:", f) if a.json: with open(a.json, "x") as f: json.dump(report, f, indent=1, default=str) print("report written:", a.json) return rc if __name__ == "__main__": sys.exit(main())