The validation package of Riemann Hypothesis over Finite Holographic Substrate: one script, rh.py, 6 blocks, Python with numpy, scipy, mpmath, sympy, matplotlib, one check per predicate of the paper's predicate ledger — 110 of 110 checks passing. The script is rendered below with an anchor per line; its markers (# 20:B4 (p20011): the predicate and its key) are the checks the ledger predicates cite, and the block index opens each block at its definition. The notebook frc-20-rh.ipynb runs one cell per predicate on Google Colab (its generator make_notebook.py). Installed from this site as a named requirement: pip install frc-20-rh --find-links https://finitering.space/pkg/.
Blocks: A (16) · B · C · D (4) · E · T (1) — the predicates each decides in parentheses.
1"""2rh.py — the validation package of "Riemann Hypothesis over Finite Holographic Substrate" (Akhtman & Voether, 2026; Preprints3doi 10.20944/preprints202606.0768.v1, whose first version is "Riemann Hypothesis over the Holographic Substrate: a Finite-Field4Dictionary and the Screen Reading"), the paper 20-rh of the FRC corpus5(finite-ring-space/src/20-rh); one script since 24 September 2026 (the registry rhcommon.py and the seven block scripts6a_shell, b_deframe, c_primeside, c_gue, c_chi, d_classification, e_resonance merged, the three prime-side scripts as one block C).7========================================================================================================================89One script, six blocks, 110 checks — 66 EXACT (integer-pinned), 41 [approx] (a floating-point observation against the value10the paper states, to a stated tolerance), 3 [chart] (a continuum reading of a finite object; notation.tex convention 6) —11Python with numpy, scipy, mpmath, sympy and matplotlib (the blocks redraw their diagnostic figures into figures/ and the paper's12three numerical figures, at the page width, into figures/paper/). Each13check names the predicate(s) of the paper's ledger it witnesses (LEDGER; predicates cited as 20:XN), and the deciding check's14`# 20:XN (<key>)` marker is what the ledger's source column links (PREDICATES; finitering.space/src/20-rh/#<key>). Since the15cut of 15 October 2026 (T05 of reports/blueprint-20261013.md) the checks of the retired rows (the staircase, the horizon-scale16resolution, the semiclassical cutoff, the prime-side spectrum, the GUE statistic, the χ-twist, ...) still run as illustrations,17marked "illustration: 20:XN retired 2026-10-15", and witness no row; since T11 of the same day (every theorem row with a python18and a Lean witness) so do the checks of the resonance, the bare main sum, the raw secular condition, the Davenport–Heilbronn19comb and the pole-corrected count (20:C2, D4, E6, F3–F5, F8). Master-ledger predicates reached through the paper20predicates: 00:D11 (the shell theorem, 20:E12) and 00:D12 (the cells and the hosted theorem, 20:D5, 20:F2).2122Block A the shell theorem and the exact shell arithmetic, the closure's instances, completeness (EXACT; 20:A2, A3, A4, B1–B10, E1, E12, E13, F7)23Block B the de-framing dictionary ([chart], [approx]; illustrations since the cut: the density, the residue, the bare main sum)24Block C the spectrum from the prime side, ζ never evaluated ([approx]; illustrations since the cut: the secular condition, the GUE statistic, the χ-twist)25Block D Turing's count on the shell, the phantom, the Davenport–Heilbronn controls, the smoothed explicit formula, the strip,26 the finite control, the cells by the argument principle ([approx], EXACT; 20:D5, D8, D9, F1; the Davenport–Heilbronn27 comb and the pole-corrected count illustrations since the cut)28Block E resonance (EXACT; illustrations since the cut: the resonance, the antipode, the horizon-scale resolution, the flatness)29Block T Turing's certificate over the laboratory Carriers' cycles: Ω = 233, 30089, 2408561, every cell clean ([approx] in30 double precision with explicit bounds; 20:F2; the largest read from its embedded record, RH_CERT_FULL=1 recomputes it)3132Everything the blocks share (the former rhcommon.py):33 * the von Mangoldt comb Λ(n) (sieved) and its raised-cosine taper in log n;34 * the archimedean scale-phase θ(T) from log Γ (never from ζ);35 * the raw comb count Ñ_N(T) = θ(T)/π + 1 + S_comb(T) (Proposition turing; Obs. trace) and the pole-corrected36 count N̂_N(T) = Ñ_N(T) − (1/π) Im[Π_N(½+iT) + log((s−1)/s)] (Proposition combformula), Π_N the explicit37 pole term of the smoothed explicit formula — the raw count carries it and has no limit in N, the corrected38 count settles;39 * the secular roots (half-integer crossings of N̂_N, or of Ñ_N for the raw reading), the colleague matrix40 and the comb-built Jacobi matrix (Obs. matrix, Def. jacobi);41 * the first Riemann heights, used only as validation markers — never as inputs;42 * the PASS/FAIL registry that every block reports into (results.json).4344 python3 rh.py every block in order, results.json written; exit 1 if a check fails (≈ 8 min: the 10⁸ comb of45 block C and the 8×10⁷ comb of block D are the slow parts; RH_FAST=1 runs reduced depths in ≈ 3 min,46 and the records then differ from the paper's)47 python3 rh.py A E the blocks named (A, B, C, D, E, T); no results.json48 from frc_20_rh import predicate; predicate("20:F1") one predicate: the blocks of the checks citing it run once per session49"""50import os, json, sys, math, cmath, time51from math import gcd52from fractions import Fraction53import numpy as np54import sympy as sp55import mpmath as mp56from scipy.special import loggamma, lambertw57from scipy.optimize import brentq58import matplotlib59matplotlib.use("Agg")60import matplotlib.pyplot as plt6162SCRIPT = os.path.splitext(os.path.basename(__file__))[0] # "rh": the one script, the name results.json and the site pages carry63FAST = os.environ.get("RH_FAST", "0") == "1" # reduced depths for a quick run (the records then differ from the paper's)64FIGDIR = os.environ.get("RH_FIGDIR", "figures")65os.makedirs(FIGDIR, exist_ok=True)66os.makedirs(f"{FIGDIR}/paper", exist_ok=True)6768# ----------------------------------------------------------------------------- the paper's figures69# The blocks draw their diagnostic figures (titles carrying the check ids) and, from the same arrays, the paper's figures:70# drawn at the text width of the manuscript (452.97 pt = 6.29 in, a4 with 1 in margins, 11 pt), so that \includegraphics scales71# nothing and the 8–9 pt type prints as drawn; no titles (the captions carry them), panel letters above the axes (16 October 2026).72TEXTWIDTH_IN = 452.97 / 72.2773PAPER_RC = {"font.size": 8.5, "axes.labelsize": 9, "axes.titlesize": 9, "xtick.labelsize": 8, "ytick.labelsize": 8,74 "legend.fontsize": 7.8, "legend.frameon": True, "legend.framealpha": 0.92, "legend.edgecolor": "0.8",75 "axes.linewidth": 0.7, "xtick.major.width": 0.6, "ytick.major.width": 0.6, "lines.linewidth": 1.1,76 "axes.grid": False, "pdf.fonttype": 42, "mathtext.fontset": "cm", "font.family": "DejaVu Sans"}77INK, RED, BLUE, ORANGE, NAVY, CRIMSON, GRAY = "#1a1a1a", "#c0392b", "#1f5fbf", "#e67e22", "#10325f", "#c8102e", "0.40"7879def panel_letter(ax, s):80 """The panel letter above the axes' top-left corner, as the captions cite it."""81 ax.text(-0.02, 1.03, s, transform=ax.transAxes, ha="right", va="bottom", fontsize=9)8283def paper_figure(name, draw, height):84 """Draw one paper figure under PAPER_RC at the text width into FIGDIR/paper/<name>.pdf (and .png), restoring the rc."""85 with plt.rc_context(PAPER_RC):86 fig, ax = plt.subplots(1, 2, figsize=(TEXTWIDTH_IN, height), gridspec_kw={"width_ratios": [1.55, 1.0]} if name != "fig_dh" else None)87 draw(fig, ax)88 fig.tight_layout(w_pad=1.2)89 fig.savefig(f"{FIGDIR}/paper/{name}.pdf", bbox_inches="tight", pad_inches=0.02)90 fig.savefig(f"{FIGDIR}/paper/{name}.png", dpi=220, bbox_inches="tight", pad_inches=0.02)91 plt.close(fig)92 print(f" wrote {FIGDIR}/paper/{name}.pdf")9394# ----------------------------------------------------------------------------- registry95RESULTS = []9697# The paper's predicate ledger (20-rh Appendix A, predicates cited as 20:XN): the predicate(s) each98# check witnesses. The ledger's source column cites these check ids in return.99LEDGER = {100 "A1": "20:B4", "A2": "20:B5", "A3": "20:B6", "A4": "20:B8", "A5": "20:B9", "A6": "20:B10", "A6b": "20:B10", "A7": "20:E1, 20:E12, 20:E13",101 "A7b": "20:E1, 20:E12, 20:E13", "A8": "", "A9": "20:B7", "A10": "20:B1", "A11": "20:B2", "A12": "20:A2, 20:A4, 20:F7", "A13": "20:A3",102 "B1": "", "B2": "", "B3": "20:B8", "B4": "",103 "C1": "", "C2": "", "C2b": "", "C2c": "", "C3": "", "C4": "", "C5": "", "C6": "", "C6b": "",104 "C7a": "", "C7b": "", "C7c": "",105 "D1": "20:F1", "D2a": "", "D2b": "", "D2c": "", "D2d": "", "D2e": "",106 "D2f": "", "D2g": "", "D3": "", "D4": "20:F1",107 "D6a": "20:D8", "D6b": "20:D8", "D6c": "20:D9", "D7a": "20:D5", "D7b": "20:D5",108 "E1a": "", "E1b": "", "E1c": "", "E2a": "", "E2b": "", "E3": "", "E4": "",109 "T1": "20:F2", "T2": "20:F2", "T3": "20:F2",110}111112BLOCK = {"A": "the shell theorem and the exact shell arithmetic (a_shell.py until 24 September 2026)", # check-id prefix -> the block (the function block_<letter> below)113 "B": "the de-framing dictionary: zero density, the reconstruction residue, the horizon-length main sum (b_deframe.py)",114 "C": "the spectrum from the prime side with ζ never evaluated, the GUE statistics of the target spectrum, the χ-twisted comb (c_primeside.py, c_gue.py, c_chi.py)",115 "D": "Turing's count on the shell, the phantom, the Davenport–Heilbronn controls, the smoothed explicit formula (d_classification.py)",116 "E": "resonance, antipode, horizon-scale resolution, square-root cancellation (e_resonance.py)",117 "T": "Turing's certificate over the laboratory Carriers' cycles: every cell clean (Thm. hosted; 15 October 2026)"}118119# the deciding check of each python-witnessed predicate: the first check the ledger's source column names (the other checks that120# cite it are corroboration, listed by predicate() from the records)121PREDICATES = {122 "20:B1": "A10", "20:B2": "A11", "20:B4": "A1", "20:B5": "A2", "20:B6": "A3", "20:B7": "A9", "20:B8": "A4", "20:B9": "A5",123 "20:B10": "A6", "20:A2": "A12", "20:A3": "A13", "20:A4": "A12", "20:F7": "A12",124 "20:E1": "A7", "20:E12": "A7", "20:E13": "A7b", "20:F1": "D1", "20:D8": "D6a", "20:D9": "D6c", "20:D5": "D7a", "20:F2": "T1",125}126_RAN = set() # blocks already run in this session (predicate() runs each once)127128def check(pid, label, ok, detail="", kind="EXACT"):129 """Record one predicate check. pid = package check id (A1, C3, ...); LEDGER[pid] = the paper's ledger predicate(s)."""130 ok = bool(ok)131 rows = LEDGER.get(pid, "")132 RESULTS.append({"id": pid, "rows": rows, "block": pid[0], "script": SCRIPT, "label": label, "ok": ok, "detail": detail, "kind": kind})133 print(f" [{'PASS' if ok else 'FAIL'}] {pid:4s} {kind:7s} [{rows}] {label}" + (f" -- {detail}" if detail else ""))134 return ok135136def summary(write=True):137 n_ok = sum(r["ok"] for r in RESULTS)138 print(f"\nSUMMARY: {n_ok}/{len(RESULTS)} checks passed" + ("" if n_ok == len(RESULTS) else " <-- FAILURES"))139 if write:140 with open("results.json", "w") as f:141 json.dump(RESULTS, f, indent=1)142 return n_ok == len(RESULTS)143144def markers():145 """predicate label -> (script file, line) of its marker `# <paper>:<label> (<key>)`: the line of the check that decides it."""146 import re147 out = {}148 for i, line in enumerate(open(os.path.abspath(__file__), encoding="utf-8"), 1):149 m = re.match(r"\s*# ((?:\d+:[A-Z]+\d+[a-z]?(?: \(p\d{5}\))?)(?:, \d+:[A-Z]+\d+[a-z]?(?: \(p\d{5}\))?)*)\s*$", line)150 if m:151 for lab in re.findall(r"\d+:[A-Z]+\d+[a-z]?", m.group(1)): out.setdefault(lab, (SCRIPT + ".py", i))152 return out153154def _run_block(letter):155 if letter not in _RAN:156 globals()[f"block_{letter}"](); _RAN.add(letter)157158def predicate(label, lines=14):159 """Verify one ledger predicate: run the blocks of the checks that cite it (each once per session), print the deciding check's160 source (from its marker) and every record that cites the predicate, and return True iff all pass."""161 pid = PREDICATES.get(label)162 if pid is None:163 print(f"{label}: no python witness (see the predicate's Lean witness or its source)"); return None164 citing = {i[0] for i, rows in LEDGER.items() if label in [t.strip() for t in rows.split(",")]}165 for b in sorted({pid[0]} | citing): _run_block(b) # the deciding block and every block whose checks cite the predicate166 mk = markers().get(label)167 if mk:168 src = open(os.path.abspath(__file__), encoding="utf-8").read().split("\n")169 print(f"— {mk[0]}:{mk[1]} (the check that decides {label}: {pid})")170 for j in range(mk[1] - 1, min(mk[1] - 1 + lines, len(src))): print(f"{j + 1:5d} {src[j]}")171 recs = [r for r in RESULTS if label in [t.strip() for t in r["rows"].split(",")]]172 ok = all(r["ok"] for r in recs)173 for r in recs:174 role = "(deciding)" if r["id"] == pid else "(corroborating)"175 print(f" [{'PASS' if r['ok'] else 'FAIL'}] {r['id']:4s} {role:16s} {r['detail'][:150]}")176 print(f"{label}: {'VERIFIED' if ok and recs else 'FAILED'} — {len(recs)} record(s)")177 return ok178179def verify_all():180 """Run every block (those already run in this session are not re-run) and print the summary; True iff every check passed."""181 for b in sorted(BLOCK): _run_block(b)182 return summary(write=False)183184# ----------------------------------------------------------------------------- the shared primitives (rhcommon.py until 24 September 2026)185# ----------------------------------------------------------------------------- primes / comb186def sieve_primes(N):187 s = np.ones(N + 1, dtype=bool); s[:2] = False188 for i in range(2, int(N ** 0.5) + 1):189 if s[i]:190 s[i * i::i] = False191 return np.nonzero(s)[0]192193def von_mangoldt_support(N):194 """Prime powers n = p^k <= N with Λ(n) = log p. Returns (n, Λ(n)) as float arrays, sorted."""195 primes = sieve_primes(N)196 ns, lam = [primes.astype(np.float64)], [np.log(primes)]197 k = 2198 while True:199 pk = primes[primes <= N ** (1.0 / k)].astype(np.float64) ** k200 if len(pk) == 0:201 break202 ns.append(pk); lam.append(np.log(pk ** (1.0 / k))); k += 1203 n_all = np.concatenate(ns); lam_all = np.concatenate(lam)204 order = np.argsort(n_all)205 return n_all[order], lam_all[order]206207class Comb:208 """The tapered von Mangoldt comb to depth N.209210 S_comb(T) = -(1/π) Σ_{n<=N} Λ(n)/(√n log n) · w(n) · sin(T log n),211 w(n) = ½(1 + cos(π log n / log N)) — the raised-cosine (Cesàro) taper of Obs. primespec.212 """213 def __init__(self, N, lam_fn=None):214 self.N = int(N)215 if lam_fn is None:216 n, lam = von_mangoldt_support(self.N)217 else: # a generalized comb (Davenport–Heilbronn)218 n, lam = lam_fn(self.N)219 self.n, self.lam = n, lam220 self.logn = np.log(n)221 self.w = 0.5 * (1.0 + np.cos(np.pi * self.logn / math.log(self.N)))222 self.amp = self.lam / (np.sqrt(n) * self.logn) * self.w223224 NBLK = 1 << 18 # comb terms per block (memory: NBLK × 64 doubles ≈ 128 MB at most)225226 def S(self, T):227 """S_comb at scalar or array T, chunked over T and over the comb (a scalar in gives a float out)."""228 T = np.atleast_1d(np.asarray(T, dtype=float))229 out = np.zeros(len(T))230 for s in range(0, len(T), 64):231 blk = T[s:s + 64]232 for a in range(0, len(self.logn), self.NBLK):233 out[s:s + 64] -= (self.amp[a:a + self.NBLK, None] * np.sin(np.outer(self.logn[a:a + self.NBLK], blk))).sum(0)234 out /= math.pi235 return out if out.size > 1 else float(out[0])236237 def sigma(self, gamma):238 """|Σ_N(γ)| = |Σ Λ(n) n^{-1/2-iγ} w(n)|, the scale-spectrum of Obs. primespec (scalar in, float out)."""239 gamma = np.atleast_1d(np.asarray(gamma, float))240 coef = self.lam * self.w / np.sqrt(self.n)241 out = np.empty(len(gamma))242 for s in range(0, len(gamma), 64):243 blk = gamma[s:s + 64]244 acc = np.zeros(len(blk), dtype=complex)245 for a in range(0, len(self.logn), self.NBLK):246 acc += (coef[a:a + self.NBLK, None] * np.exp(-1j * np.outer(self.logn[a:a + self.NBLK], blk))).sum(0)247 out[s:s + 64] = np.abs(acc)248 return out if out.size > 1 else float(out[0])249250def theta(T):251 """Riemann–Siegel θ(T) = Im log Γ(¼ + iT/2) − (T/2) log π — the archimedean phase, from Γ only."""252 T = np.asarray(T, dtype=float)253 return np.imag(loggamma(0.25 + 0.5j * T)) - 0.5 * T * math.log(math.pi)254255def count_raw(comb, T, const=1.0, theta_fn=theta):256 """Ñ_N(T) = θ(T)/π + const + S_comb(T): the raw comb count (Proposition turing; uncertified)."""257 T = np.atleast_1d(np.asarray(T, float))258 val = theta_fn(T) / math.pi + const + comb.S(T)259 return val if val.size > 1 else float(val[0])260261# ----------------------------------------------------------------------------- the pole term (Prop. combformula)262def pole_term(N, T):263 """Π_N(½ + iT) = ∫_0^L W(u/L) (e^{(½−iT)u} − e^{−(½+iT)u}) du/u, L = log N — the pole term of the smoothed264 explicit formula for the tapered log-ζ series, Σ_{n≤N} Λ(n) W(log n/L) n^{−s}/log n = Π_N(s) + J_N(s) (Prop.265 combformula (i), under the hypothesis of no zero to the right of the line), with J_N(s) → log((s−1)ζ(s)/s) for266 s not a zero (Prop. combformula (iv)). Written as ∫_0^L g(u) e^{−iTu} du, g(u) = W(u/L)·2 sinh(u/2)/u, and integrated with267 QUADPACK's oscillatory weights. Returns the complex value. Its size is (π²/2) √N/(|½−iT| L)³ (1 + O(1/L)):268 the raw count carries it and has no limit as N → ∞."""269 from scipy.integrate import quad270 L = math.log(N)271 g = lambda u: (0.5 * (1.0 + math.cos(math.pi * u / L))) * (2.0 * math.sinh(0.5 * u) / u if u > 1e-12 else 1.0)272 re = quad(g, 0.0, L, weight="cos", wvar=T, limit=1000)[0]273 im = -quad(g, 0.0, L, weight="sin", wvar=T, limit=1000)[0]274 return complex(re, im)275276def pole_count_term(N, T):277 """(1/π) Im[Π_N(½+iT) + log((s−1)/s)] — the term the raw count carries beyond the exact count (Prop. combformula);278 Im log((s−1)/s) = π − 2 arctan(2T)."""279 return pole_term(N, T).imag / math.pi + 1.0 - (2.0 / math.pi) * math.atan(2.0 * T)280281_POLE_CACHE = {}282def count_corrected(comb, T):283 """N̂_N(T) = Ñ_N(T) − (1/π) Im[Π_N + log((s−1)/s)]: the pole-corrected comb count (Prop. combformula; Obs. trace).284 For ζ only (const = 1, the archimedean phase θ); a Dirichlet L-function has no pole and no correction."""285 T = np.atleast_1d(np.asarray(T, float))286 corr = np.array([_POLE_CACHE.setdefault((comb.N, float(t)), pole_count_term(comb.N, float(t))) for t in T])287 val = count_raw(comb, T) - corr288 return val if val.size > 1 else float(val[0])289290def secular_roots(comb, windows, const=1.0, theta_fn=theta, corrected=False):291 """The secular roots: solutions of N̂_N(T) = n − ½ (corrected=True) or Ñ_N(T) = n − ½ (raw) in each (n, lo, hi)292 window (Obs. trace)."""293 roots = []294 for n, lo, hi in windows:295 if corrected:296 f = lambda T: count_corrected(comb, T) - (n - 0.5)297 else:298 f = lambda T: count_raw(comb, T, const, theta_fn) - (n - 0.5)299 roots.append(brentq(f, lo, hi, xtol=1e-12))300 return np.array(roots)301302def crossings(count_fn, lo, hi, step=0.1):303 """Brackets of every half-integer crossing of a count on [lo, hi], from the count alone: (level, a, b, direction).304 No reference height enters; the heights are compared only afterwards."""305 ts = np.arange(lo, hi + 1e-9, step); v = np.asarray(count_fn(ts), float)306 out = []307 for i in range(len(ts) - 1):308 for h in np.arange(math.floor(min(v[i], v[i + 1])) + 0.5, max(v[i], v[i + 1]), 1.0):309 if (v[i] - h) * (v[i + 1] - h) < 0:310 out.append((float(h), float(ts[i]), float(ts[i + 1]), int(np.sign(v[i + 1] - v[i]))))311 return out312313def secular_roots_scan(comb, lo, hi, const=1.0, theta_fn=theta, corrected=False, step=0.1):314 """The secular roots on [lo, hi] with brackets from the count itself (Obs. trace): every upward half-integer315 crossing refined by Brent's method. Returns (roots, levels, extra) — extra lists any downward or repeated crossing."""316 f = (lambda T: count_corrected(comb, T)) if corrected else (lambda T: count_raw(comb, T, const, theta_fn))317 cr = crossings(f, lo, hi, step)318 up = [(h, a, b) for h, a, b, d in cr if d > 0]319 levels = [h for h, _, _ in up]320 extra = [(h, a) for h, a, b, d in cr if d < 0] + [(h, a) for h, a, b in up if levels.count(h) > 1]321 roots = np.array([brentq(lambda T: float(np.atleast_1d(f(T))[0]) - h, a, b, xtol=1e-12) for h, a, b in up])322 return roots, levels, extra323324# ----------------------------------------------------------------------------- matrices325def colleague_roots(comb, TA, TB, deg, corrected=False):326 """Eigenvalues of the colleague matrix of the Chebyshev fit of the secular determinant Ξ(T) = cos(π N̂_N(T))327 (corrected=True; Obs. matrix) or cos(π Ñ_N(T)) (raw) on [TA, TB]. Returns the real roots inside the window."""328 from numpy.polynomial import chebyshev as Ch329 xs = np.cos(np.pi * (np.arange(deg + 1) + 0.5) / (deg + 1))330 T = 0.5 * (TB + TA) + 0.5 * (TB - TA) * xs331 cnt = count_corrected(comb, T) if corrected else count_raw(comb, T)332 c = Ch.chebfit(xs, np.cos(math.pi * cnt), deg)333 n = len(c) - 1334 a = c[:-1] / c[-1]335 M = np.zeros((n, n)); M[0, 1] = 1.0336 for k in range(1, n - 1):337 M[k, k - 1] = 0.5; M[k, k + 1] = 0.5338 M[n - 1, n - 2] = 0.5339 M[-1, :] -= 0.5 * a340 ev = np.linalg.eigvals(M)341 ev = ev[np.abs(ev.imag) < 1e-5].real342 ev = ev[(ev > -1) & (ev < 1)]343 return np.sort(0.5 * (TB + TA) + 0.5 * (TB - TA) * ev)344345def jacobi_from_points(x):346 """Real-symmetric tridiagonal Jacobi matrix of the unit-weight measure Σ δ_{x_i} (Stieltjes–Lanczos,347 full reorthogonalisation). Its eigenvalues are exactly the x_i (Def. jacobi)."""348 x = np.asarray(x, float); K = len(x)349 A = np.diag(x); V = np.zeros((K, K))350 v = np.ones(K) / math.sqrt(K); V[:, 0] = v351 alpha = np.zeros(K); beta = np.zeros(max(K - 1, 0))352 r = A @ v; alpha[0] = v @ r; r = r - alpha[0] * v353 for j in range(1, K):354 beta[j - 1] = np.linalg.norm(r)355 vj = r / beta[j - 1]356 for k in range(j):357 vj -= (V[:, k] @ vj) * V[:, k]358 vj /= np.linalg.norm(vj); V[:, j] = vj359 r = A @ vj; alpha[j] = vj @ r360 r = r - alpha[j] * vj - beta[j - 1] * V[:, j - 1]361 return np.diag(alpha) + np.diag(beta, 1) + np.diag(beta, -1)362363# ----------------------------------------------------------------------------- validation markers364# The first thirty nontrivial zero heights of ζ (Odlyzko's tables). Validation markers only:365# no construction in this package takes them as input.366HEIGHTS = np.array([367 14.134725141734693, 21.022039638771555, 25.010857580145688, 30.424876125859513, 32.935061587739189,368 37.586178158825671, 40.918719012147495, 43.327073280914999, 48.005150881167159, 49.773832477672302,369 52.970321477714460, 56.446247697063394, 59.347044002602353, 60.831778524609809, 65.112544048081606,370 67.079810529494173, 69.546401711173979, 72.067157674481907, 75.704690699083933, 77.144840068874805,371 79.337375020249367, 82.910380854086030, 84.735492980517050, 87.425274613125229, 88.809111207634465,372 92.491899270558484, 94.651344040519886, 95.870634228245309, 98.831194218193692, 101.317851005731392])373374def riemann_zeros(n, dps=18):375 """The first n zero heights γ_k of ζ from mpmath's zetazero, cached on disk. Validation markers only."""376 path = os.path.join(FIGDIR, f".zeros_{n}.npy")377 if os.path.exists(path):378 return np.load(path)379 with mp.workdps(dps):380 g = np.array([float(mp.im(mp.zetazero(k))) for k in range(1, n + 1)])381 np.save(path, g)382 return g383384def nearest_errors(roots, heights):385 rec = np.array([roots[np.argmin(np.abs(roots - g))] for g in heights])386 return np.abs(rec - heights)387388class Timer:389 def __init__(self, name): self.name = name390 def __enter__(self): self.t = time.time(); return self391 def __exit__(self, *a): print(f" ({self.name}: {time.time() - self.t:.1f} s)")392393# ======================================================================================================================394# block A — a_shell.py (until 24 September 2026), its former docstring:395# a_shell.py — block A: the shell theorem and the exact shell arithmetic (EXACT, integer-pinned)396# =============================================================================================397# Master ledger predicate 00:D11 (the shell theorem, 20-rh Thm. hp); paper ledger predicates 20:B1–B11, 20:E1, 20:E9, 20:E12–E13. Package checks:398#399# A1 Thm. zeroslot zero-slot: Σ_{x∈F_p^×} x^k = 0 for every nonterminal k, = −1 on the terminal slot400# A2 Thm. compl slot complementarity: Φ(k) = −g^k bijects the nonterminal slots onto F_p^× \ {−1}401# A3 Lem. K Hermitian phase calculus on F_{p²} = F_p(η), η² = ν a non-residue: |U_{p+1}| = p+1,402# Frobenius = inversion on U_{p+1}, i = √−1 Frobenius-fixed and off the circle (N(i) = −1)403# A4 Prop. critical the finite critical line: Tr z = 1 exactly on z = 2⁻¹ + bη; |L_{1/2}| = p;404# N(z) = ¼ − νb²; 2⁻¹ = 2κ+1 = −π (capacity-first)405# A5 Thm. agree the Klein four-group ⟨φ, ρ⟩ and its fixed loci F_p, L_{1/2}, {2⁻¹}; F_p ∩ L_{1/2} = {2⁻¹}406# A6 §1.3 (register) the Subject constants: π = 2κ, 2π ≡ −1, i = g^{−κ}, i² ≡ −1, e = g^i on the odd407# lift, g^π ≡ −1, e^{iπ} ≡ −1 (F_13: g = 2, κ = 3, i = 5, e = 6)408# A6b Rem. frame the odd-lift rule on every frame p ≡ 1 (mod 4) below 3000, least primitive root: with r the409# least residue of i, the lift is r if odd else r + p; on it e^{iπ} ≡ −1, on the even member410# e^{r'π} ≡ +1, so the lift matters in an exponent (F_17: g = 3, i = 4, lift 21, e = 5)411# A7 Thm. hp (iii) EXACT: the scale-shift on the power characters, S x^k = g^k x^k in F_p, and its trace412# Tr S^r = (p−1)·[(p−1) | r] (§10.2: the shell's trace carries no prime data) — integer arithmetic413# A7b Thm. hp (i),(iii) [approx]: the complex reading — the constant mode carries the mean v̄ = ψ(p−1)/(p−1),414# the nontrivial characters carry v − v̄·1 (Rem. parseval), the normalised characters χ_j/√(p−1) are orthonormal and415# S χ_j = ω^j χ_j on ℓ²(F_p^×) — floating-point roots of unity, tolerance 10⁻¹⁰416# A8 Def. jacobi [approx]: self-adjointness is free: any real multiset is the spectrum of a real-symmetric417# tridiagonal matrix (Def. jacobi), checked on the first ten heights in floating point (10⁻¹¹)418# A9 Prop. ground EXACT: flat ground state, the Ramanujan sum c_p(n) = −1 for every n ≢ 0 (mod p), computed in419# F_q with a primitive p-th root of unity (q the least prime ≡ 1 mod p): Σ_{a=1}^{p−1} ω^{an} ≡ −1420# A10 Def. shells the shared structure of two shells is the quarter-turn core Q₄ (4 | p−1, 4 | Ω−1);421# on the laboratory pair (13, 233) the cycle projection C_{Ω−1} → C_{p−1} does not422# exist (12 ∤ 232) — the negative check behind the round-02 chronon paragraph423# A11 Prop. coincide frame coincidence below the horizon √p: residues, window products and primality424# (trial division inside the window) agree between F_p and the Carrier chart F_Ω,425# Ω the least prime ≡ 1 (mod 4) above p²426# A12 Thm. closure the closure of the heights on the laboratory Carriers Ω = 13, 233, 30089, 2408561: the427# largest embedded Subject p* = 3, 13, 173, 1549 (p*² < Ω < p'²), the cells of the cycle as its428# residues, the embedded cycle shorter, the wrap of the height reading, the pigeonhole of429# 00:A2/00:A4 on the cells (rows A2, A4, F7)430# A13 Completeness row A3 (00:A3) on the laboratory pairs (p, Ω): the registration of the cycle is onto its registered431# Universe (every residue of it has a member, whole and under the Subject's reading), the Subject's cycle432# registers within the Carrier's, no reading of the Carrier's cycle within the Subject's is injective, one433# count both ways (FRC.Foundation.registration_onto, mutual_registration_count)434#435# Shells: p = 13, 17, 29, 37, 41 in full (all p² points of F_{p²}); 173 for A1, A4, A6; 1009 and 10009 for A11.436# Kinds: A1–A7, A6b, A9–A13 are integer arithmetic (EXACT); A7b and A8 use floating-point roots of unity and eigenvalues ([approx]).437438SHELLS = [13, 17, 29, 37, 41]439BIG = 173440441def prime_factors(n):442 f, d = set(), 2443 while d * d <= n:444 while n % d == 0:445 f.add(d); n //= d446 d += 1447 if n > 1: f.add(n)448 return f449450def generator(p):451 for g in range(2, p):452 if all(pow(g, (p - 1) // q, p) != 1 for q in prime_factors(p - 1)):453 return g454455def nonresidue(p):456 sq = {(x * x) % p for x in range(1, p)}457 return next(v for v in range(2, p) if v not in sq)458459# F_{p^2} = F_p(η), η² = ν: elements (a, b) = a + bη460def mul(z, w, nu, p):461 return ((z[0] * w[0] + nu * z[1] * w[1]) % p, (z[0] * w[1] + z[1] * w[0]) % p)462def frob(z, p): # z ↦ z^p : a + bη ↦ a − bη463 return (z[0], (-z[1]) % p)464def trace(z, p): return (2 * z[0]) % p465def norm(z, nu, p): return (z[0] * z[0] - nu * z[1] * z[1]) % p466def zpow(z, k, nu, p):467 r, b = (1, 0), z468 while k:469 if k & 1: r = mul(r, b, nu, p)470 b = mul(b, b, nu, p); k >>= 1471 return r472473def block_A():474 """Block A — the shell theorem and the exact shell arithmetic: 62 checks on the shells p = 13, 17, 29, 37, 41 (173, 1009, 10009 and the frames below 3000 for single checks), 56 integer-exact and 6 [approx]."""475 print("\n== block A: the shell theorem (00:D11) ==")476 for p in SHELLS + [BIG]:477 g, nu = generator(p), nonresidue(p)478 kappa = (p - 1) // 4; assert 4 * kappa + 1 == p479 # A1 zero-slot480 ok = all(sum(pow(x, k, p) for x in range(1, p)) % p == 0 for k in range(1, p - 1)) \481 and sum(pow(x, p - 1, p) for x in range(1, p)) % p == p - 1482 # 20:B4 (p20011)483 check("A1", f"zero-slot on F_{p}", ok, f"Σx^k=0 for k=1..{p-2}; = -1 at k={p-1}")484 # A6 the Subject constants, capacity-first485 pi_ = 2 * kappa486 i_ = pow(g, p - 1 - kappa, p) # g^{-κ}487 assert i_ == (-pow(g, kappa, p)) % p488 i_odd = i_ if i_ % 2 == 1 else i_ + p # the odd integer lift of i (Remark frame)489 e_ = pow(g, i_odd, p)490 ok = ((2 * pi_) % p == p - 1 and (i_ * i_) % p == p - 1 and pow(g, pi_, p) == p - 1491 and pow(e_, i_odd * pi_, p) == p - 1 and (2 * (2 * kappa + 1)) % p == 1492 and (2 * kappa + 1) % p == (-pi_) % p)493 # 20:B10 (p20017)494 check("A6", f"Subject constants on F_{p}: π=2κ, i=g^-κ, e=g^i (odd lift)", ok,495 f"g={g} κ={kappa} π={pi_} i={i_} e={e_}: 2π≡-1, i²≡-1, g^π≡-1, e^(iπ)≡-1, 2⁻¹=2κ+1=-π")496 # A4 the finite critical line497 inv2 = (2 * kappa + 1) % p498 line = [(inv2, b) for b in range(p)]499 ok = all(trace(z, p) == 1 for z in line) and len(set(line)) == p \500 and all(norm(z, nu, p) == (pow(4, p - 2, p) - nu * z[1] * z[1]) % p for z in line)501 if p != BIG: # the converse on the full plane: Tr z = 1 only on the line502 ok = ok and all(((2 * a) % p == 1) == (a == inv2) for a in range(p))503 # 20:B8 (p20015)504 check("A4", f"finite critical line on F_{p}: Tr z = 1 ⟺ z = 2⁻¹ + bη, |L| = p, N = ¼ − νb²", ok,505 f"ν={nu}, 2⁻¹={inv2}")506 if p == BIG:507 continue508 # A2 slot complementarity509 image = {(-pow(g, k, p)) % p for k in range(1, p - 1)}510 # 20:B5 (p20012)511 check("A2", f"Φ(k) = −g^k bijects the nonterminal slots onto F_{p}^× \\ {{−1}}", image == set(range(1, p - 1)),512 f"|image| = {len(image)} = p−2")513 # A3 Hermitian phase calculus514 K = [(a, b) for a in range(p) for b in range(p)]515 U = [z for z in K if norm(z, nu, p) == 1]516 inv_ok = all(frob(z, p) == zpow(z, p, nu, p) and mul(z, frob(z, p), nu, p) == (1, 0) for z in U)517 i2 = [(a, 0) for a in range(p) if (a * a) % p == p - 1]518 ok = len(U) == p + 1 and inv_ok and len(i2) == 2 and all(frob(z, p) == z and norm(z, nu, p) == p - 1 for z in i2)519 # 20:B6 (p20013)520 check("A3", f"F_{p}²: |U_(p+1)| = p+1, Frobenius = inversion on it, i Frobenius-fixed with N(i) = −1", ok)521 # A5 Klein four-group522 phi = lambda z: frob(z, p)523 rho = lambda z: ((1 - z[0]) % p, (-z[1]) % p)524 sig = lambda z: rho(phi(z))525 fix_phi = {z for z in K if phi(z) == z}; fix_rho = {z for z in K if rho(z) == z}; fix_sig = {z for z in K if sig(z) == z}526 ok = (all(phi(phi(z)) == z and rho(rho(z)) == z and phi(rho(z)) == rho(phi(z)) for z in K)527 and fix_phi == {(a, 0) for a in range(p)} and fix_sig == set(line) and fix_rho == {(inv2, 0)}528 and fix_phi & fix_sig == {(inv2, 0)})529 # 20:B9 (p20016)530 check("A5", f"Klein four-group on F_{p}²: Fix φ = F_p, Fix σ = L_1/2, Fix ρ = {{2⁻¹}}, F_p ∩ L_1/2 = {{2⁻¹}}", ok)531 # A9 flat ground state, in exact arithmetic: a primitive p-th root of unity ω in F_q, q ≡ 1 (mod p) prime532 q = p + 1533 while not (q % p == 1 and sp.isprime(q)): q += p534 gq = sp.primitive_root(q); om = pow(gq, (q - 1) // p, q)535 assert pow(om, p, q) == 1 and om != 1536 cp = [sum(pow(om, a * n, q) for a in range(1, p)) % q for n in range(1, p)]537 # 20:B7 (p20014)538 check("A9", f"c_{p}(n) ≡ −1 (mod {q}) for n ≠ 0: Σ_a ω^{{an}} in F_{q}, ω of order {p}", all(c == q - 1 for c in cp))539 # A7 spectral content and the two readings of the characters540 lam = {}541 for n in range(1, p):542 m, q = n, None543 for d in range(2, n + 1):544 if n % d == 0:545 q = d; break546 if q is None: lam[n] = 0.0; continue547 while m % q == 0: m //= q548 lam[n] = math.log(q) if m == 1 else 0.0549 v = np.array([lam[n] for n in range(1, p)])550 dlog = {pow(g, m, p): m for m in range(p - 1)}551 w = cmath.exp(2j * math.pi / (p - 1))552 chi = np.array([[w ** (j * dlog[n]) for n in range(1, p)] for j in range(p - 1)])553 coef = chi.conj() @ v / (p - 1)554 mean = v.mean()555 recon_nontriv = (coef[1:] @ chi[1:]).real556 ok_mean = abs(coef[0].real - mean) < 1e-12 and np.allclose(recon_nontriv, v - mean, atol=1e-10) \557 and np.allclose((coef @ chi).real, v, atol=1e-10)558 gram = chi.conj() @ chi.T / (p - 1)559 ok_orth = np.allclose(gram, np.eye(p - 1), atol=1e-10)560 # S f(x) = f(g x): power characters in F_p, complex characters on ℓ²561 ok_fp = all(pow((g * x) % p, k, p) == (pow(g, k, p) * pow(x, k, p)) % p for k in range(p - 1) for x in range(1, p))562 shifted = lambda j: np.array([w ** (j * dlog[(g * n) % p]) for n in range(1, p)]) # (S χ_j)(n) = χ_j(g n)563 ok_c = all(np.allclose(shifted(j), w ** j * chi[j]) for j in range(p - 1))564 tr = [sum(1 for x in range(1, p) if (pow(g, r, p) * x) % p == x) for r in range(1, p)]565 ok_tr = all(tr[r - 1] == ((p - 1) if r % (p - 1) == 0 else 0) for r in range(1, p))566 # 20:E1 (p20032), 20:E12 (p20043)567 check("A7", f"F_{p}: S x^k = g^k x^k on the power characters (integer arithmetic); Tr S^r = (p−1)[(p−1)|r]", ok_fp and ok_tr)568 # 20:E13 (p20044)569 check("A7b", f"ℓ²(F_{p}^×): mean v̄ = ψ(p−1)/(p−1) on the trivial mode, v − v̄·1 on the nontrivial characters, χ_j/√(p−1) orthonormal; Sχ_j = ω^j χ_j",570 ok_mean and ok_orth and ok_c, f"v̄ = {mean:.6f} (= log 27720/12 at p = 13)" if p == 13 else "", kind="[approx]")571 # A10 the shared quarter-turn core; no cycle projection onto (13, 233)572 Om = 233573 q4 = {pow(i_, a, p) for a in range(4)}574 ok = q4 == {1, i_, p - 1, (-i_) % p} and (p - 1) % 4 == 0 and (Om - 1) % 4 == 0575 if p == 13:576 ok = ok and (Om - 1) % (p - 1) != 0577 # 20:B1 (p20008)578 check("A10", "Q₄ ⊂ both cycles; C_232 ↠ C_12 does not exist (12 ∤ 232) on the pair (13, 233)", ok, "the shells share the core, not a clock")579 else:580 check("A10", f"Q₄ = {{1, i, −1, −i}} ⊂ F_{p}^×, 4 | p−1", ok)581 # A11 frame coincidence below the horizon (Prop. coincide): the Carrier chart Ω is the least prime ≡ 1 (mod 4) above p²582583 for p in SHELLS + [1009, 10009]:584 Om = int(sp.nextprime(p * p))585 while Om % 4 != 1:586 Om = int(sp.nextprime(Om))587 H = math.isqrt(p)588 same = all(n % p == n % Om == n for n in range(1, H + 1))589 prods = all((a * b) % p == (a * b) % Om == a * b for a in range(1, H + 1) for b in range(1, H + 1))590 irreducible = {n for n in range(2, H + 1) if all(n % d for d in range(2, math.isqrt(n) + 1))} # trial division inside the window591 primes = {n for n in range(2, H + 1) if sp.isprime(n)}592 # 20:B2 (p20009)593 check("A11", f"frame coincidence below √{p} = {H} on the pair ({p}, {Om}): residues, window products and primality agree", same and prods and irreducible == primes,594 f"Π_p = {sorted(primes)}")595 # A6b the odd-lift rule (Remark frame) on every admissible frame below 3000, least primitive root596 frames = [q for q in range(5, 3000) if q % 4 == 1 and sp.isprime(q)]597 bad = []598 for q in frames:599 gq = int(sp.primitive_root(q)); kq = (q - 1) // 4; piq = 2 * kq600 r = pow(gq, q - 1 - kq, q) # the least residue of i = g^{-κ}601 odd, even = (r, r + q) if r % 2 == 1 else (r + q, r)602 e_odd = pow(gq, odd, q)603 if not (pow(e_odd, odd * piq, q) == q - 1 and pow(e_odd, even * piq, q) == 1 and (r * r) % q == q - 1):604 bad.append(q)605 f17 = (pow(3, 16 - 4, 17), pow(3, 21, 17)) # F_17 with g = 3: i = 3^{-4} = 4 (even), lift 21, e = 3^21 = 5606 ok = not bad and f17 == (4, 5)607 check("A6b", f"the odd-lift rule on all {len(frames)} frames p ≡ 1 (mod 4) below 3000: e^{{iπ}} ≡ −1 on the odd lift, ≡ +1 on the even member; F_17: g=3, i=4, lift 21, e=5", ok,608 f"frames {len(frames)}, failures {bad}; F_17 (i, e) = {f17}")609 # A8 self-adjointness is free610 J = jacobi_from_points(HEIGHTS[:10])611 ev = np.sort(np.linalg.eigvalsh(J))612 ok = np.allclose(J, J.T) and np.max(np.abs(np.triu(J, 2))) == 0 and np.max(np.abs(ev - HEIGHTS[:10])) < 1e-11613 check("A8", "any real multiset is the spectrum of a real-symmetric tridiagonal matrix (Lanczos on ten heights, floating point)", ok,614 f"max|eig(J) − height| = {np.max(np.abs(ev - HEIGHTS[:10])):.1e}", kind="[approx]")615616 # A12 the closure of the heights (Thm. closure; rows A2, A4, F7), on the four laboratory Carriers Ω = 13, 233, 30089,617 # 2408561: the largest embedded Subject p* with p*² < Ω < p'² (p' the next prime); the cells of a cycle are its residues618 # 0..q−1 and a Subject of cardinality p reads the cells y < p, its own residues, so an embedded Subject's cycle is shorter619 # than the Carrier's and the Carrier's top cell is read from labels up to √Ω, beyond every embedded window √p*; the620 # wrap of the height reading, (T + Ω) ≡ T; the pigeonhole of 00:A2/00:A4 on the cells — a reading of the Ω cells into621 # o < Ω records identifies two cells (found by search for o = p* and o = Ω − 1), and an injective reading needs o ≥ Ω622 for Om in (13, 233, 30089, 2408561):623 ps = int(sp.prevprime(math.isqrt(Om) + 1)); pn = int(sp.nextprime(ps))624 ok = sp.isprime(Om) and Om % 4 == 1 and sp.isprime(ps) and sp.isprime(pn) and ps * ps < Om < pn * pn625 ok &= all(not sp.isprime(q) for q in range(ps + 1, pn)) # p* the largest prime below √Ω626 ok &= ps < Om and math.isqrt(ps) < math.isqrt(Om) # the embedded cycle shorter; the top cell's labels beyond √p*627 ok &= all((T + Om) % Om == T % Om for T in (0, 1, ps, Om - 1, Om, 7 * Om + 3)) # the wrap of the height reading628 collisions = {}629 for o in (ps, Om - 1):630 seen = {}631 for y in range(Om): # the cells of the Carrier's cycle, read into o records632 r = y % o633 if r in seen:634 collisions[o] = (seen[r], y); break635 seen[r] = y636 ok &= o in collisions and collisions[o][0] != collisions[o][1] and collisions[o][0] % o == collisions[o][1] % o637 ok &= len({y for y in range(Om)}) == Om # the identity reading into Ω records is injective: Ω ≤ Ω638 # 20:F7 (p20079)639 check("A12", f"closure of the heights on Ω = {Om}: p* = {ps} ({ps}² < Ω < {pn}²), the cycle's {Om} cells, an embedded Subject reads {ps} of them; readings into {ps} and {Om - 1} records collide", ok,640 f"collisions {collisions}; √Ω = {math.isqrt(Om)}, √p* = {math.isqrt(ps)}")641642 # A13 Completeness (row A3; 00:A3), on the laboratory pairs (p, Ω) = (13, 233), (173, 30089), (1549, 2408561): the643 # registration of a totality is onto its registered Universe — the registered Universe is the set of residues at which some644 # member registers, read back as a member for every residue of it (the Carrier's cells registered whole, ρ(y) = y, and the645 # Subject's reading of them, ρ_p(y) = y mod p: Ω and p residues, each with a member); and two totalities each registered646 # whole within the other have one count — the identity registers the Subject's cycle within the Carrier's, every reading of647 # the Carrier's Ω cells into the Subject's p records identifies two cells (y mod p, an affine reading, the square; the648 # pigeonhole of 00:A2 both ways is FRC.Foundation.mutual_registration_count), and only cycles of one count embed both ways649 for p, Om in ((13, 233), (173, 30089), (1549, 2408561)):650 members = {}651 for y in range(Om): # the Carrier's cycle registered whole652 members.setdefault(y % Om, y)653 ok = len(members) == Om and all(members[r] % Om == r for r in members) # onto: every residue has a member654 membersp = {}655 for y in range(Om): # the Subject's reading of the Carrier's cells656 membersp.setdefault(y % p, y)657 ok &= len(membersp) == p and all(membersp[r] % p == r for r in membersp)658 ok &= len({x % Om for x in range(p)}) == p # the identity registers Shell p within Shell Ω659 readings = {"y mod p": lambda y: y % p, "3y+5 mod p": lambda y: (3 * y + 5) % p, "y² mod p": lambda y: (y * y) % p}660 collisions = {}661 for name, f in readings.items(): # no reading of Shell Ω within Shell p is injective662 seen = {}663 for y in range(Om):664 r = f(y)665 if r in seen:666 collisions[name] = (seen[r], y); break667 seen[r] = y668 ok &= name in collisions and collisions[name][0] != collisions[name][1] and f(collisions[name][0]) == f(collisions[name][1])669 ok &= p < Om and len({x % p for x in range(p)}) == p # one count: the identity both ways on Shell p670 # 20:A3 (p00003)671 check("A13", f"completeness on the pair ({p}, {Om}): the registration of the cycle is onto its registered Universe ({Om} residues; {p} under the Subject's reading); F_{p} registers within F_{Om}, every reading of the {Om} cells into {p} records collides, one count only both ways", ok,672 "collisions " + ", ".join(f"{k}: {v}" for k, v in collisions.items()))673674# ======================================================================================================================675# block B — b_deframe.py (until 24 September 2026), its former docstring:676# b_deframe.py — block B: the de-framing dictionary and its consequences ([chart], [approx])677# =========================================================================================678# The dictionary of Definition deframe: a shell of cardinality p reads ζ on the line to height679# T = 2πp with the Riemann–Siegel main sum of length √p = √(T/2π); a height is the scale680# coordinate of the shell that reads it (T/2π = p). Paper-local predicates:681#682# B1 Prop. density [chart] smooth zero density at T = 2πp equals the shell scale-depth (1/2π) log p:683# 0.6226, 0.9891, 1.3556, 1.7220 at p = 50, 500, 5·10³, 5·10⁴ (a substitution684# into the Riemann–von Mangoldt formula; a consistency of the dictionary)685# B2 Obs. residue [approx] the finite-horizon reconstruction error |Z_horizon − Z| scales as p^{-1/4}686# (unconditional: truncation of the main sum), fitted slope ≈ −1/4687# B3 Prop. critical [chart] the de-framing limit of the critical real part: 2⁻¹/p = (2κ+1)/p → ½688# B4 Def. deframe [approx] the horizon-length main sum (length ⌊√(t/2π)⌋, weight n^{-1/2}) tracks Z(t)689# through the first ten zeros: sign changes of the main sum on [10, 55]690# bracket every one of γ_1..γ_10 (the on-line reading of Proposition turing)691#692# The Hardy Z function Z(t) enters here as the classical side of the correspondence (the paper's693# "ζ never evaluated" discipline binds the prime-side realisation of block C, not this dictionary694# check). Figure: fig_deframing.pdf.695696697def Z_true(t):698 return float(mp.siegelz(t))699700def Z_horizon(t):701 """Riemann–Siegel main sum of length ⌊√(t/2π)⌋ = the shell horizon √p at T = 2πp; no remainder."""702 N = int(mp.floor(mp.sqrt(t / (2 * mp.pi))))703 th = mp.siegeltheta(t)704 return float(2 * mp.fsum(mp.cos(th - t * mp.log(n)) / mp.sqrt(n) for n in range(1, N + 1)))705706def block_B():707 """Block B — the de-framing dictionary: the zero density as shell scale-depth, the reconstruction residue, the horizon-length main sum through the first zeros; 4 checks ([chart], [approx])."""708 mp.mp.dps = 18709 print("\n== block B: the de-framing dictionary (Def. deframe) ==")710 # B1 density711 stated = {50: 0.6226, 500: 0.9891, 5000: 1.3556, 50000: 1.7220}712 got = {p: math.log(p) / (2 * math.pi) for p in stated}713 ok = all(abs(got[p] - stated[p]) < 6e-5 for p in stated)714 # illustration: 20:D2 retired 2026-10-15 (T05)715 check("B1", "smooth zero density at T = 2πp is (1/2π) log p", ok,716 " ".join(f"p={p}: {got[p]:.4f}" for p in stated), kind="[chart]")717 # B3 de-framing limit718 vals = [(p, (2 * ((p - 1) // 4) + 1) / p) for p in [13, 53, 101, 409, 1009, 10009]]719 ok = all(v > 0.5 for _, v in vals) and all(vals[i][1] > vals[i + 1][1] for i in range(len(vals) - 1)) and abs(vals[-1][1] - 0.5) < 1e-3720 check("B3", "de-framing limit (2κ+1)/p → ½ from above", ok, " ".join(f"{p}: {v:.5f}" for p, v in vals), kind="[chart]")721 # B4 main sum brackets the first ten zeros722 with Timer("horizon sums on [10,55]"):723 ts = np.linspace(10, 55, 700 if not FAST else 350)724 zt = np.array([Z_true(t) for t in ts]); zh = np.array([Z_horizon(t) for t in ts])725 sc = ts[:-1][zh[:-1] * zh[1:] < 0]726 inwin = HEIGHTS[(HEIGHTS > 10) & (HEIGHTS < 55)] # γ_1..γ_11727 dev = [float(np.min(np.abs(sc - g))) for g in inwin]728 ok = len(sc) == len(inwin) and max(dev) < 0.5729 # illustration: 20:D4 retired 2026-10-15 (T11)730 check("B4", "the bare horizon-length main sum has one sign change per zero on [10,55], each within 0.5 of γ_n", ok,731 f"{len(sc)} sign changes for {len(inwin)} zeros; max deviation {max(dev):.2f} (1–2 terms only at these heights)", kind="[approx]")732 # B2 residue scaling733 with Timer("residue scan"):734 T0 = np.array([40, 80, 160, 320, 640, 1280, 2560, 5120] if not FAST else [40, 80, 160, 320, 640, 1280])735 err = []736 for t0 in T0:737 e = [abs(Z_horizon(float(t0) + 0.123 * i) - Z_true(float(t0) + 0.123 * i)) for i in range(1, 16)]738 err.append(float(np.median(e)))739 slope = np.polyfit(np.log(T0), np.log(err), 1)[0]740 # illustration: 20:D3 retired 2026-10-15 (T05)741 check("B2", "reconstruction residue ∝ p^{-1/4}: fitted log-log slope in [−0.40, −0.10]", -0.40 < slope < -0.10,742 f"slope = {slope:.3f}", kind="[approx]")743 # figure744 fig, ax = plt.subplots(1, 3, figsize=(15, 4.3))745 ax[0].axhline(0, color="#ccc", lw=0.8)746 for g in HEIGHTS[:10]: ax[0].axvline(g, color="#ccc", lw=0.8, ls=":")747 ax[0].plot(ts, zt, color="#1a1a1a", lw=1.5, label=r"$Z(t)$")748 ax[0].plot(ts, zh, color="#c0392b", lw=1.3, ls="--", label=r"horizon sum, $N=\lfloor\sqrt{t/2\pi}\rfloor$")749 ax[0].plot(HEIGHTS[:10], [0] * 10, "o", color="#1f5fbf", ms=5, zorder=5, label=r"zeros $\gamma_n$")750 ax[0].set_xlabel(r"height $t$ ($=2\pi p$ at shell scale $p$)"); ax[0].set_ylabel(r"$Z(t)$"); ax[0].set_ylim(-6, 6)751 ax[0].set_title("B4: de-framing = finite Riemann–Siegel"); ax[0].legend(fontsize=8)752 ps = np.logspace(1.3, 5, 30)753 ax[1].plot(ps, np.log(ps) / (2 * np.pi), color="#1a1a1a", lw=1.8, label=r"$\frac{1}{2\pi}\log\frac{T}{2\pi}$, $T=2\pi p$")754 ax[1].plot(list(stated), [got[p] for p in stated], "o", color="#c0392b", label=r"shell scale-depth $\frac{1}{2\pi}\log p$")755 ax[1].set_xscale("log"); ax[1].set_xlabel("shell cardinality $p$"); ax[1].set_ylabel("zeros per unit height")756 ax[1].set_title("B1: zero density = shell scale-depth [chart]"); ax[1].legend(fontsize=8)757 ax[2].loglog(T0, err, "o-", color="#c0392b", lw=1.4, label=r"median $|Z_{\rm horizon}-Z|$")758 ax[2].loglog(T0, T0 ** -0.25, color="#1a1a1a", ls="--", label=r"$t^{-1/4}=p^{-1/4}$")759 ax[2].set_xlabel(r"height $t$ ($=2\pi p$)"); ax[2].set_ylabel("reconstruction error")760 ax[2].set_title(f"B2: residue $\\propto p^{{-1/4}}$ (slope {slope:.2f})"); ax[2].legend(fontsize=8)761 plt.tight_layout(); plt.savefig(f"{FIGDIR}/fig_deframing.pdf", bbox_inches="tight"); plt.savefig(f"{FIGDIR}/fig_deframing.png", dpi=110, bbox_inches="tight"); plt.close()762 print(f" wrote {FIGDIR}/fig_deframing.pdf")763 # the paper's Figure "the reading of ζ on the line": (a) the main sum of the window through the first ten zeros, (b) the residue764 def _reading(fig, ax):765 a = ax[0]; a.axhline(0, color="0.8", lw=0.6)766 for g in HEIGHTS[:10]: a.axvline(g, color="0.85", lw=0.6, ls=":")767 a.plot(ts, zt, color=INK, lw=1.2, label=r"$Z(t)$")768 a.plot(ts, zh, color=RED, lw=1.1, ls="--", label=r"main sum to $\lfloor\sqrt{t/2\pi}\rfloor$")769 a.plot(HEIGHTS[:10], [0] * 10, "o", color=BLUE, ms=3.6, zorder=5, label=r"zeros $\gamma_n$")770 a.set_xlim(10, 55); a.set_ylim(-6.3, 7.6); a.set_yticks([-6, -4, -2, 0, 2, 4, 6]); a.set_xlabel(r"height $t$"); a.set_ylabel(r"$Z(t)$")771 a.legend(loc="upper center", ncol=3, columnspacing=1.0, handlelength=1.6); panel_letter(a, "(a)")772 b = ax[1]; b.loglog(T0, err, "o-", color=RED, lw=1.1, ms=3.6, label=r"median $|Z_M-Z|$")773 b.loglog(T0, T0 ** -0.25, color=INK, ls="--", lw=1.0, label=r"$t^{-1/4}$")774 b.set_xlabel(r"height $t$"); b.set_ylabel("residue"); b.legend(loc="lower left"); panel_letter(b, "(b)")775 paper_figure("fig_reading", _reading, 2.45)776777# ======================================================================================================================778# block C — c_primeside.py (until 24 September 2026), its former docstring:779# c_primeside.py — block C: the spectrum from the prime side, ζ never evaluated ([approx])780# =======================================================================================781# Every construction here uses only {Λ(n)} (sieved), their logarithms, and the archimedean782# phase θ(T) from log Γ. The Riemann heights enter as validation markers only.783# Paper-local predicates:784#785# C1 Obs. primespec the scale-spectrum |Σ_N(γ)| of the tapered comb, N = 10⁶, peaks at the first six786# heights (paper: 14.14, 21.02, 25.02, 30.40, 32.96, 37.57; the peak is the taper's, of width ~2π/log N)787# (Root search: the brackets of every root are the half-integer crossings of the count itself on [10, 52], found by a788# 0.1 scan and refined by Brent's method; the reference heights enter only in the errors reported.)789# C2 Obs. trace the raw secular condition Ñ_N(T) = n − ½ recovers the first ten heights: mean error790# 4.4×10⁻⁵, maximum 1.5×10⁻⁴ (at γ₁) at N = 10⁶; 3.6×10⁻⁴ at 10³; 6×10⁻⁵ at 10⁵791# C2b Obs. trace the pole-corrected secular condition N̂_N(T) = n − ½ (Prop. combformula): mean error 3.4×10⁻⁵792# (max 6.6×10⁻⁵) at 10⁶, 2.3×10⁻⁵ at 10⁷, 1.3×10⁻⁵ (max 2.8×10⁻⁵) at 10⁸ — falling with depth793# C2c Obs. trace the raw condition does not sharpen past 10⁷: mean 2.7×10⁻⁵ at 10⁷ and 4.2×10⁻⁵ at 10⁸, the γ₁794# error growing from 1.5×10⁻⁴ (10⁶) to 2.6×10⁻⁴ (10⁸) — the pole term, largest at the lowest height795# C3 Obs. matrix the colleague matrix of the Chebyshev fit of Ξ(T) = cos(π N̂_N(T)) on [10, 52] at N = 10⁶:796# mean error 4.8×10⁻⁵ at dimension 520, 3.5×10⁻⁵ by 620, converging to the corrected roots' own 3.4×10⁻⁵797# C4 Obs. jacobi the comb-built Jacobi matrix of the corrected secular roots: real-symmetric tridiagonal, eigenvalues =798# the roots to 10⁻¹², the heights to 3.4×10⁻⁵; the ten recurrence coefficients a, b as stated799# C5 Prop. gauge the additive injection −i d/du + V_comb is gauge-trivial: spacing standard800# deviation 0.08 (absolute; 0.11 of the mean spacing 0.72) against 0.39 of the mean801# for the heights (package finding: the paper's two figures use two conventions)802#803# Figures: fig_prime_spectrum.pdf (C1), fig_operator_spectrum.pdf (C2, C3, C5).804805# the paper's stated numbers806PEAKS_STATED = [14.14, 21.02, 25.02, 30.40, 32.96, 37.57]807A_STATED = [34.314, 31.021, 31.179, 31.763, 32.812, 36.227, 35.012, 37.157, 36.092, 37.564]808B_STATED = [11.178, 10.431, 9.579, 8.263, 7.008, 8.415, 7.048, 4.093, 4.110]809810LO, HI = 10.0, 52.0 # the shell window of Obs. trace/matrix; the ten roots are the crossings of ½ … 9½ found on it811812def roots_on_window(comb, corrected):813 """Ten secular roots from the count alone (brackets by a 0.1 scan of the window, Brent refinement); the reference814 heights enter only in the error reported afterwards. Fails loudly if the crossings are not exactly ½ … 9½, once each."""815 roots, levels, extra = secular_roots_scan(comb, LO, HI, corrected=corrected)816 assert levels == [k + 0.5 for k in range(10)] and not extra, f"crossings {levels}, extra {extra}"817 return roots818819def additive_spectrum(U=15.0, M=2048, potential_scale=5.0):820 """H = −i d/du + V_comb on a periodic log-grid of length U (Prop. gauge)."""821 du = U / M822 k = np.fft.fftfreq(M, d=du) * 2 * np.pi823 F = np.fft.fft(np.eye(M), axis=0) / np.sqrt(M)824 D = F.conj().T @ np.diag(k) @ F825 D = 0.5 * (D + D.conj().T)826 Np = int(math.e ** U) + 1827 n, lam = von_mangoldt_support(Np)828 V = np.zeros(M)829 for nn, ll in zip(n, lam):830 j = int(round(math.log(nn) / du)) % M831 V[j] += ll / math.sqrt(nn) / du832 H = D + np.diag(potential_scale * V)833 H = 0.5 * (H + H.conj().T)834 return np.sort(np.linalg.eigvalsh(H))835836def _C_primeside():837 print("\n== block C: the prime-side realisation, ζ never evaluated ==")838 N = 10 ** 6839 with Timer(f"comb to {N}"):840 comb = Comb(N)841 # C1 peaks of the scale-spectrum842 with Timer("scale-spectrum"):843 gs = np.arange(10.0, 40.0, 0.005)844 sig = comb.sigma(gs)845 peaks = [gs[i] for i in range(1, len(gs) - 1) if sig[i] > sig[i - 1] and sig[i] > sig[i + 1]]846 top = sorted(sorted(peaks, key=lambda g: -sig[np.argmin(np.abs(gs - g))])[:6])847 err1 = [abs(t - h) for t, h in zip(top, HEIGHTS[:6])]848 err_stated = [abs(t - h) for t, h in zip(top, PEAKS_STATED)]849 base = np.median(sig)850 ratio = min(sig[np.argmin(np.abs(gs - t))] for t in top) / base851 # illustration: 20:E5 retired 2026-10-15 (T05)852 check("C1", "scale-spectrum peaks at the stated 14.14 … 37.57 (within 0.015) and at γ_1..γ_6 (within 0.05, the taper's width); each peak ≥ 4× the baseline",853 max(err_stated) < 0.015 and max(err1) < 0.05 and ratio > 4,854 "peaks " + ", ".join(f"{t:.2f}" for t in top) + f"; max |peak − γ| {max(err1):.3f}; min peak/baseline {ratio:.1f}", kind="[approx]")855 # C2 secular roots at three depths856 errs = {}857 for depth in ([10 ** 3, 10 ** 5, 10 ** 6] if not FAST else [10 ** 5, 10 ** 6]):858 c = comb if depth == N else Comb(depth)859 roots = roots_on_window(c, corrected=False)860 e = np.abs(roots - HEIGHTS[:10]); errs[depth] = (e.mean(), e.max())861 e6 = errs[10 ** 6]862 ok = abs(e6[0] - 4.4e-5) < 0.6e-5 and abs(e6[1] - 1.46e-4) < 0.15e-4 and (FAST or errs[10 ** 3][0] < 6e-4) and errs[10 ** 5][0] < 1e-4863 # illustration: 20:E6 retired 2026-10-15 (T11)864 check("C2", "raw secular condition: ten heights to mean 4.4e-5, max 1.5e-4 (γ₁) at N = 10⁶; 3.6e-4 at 10³, 6e-5 at 10⁵", ok,865 "; ".join(f"N={d:.0e}: mean {m:.2e} max {x:.2e}" for d, (m, x) in errs.items()), kind="[approx]")866 # C2b, C2c: the pole-corrected secular condition against the raw one, with depth (Prop. combformula)867 deep = [10 ** 6, 10 ** 7] + ([10 ** 8] if not FAST else [])868 raw_e, cor_e = {}, {}869 with Timer(f"secular roots, raw and corrected, to {deep[-1]:.0e}"):870 for depth in deep:871 c = comb if depth == N else Comb(depth)872 er = np.abs(roots_on_window(c, corrected=False) - HEIGHTS[:10]); ec = np.abs(roots_on_window(c, corrected=True) - HEIGHTS[:10])873 raw_e[depth] = (er.mean(), er.max(), er[0]); cor_e[depth] = (ec.mean(), ec.max(), ec[0])874 if depth != N: del c875 c6 = cor_e[10 ** 6]876 ok = (abs(c6[0] - 3.4e-5) < 0.4e-5 and abs(c6[1] - 6.6e-5) < 0.8e-5 and abs(cor_e[10 ** 7][0] - 2.3e-5) < 0.4e-5877 and (FAST or (abs(cor_e[10 ** 8][0] - 1.3e-5) < 0.3e-5 and cor_e[10 ** 8][1] < 3.2e-5))878 and all(cor_e[deep[i + 1]][0] < cor_e[deep[i]][0] for i in range(len(deep) - 1)))879 check("C2b", "pole-corrected secular condition: mean 3.4e-5 (max 6.6e-5) at 10⁶, 2.3e-5 at 10⁷, 1.3e-5 (max 2.8e-5) at 10⁸; falling with depth", ok,880 "; ".join(f"N={d:.0e}: mean {m:.2e} max {x:.2e} γ₁ {g:.2e}" for d, (m, x, g) in cor_e.items()), kind="[approx]")881 ok = (abs(raw_e[10 ** 7][0] - 2.7e-5) < 0.4e-5 and (FAST or (abs(raw_e[10 ** 8][0] - 4.2e-5) < 0.5e-5 and raw_e[10 ** 8][0] > raw_e[10 ** 7][0]882 and abs(raw_e[10 ** 8][2] - 2.6e-4) < 0.3e-4 and raw_e[10 ** 8][2] > raw_e[10 ** 6][2])))883 check("C2c", "raw secular condition does not sharpen past 10⁷: mean 2.7e-5 at 10⁷, 4.2e-5 at 10⁸; γ₁ error 1.5e-4 (10⁶) → 2.6e-4 (10⁸): the pole term", ok,884 "; ".join(f"N={d:.0e}: mean {m:.2e} max {x:.2e} γ₁ {g:.2e}" for d, (m, x, g) in raw_e.items()), kind="[approx]")885 roots6 = roots_on_window(comb, corrected=True)886 # C4 Jacobi (of the corrected secular roots)887 J = jacobi_from_points(roots6)888 ev = np.sort(np.linalg.eigvalsh(J))889 a, b = np.diag(J), np.diag(J, 1)890 ok = (np.allclose(J, J.T) and np.max(np.abs(np.triu(J, 2))) == 0 and np.max(np.abs(ev - roots6)) < 1e-12891 and np.max(np.abs(a - A_STATED)) < 2e-3 and np.max(np.abs(b - B_STATED)) < 2e-3892 and abs(np.mean(np.abs(ev - HEIGHTS[:10])) - 3.4e-5) < 0.4e-5)893 # illustration: 20:E9 retired 2026-10-15 (T05)894 check("C4", "Jacobi matrix of the corrected roots: real-symmetric tridiagonal; eig = roots to 1e-12; a, b as stated; heights to 3.4e-5", ok,895 f"max|eig−roots| = {np.max(np.abs(ev - roots6)):.1e}; a[0..2] = {a[0]:.3f}, {a[1]:.3f}, {a[2]:.3f}; b[0] = {b[0]:.3f}", kind="[approx]")896 # C3 colleague matrix of the corrected secular determinant897 with Timer("colleague matrices"):898 degs = [180, 260, 340, 420, 520, 620] if not FAST else [340, 520]899 conv = {d: float(np.mean(nearest_errors(colleague_roots(comb, 10.0, 52.0, d, corrected=True), HEIGHTS[:10]))) for d in degs}900 ok = abs(conv[520] - 4.8e-5) < 0.6e-5 and (FAST or abs(conv[620] - 3.5e-5) < 0.5e-5)901 # illustration: 20:E8 retired 2026-10-15 (T05)902 check("C3", "colleague matrix of cos(π N̂_N) on [10,52], N = 10⁶: mean error 4.8e-5 at dim 520, 3.5e-5 by 620 (the roots' own 3.4e-5)", ok,903 "; ".join(f"dim {d}: {v:.2e}" for d, v in conv.items()), kind="[approx]")904 eig520 = colleague_roots(comb, 10.0, 52.0, 520, corrected=True)905 # C5 additive injection906 with Timer("additive operator"):907 ev_add = additive_spectrum()908 eva = ev_add[(ev_add > 10) & (ev_add < 50)]909 da, dh = np.diff(eva), np.diff(HEIGHTS[:10])910 cv_add, cv_h = np.std(da) / np.mean(da), np.std(dh) / np.mean(dh)911 # illustration: 20:E10 retired 2026-10-15 (T05)912 check("C5", "additive −i d/du + V_comb has a uniform spectrum: spacing std 0.08 (0.11 of its mean) against 0.39 of the mean for the heights",913 abs(np.std(da) - 0.08) < 0.01 and cv_add < 0.15 and abs(cv_h - 0.39) < 0.02,914 f"additive: std {np.std(da):.3f} on mean spacing {np.mean(da):.3f} (ratio {cv_add:.3f}); heights: std {np.std(dh):.2f} on mean {np.mean(dh):.2f} (ratio {cv_h:.3f}) — the paper's 0.08 is the absolute std, its 0.39 the ratio", kind="[approx]")915 # ---- figures916 fig, ax = plt.subplots(1, 2, figsize=(12, 3.8))917 ax[0].plot(gs, sig, color="#10325f", lw=1.0)918 for h in HEIGHTS[:6]: ax[0].axvline(h, color="#c0392b", ls="--", lw=0.8)919 ax[0].plot(top, [sig[np.argmin(np.abs(gs - t))] for t in top], "o", color="#e67e22")920 ax[0].set_xlabel(r"$\gamma$"); ax[0].set_ylabel(r"$|\Sigma_N(\gamma)|$"); ax[0].set_title(f"C1: scale-spectrum of the comb, N = 10⁶ (peaks vs heights)")921 g2 = np.arange(13.0, 15.3, 0.002); sharp = {}922 for Ndepth, col in [(10 ** 4, "#9ec5e8"), (10 ** 5, "#5b9bd5"), (10 ** 6, "#10325f")]:923 c = comb if Ndepth == N else Comb(Ndepth)924 sharp[Ndepth] = (col, c.sigma(g2) / c.sigma(HEIGHTS[0]))925 ax[1].plot(g2, sharp[Ndepth][1], color=col, lw=1.1, label=f"N = 10^{int(math.log10(Ndepth))}")926 ax[1].axvline(HEIGHTS[0], color="#c0392b", ls="--", lw=0.8); ax[1].set_xlabel(r"$\gamma$"); ax[1].legend(fontsize=8)927 ax[1].set_title(r"peak at $\gamma_1$ sharpens with depth, width $\sim 2\pi/\log N$")928 plt.tight_layout(); plt.savefig(f"{FIGDIR}/fig_prime_spectrum.pdf", bbox_inches="tight"); plt.savefig(f"{FIGDIR}/fig_prime_spectrum.png", dpi=110, bbox_inches="tight"); plt.close()929 print(f" wrote {FIGDIR}/fig_prime_spectrum.pdf")930 # the paper's Figure "the zero heights from the prime side": (a) the scale-spectrum with its peaks, (b) the peak at γ₁ sharpening931 def _primespec(fig, ax):932 a = ax[0]933 for h in HEIGHTS[:6]: a.axvline(h, color=RED, ls="--", lw=0.7)934 a.plot(gs, sig, color=NAVY, lw=0.9, label=r"$|\Sigma_N(\gamma)|$, $N=10^{6}$")935 a.plot(top, [sig[np.argmin(np.abs(gs - t))] for t in top], "o", color=ORANGE, ms=4, label="peaks, from the primes alone")936 a.plot([], [], color=RED, ls="--", lw=0.7, label=r"zero heights $\gamma_n$")937 a.set_xlim(10, 40); a.set_ylim(0, 10.2); a.set_xlabel(r"scale frequency $\gamma$"); a.set_ylabel("amplitude"); a.legend(loc="upper right"); panel_letter(a, "(a)")938 b = ax[1]939 for Ndepth, (col, curve) in sharp.items():940 b.plot(g2, curve, color=col, lw=1.1, label=rf"$N=10^{{{int(math.log10(Ndepth))}}}$")941 b.axvline(HEIGHTS[0], color=RED, ls="--", lw=0.7); b.set_xlim(13, 15.3); b.set_ylim(0, 1.32); b.set_xlabel(r"$\gamma$"); b.set_ylabel(r"amplitude at $\gamma_1$, normalised")942 b.legend(loc="upper left"); panel_letter(b, "(b)")943 paper_figure("fig_primespec", _primespec, 2.45)944 Tg = np.arange(11.0, 40.0, 0.02)945 NB = theta(Tg) / math.pi + 1.0; NC = count_raw(comb, Tg); NH = count_corrected(comb, Tg)946 fig, ax = plt.subplots(1, 3, figsize=(13.5, 4.2))947 ax[0].plot(Tg, NC, color="#10325f", lw=1.4, label=r"raw $\tilde{N}_N(T)=\theta/\pi+1+S_{\rm comb}$")948 ax[0].plot(Tg, NH, color="#1f6b4a", lw=1.0, ls="--", label=r"pole-corrected $\hat{N}_N(T)$")949 ax[0].plot(Tg, NB, color="#9ec5e8", lw=1.3, label=r"$\theta(T)/\pi+1$")950 for n in range(2, 9): ax[0].axhline(n - 0.5, color="#bbb", lw=0.6, ls=":")951 for g in HEIGHTS[:6]: ax[0].axvline(g, color="#c0392b", lw=0.8, ls="--", alpha=0.7)952 ax[0].set_xlim(11, 40); ax[0].set_ylim(1, 9); ax[0].set_xlabel("$T$"); ax[0].set_ylabel("raw comb count")953 ax[0].set_title(r"C2/C2b: secular condition, count $= n-\frac{1}{2}$ at the heights", fontsize=10); ax[0].legend(fontsize=8, loc="upper left")954 ax[1].vlines(HEIGHTS[:10], 0, 1, color="#c0392b", ls="--", lw=1.1, label="heights (markers)")955 ax[1].plot(eig520, 0.5 * np.ones_like(eig520), "o", color="#10325f", ms=7, label="colleague eigenvalues")956 ax[1].set_xlim(10, 52); ax[1].set_ylim(0, 1); ax[1].set_yticks([]); ax[1].set_xlabel(r"$\gamma$")957 ax[1].set_title(f"C3: colleague matrix (corrected), dim 520, mean err {conv[520]:.1e}", fontsize=10); ax[1].legend(fontsize=8, loc="lower center")958 axin = ax[1].inset_axes([0.10, 0.58, 0.42, 0.32]); axin.semilogy(list(conv), list(conv.values()), "o-", color="#1f6b4a", ms=3.5)959 axin.set_title("mean error vs dimension", fontsize=7); axin.tick_params(labelsize=6)960 ax[2].vlines(eva, 0.55, 0.95, color="#1f6b4a", lw=1.0); ax[2].vlines(HEIGHTS[(HEIGHTS > 10) & (HEIGHTS < 50)], 0.05, 0.45, color="#c0392b", lw=1.4)961 ax[2].set_xlim(10, 50); ax[2].set_ylim(0, 1); ax[2].set_yticks([0.25, 0.75]); ax[2].set_yticklabels(["heights", "dilation\n+ comb"], fontsize=8.5)962 ax[2].set_xlabel(r"$\gamma$"); ax[2].set_title("C5: additive injection is gauge-trivial", fontsize=10)963 plt.tight_layout(); plt.savefig(f"{FIGDIR}/fig_operator_spectrum.pdf", bbox_inches="tight"); plt.savefig(f"{FIGDIR}/fig_operator_spectrum.png", dpi=110, bbox_inches="tight"); plt.close()964 print(f" wrote {FIGDIR}/fig_operator_spectrum.pdf")965 return comb966967# ======================================================================================================================968# block C, continued — c_gue.py (until 24 September 2026), its former docstring:969# c_gue.py — block C, continued: random-matrix statistics of the target spectrum ([approx])970# =========================================================================================971# The zero heights read as the target spectrum of the scale-evolution operator (Definition972# specmap: ρ ↦ −i(ρ − ½)). The heights here are computed independently with mpmath's zetazero973# and serve as the validation markers of Numerical Observation gue; no construction of the974# package takes them as input. Paper-local predicates:975#976# C6 Obs. gue [approx] unfolded spacings of the first 240 zeros show GUE level repulsion, not977# Poisson: P(s < ½) = 0.05 (GUE ≈ 0.12, Poisson ≈ 0.39), variance 0.13978# (GUE ≈ 0.18, Poisson 1)979# C6b Fig. hp [chart] the eigenvalue count N(T) lies on the Berry–Keating semiclassical980# (T/2π)(log(T/2π) − 1) + 7/8: density (1/2π) log p at T = 2πp, the shell981# scale-depth (Proposition density, block B)982#983# Figure: fig_hilbert_polya.pdf.984985def unfold(g):986 """Smooth counting function N̄(γ) = (γ/2π) log(γ/2πe) + 7/8: unfolded ordinates with unit mean spacing."""987 return g / (2 * math.pi) * np.log(g / (2 * math.pi * math.e)) + 7.0 / 8.0988989def gue_wigner(s):990 return 32 / math.pi ** 2 * s ** 2 * np.exp(-4 * s ** 2 / math.pi)991992def _C_gue():993 print("\n== block C (continued): GUE level repulsion of the target spectrum (Obs. gue) ==")994 n = 240 if not FAST else 120995 with Timer(f"first {n} zeros (mpmath zetazero)"):996 g = riemann_zeros(n)997 u = unfold(g)998 s = np.diff(u); s = s / s.mean()999 p_half = float(np.mean(s < 0.5)); var = float(np.var(s))1000 # reference values of the two ensembles (Wigner surmise for GUE; exponential for Poisson)1001 from scipy.integrate import quad1002 gue_half = quad(gue_wigner, 0, 0.5)[0]; gue_var = 3 * math.pi / 8 - 11003 poi_half = 1 - math.exp(-0.5); poi_var = 1.01004 tol = 0.02 if not FAST else 0.041005 ok = abs(p_half - 0.05) < tol and abs(var - 0.13) < 0.03 and p_half < 0.5 * poi_half and var < 0.5 * poi_var1006 # illustration: 20:E11 retired 2026-10-15 (T05)1007 check("C6", "unfolded spacings of the first 240 zeros: P(s<½) ≈ 0.05, variance ≈ 0.13 (repulsion; far from Poisson)", ok,1008 f"P(s<½) = {p_half:.3f} (GUE {gue_half:.3f}, Poisson {poi_half:.3f}); var = {var:.3f} (GUE {gue_var:.3f}, Poisson 1)", kind="[approx]")1009 # C6b: the count against Berry–Keating1010 Ts = np.linspace(10, g[-1], 200)1011 Ncount = np.searchsorted(g, Ts)1012 Nsemi = Ts / (2 * math.pi) * (np.log(Ts / (2 * math.pi)) - 1) + 7.0 / 8.01013 dev = np.abs(Ncount - Nsemi)1014 # illustration: 20:E3 retired 2026-10-15 (T05)1015 check("C6b", "eigenvalue count N(T) within 1.5 of the Berry–Keating smooth count on [10, γ_240]", dev.max() < 1.5,1016 f"max |N − N̄| = {dev.max():.2f}; mean = {dev.mean():.2f}", kind="[chart]")1017 # figure1018 fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))1019 A = ax[0]1020 A.hist(s, bins=np.linspace(0, 3, 16), density=True, color="#1f5fbf", alpha=0.55, edgecolor="white", label=f"ζ zeros (first {n})")1021 xx = np.linspace(0, 3, 400)1022 A.plot(xx, gue_wigner(xx), color="#16a085", lw=2.2, label="GUE (Wigner surmise)")1023 A.plot(xx, np.exp(-xx), color="#1a1a1a", lw=1.4, ls="--", label="Poisson")1024 A.set_xlabel("normalised spacing $s$"); A.set_ylabel("$P(s)$"); A.set_xlim(0, 3)1025 A.set_title(f"C6: level repulsion — $P(s<\\frac{{1}}{{2}})={p_half:.2f}$, var $={var:.2f}$", fontsize=10); A.legend(fontsize=8.5)1026 B = ax[1]1027 B.step(np.concatenate([[10], g]), np.arange(n + 1), where="post", color="#1f5fbf", lw=1.4, label=r"$N(T)=\#\{\gamma_k\leq T\}$")1028 B.plot(Ts, Nsemi, color="#c0392b", lw=2.0, ls="--", label=r"Berry–Keating $\frac{T}{2\pi}(\log\frac{T}{2\pi}-1)+\frac{7}{8}$")1029 B.set_xlabel(r"$T$ ($=2\pi p$ at shell scale $p$)"); B.set_ylabel("$N(T)$"); B.set_xlim(10, g[-1])1030 B.set_title(r"C6b: density $\frac{1}{2\pi}\log\frac{T}{2\pi}=\frac{1}{2\pi}\log p$, the shell scale-depth", fontsize=10)1031 B.legend(fontsize=8.5, loc="upper left")1032 plt.tight_layout(); plt.savefig(f"{FIGDIR}/fig_hilbert_polya.pdf", bbox_inches="tight"); plt.savefig(f"{FIGDIR}/fig_hilbert_polya.png", dpi=110, bbox_inches="tight"); plt.close()1033 print(f" wrote {FIGDIR}/fig_hilbert_polya.pdf")10341035# ======================================================================================================================1036# block C, continued — c_chi.py (until 24 September 2026), its former docstring:1037# c_chi.py — block C, continued: the χ-twisted comb and the generalized hypothesis ([approx])1038# ==========================================================================================1039# Definition chi states that the construction goes through verbatim for a real (quadratic)1040# character; Numerical Observation chi supports it for χ_{−4} (conductor 4, odd, root number 1).1041# Two arms:1042#1043# * the construction arm — the twisted secular condition N_χ(T) = θ_χ(T)/π + S^χ_comb(T), with1044# θ_χ(t) = Im log Γ(¾ + it/2) + (t/2) log(4/π) and the χ-weighted tapered comb Λ_χ(n) = χ(n)Λ(n),1045# with no pole term; L never evaluated;1046# * the validation arm — L(s, χ_{−4}) = 4^{−s}[ζ(s, ¼) − ζ(s, ¾)] via Hurwitz zeta (mpmath), its1047# completed form real on the line, its zeros located as sign changes.1048#1049# Paper-local predicates:1050#1051# C7a Obs. chi [approx] N_χ evaluated at the first six zeros of L(s, χ_{−4}) (validation arm) reads1052# 0.4997, 1.4998, 2.4999, 3.4999, 4.4999, 5.49981053# C7b Obs. chi [approx] solving the twisted secular condition with the N = 10⁶ comb recovers the six1054# heights 6.0209, 10.2438, 12.9881, 16.3426, 18.2920, 21.4506 to mean error1055# 8.6×10⁻⁵ (the floor of the untwisted case, Obs. trace)1056# C7c Obs. chi [approx] the complex character χ mod 5 with χ(2) = i (odd, conductor 5): the half-phase1057# W(χ)^{−1/2} Λ(½ + it, χ) is real on the line (validation arm, Hurwitz zeta); the1058# twisted count (1/π)[θ_χ(T) − θ_χ(0)] + (1/π) Im[Σ^χ_w(½+iT) − Σ^χ_w(½)], with1059# θ_χ(t) = Im log Γ(¾ + it/2) + (t/2) log(5/π) and no pole term, reads the exact1060# integers 1 … 14 midway between the fifteen zeros below 40; its half-integer1061# crossings, bracketed by the count alone, recover the fifteen heights with the1062# N = 10⁶ comb to mean error 8.0×10⁻⁵ (max 2.0×10⁻⁴), L never evaluated (ledger E16)10631064HEIGHTS_CHI_STATED = [6.0209, 10.2438, 12.9881, 16.3426, 18.2920, 21.4506]1065NCHI_STATED = [0.4997, 1.4998, 2.4999, 3.4999, 4.4999, 5.4998]10661067def chi4(n):1068 """χ_{−4}(n): 0 on even n, +1 on n ≡ 1, −1 on n ≡ 3 (mod 4)."""1069 n = np.asarray(n)1070 return np.where(n % 2 == 0, 0.0, np.where(n % 4 == 1, 1.0, -1.0))10711072def twisted_support(N):1073 """(n, χ(n)Λ(n)) on the prime powers ≤ N with χ(n) ≠ 0."""1074 n, lam = von_mangoldt_support(N)1075 c = chi4(n.astype(np.int64))1076 keep = c != 01077 return n[keep], (c * lam)[keep]10781079def theta_chi(T):1080 T = np.asarray(T, float)1081 return np.imag(loggamma(0.75 + 0.5j * T)) + 0.5 * T * math.log(4 / math.pi)10821083def Lambda_chi_line(t):1084 """Completed Λ(½ + it, χ_{−4}) = (4/π)^{(s+1)/2} Γ((s+1)/2) L(s, χ); real for this character (root number 1)."""1085 s = mp.mpf(1) / 2 + mp.mpc(0, 1) * t1086 L = mp.power(4, -s) * (mp.zeta(s, mp.mpf(1) / 4) - mp.zeta(s, mp.mpf(3) / 4))1087 return mp.power(4 / mp.pi, (s + 1) / 2) * mp.gamma((s + 1) / 2) * L10881089def _C_chi():1090 mp.mp.dps = 201091 print("\n== block C (continued): the χ-twisted comb, χ = χ_{−4} (Def. chi, Obs. chi) ==")1092 # validation arm: the completed L is real on the line; its first six zeros by sign change + refinement1093 with Timer("Hurwitz-zeta validation arm"):1094 ts = np.arange(1.0, 23.0, 0.01)1095 vals = [Lambda_chi_line(float(t)) for t in ts]1096 im_max = max(abs(mp.im(v)) for v in vals)1097 re = np.array([float(mp.re(v)) for v in vals])1098 idx = np.nonzero(re[:-1] * re[1:] < 0)[0]1099 zeros = np.array([brentq(lambda t: float(mp.re(Lambda_chi_line(t))), ts[i], ts[i + 1], xtol=1e-10) for i in idx])[:6]1100 ok_real = im_max < 1e-151101 # construction arm1102 N = 10 ** 61103 with Timer(f"twisted comb to {N}"):1104 comb = Comb(N, lam_fn=twisted_support)1105 nchi = count_raw(comb, zeros, const=0.0, theta_fn=theta_chi)1106 ok = ok_real and len(zeros) == 6 and np.max(np.abs(zeros - HEIGHTS_CHI_STATED)) < 6e-5 and np.max(np.abs(nchi - NCHI_STATED)) < 6e-51107 # illustration: 20:E14, 20:E15 retired 2026-10-15 (T05)1108 check("C7a", "N_χ at the first six zeros of L(s, χ_{−4}) reads 0.4997 … 5.4998 (completed L real on the line)", ok,1109 "zeros " + ", ".join(f"{z:.4f}" for z in zeros) + "; N_χ " + ", ".join(f"{v:.4f}" for v in nchi) + f"; max|Im Λ| = {float(im_max):.1e}", kind="[approx]")1110 # construction arm: the brackets come from the twisted count itself (its half-integer crossings on [3, 23]); the1111 # reference zeros of the validation arm enter only in the error reported1112 roots, levels, extra = secular_roots_scan(comb, 3.0, 23.0, const=0.0, theta_fn=theta_chi)1113 assert levels == [k + 0.5 for k in range(6)] and not extra, f"crossings {levels}, extra {extra}"1114 err = np.abs(roots - zeros)1115 check("C7b", "twisted secular condition, N = 10⁶ comb, L never evaluated: six heights to mean error 8.6e-5", abs(err.mean() - 8.6e-5) < 1.5e-5,1116 "roots " + ", ".join(f"{r:.4f}" for r in roots) + f"; mean error {err.mean():.2e}, max {err.max():.2e}", kind="[approx]")1117 _C_chi_complex()11181119# ----------------------------------------------------------------------------- C7c: the complex character mod 51120CHI5 = np.array([0, 1, 1j, -1j, -1], dtype=complex) # χ(n) by n mod 5: χ(1)=1, χ(2)=i, χ(3)=−i, χ(4)=−11121MP_CHI5 = {1: mp.mpc(1), 2: mp.mpc(0, 1), 3: mp.mpc(0, -1), 4: mp.mpc(-1)}11221123class ComplexTwistedComb:1124 """The χ-weighted tapered comb for a complex character: Σ^χ_w(½ + iT) = Σ_{n≤N} χ(n)Λ(n) w(n) n^{−½−iT}/log n, and1125 S(T) = (1/π) Im[Σ^χ_w(½+iT) − Σ^χ_w(½)] — the comb side of the twisted count, the difference from height 01126 removing the constant of the complex case (Obs. chi). Exposes N and S(T) so that the scan of the shared primitives applies."""1127 def __init__(self, N, chi_table):1128 self.N = int(N)1129 n, lam = von_mangoldt_support(self.N)1130 logn = np.log(n); L = math.log(self.N)1131 w = 0.5 * (1.0 + np.cos(np.pi * logn / L))1132 self.logn = logn1133 self.coef = lam * w / (np.sqrt(n) * logn) * chi_table[(n % len(chi_table)).astype(int)]1134 self.sigma0 = np.sum(self.coef)1135 def sigma_w(self, T):1136 T = np.atleast_1d(np.asarray(T, float))1137 out = np.empty(len(T), dtype=complex)1138 for s in range(0, len(T), 64):1139 blk = T[s:s + 64]1140 out[s:s + 64] = (self.coef[:, None] * np.exp(-1j * np.outer(self.logn, blk))).sum(0)1141 return out1142 def S(self, T):1143 T = np.atleast_1d(np.asarray(T, float))1144 val = (self.sigma_w(T) - self.sigma0).imag / math.pi1145 return val if val.size > 1 else float(val[0])11461147def theta_chi5(T):1148 T = np.asarray(T, float)1149 return np.imag(loggamma(0.75 + 0.5j * T)) + 0.5 * T * math.log(5 / math.pi)11501151def theta_chi5_from0(T):1152 return theta_chi5(T) - theta_chi5(0.0)11531154def L_chi5(s):1155 return mp.power(5, -s) * sum(MP_CHI5[a] * mp.zeta(s, mp.mpf(a) / 5) for a in (1, 2, 3, 4))11561157def Z_chi5_factory():1158 """Z_χ(t) = W(χ)^{−1/2} Λ(½ + it, χ) with Λ(s, χ) = (5/π)^{(s+1)/2} Γ((s+1)/2) L(s, χ) and W(χ) = τ(χ)/(i√5): real on the line."""1159 I = mp.mpc(0, 1)1160 tau = sum(MP_CHI5[a] * mp.expjpi(mp.mpf(2 * a) / 5) for a in (1, 2, 3, 4))1161 W = tau / (I * mp.sqrt(5)); Whalf = mp.sqrt(W)1162 def Z(t):1163 s = mp.mpf(1) / 2 + I * t1164 return mp.power(5 / mp.pi, (s + 1) / 2) * mp.gamma((s + 1) / 2) * L_chi5(s) / Whalf1165 return Z, W11661167def _C_chi_complex():1168 print("\n== block C (continued): the complex character χ mod 5, χ(2) = i (Def. chi, Obs. chi; ledger E16) ==")1169 with Timer("Hurwitz-zeta validation arm, χ mod 5"):1170 Z, W = Z_chi5_factory()1171 ts = np.arange(0.0, 40.0, 0.02)1172 zv = [Z(float(t)) for t in ts]1173 im_max = float(max(abs(mp.im(z)) for z in zv))1174 re = np.array([float(mp.re(z)) for z in zv])1175 idx = np.nonzero(re[:-1] * re[1:] < 0)[0]1176 zeros = np.array([brentq(lambda t: float(mp.re(Z(t))), ts[i], ts[i + 1], xtol=1e-10) for i in idx])1177 N = 10 ** 61178 with Timer(f"complex twisted comb to {N}"):1179 comb = ComplexTwistedComb(N, CHI5)1180 mids = 0.5 * (zeros[1:] + zeros[:-1])1181 cnt = count_raw(comb, mids, const=0.0, theta_fn=theta_chi5_from0)1182 ok_count = len(zeros) == 15 and np.max(np.abs(cnt - np.arange(1, len(zeros)))) < 0.051183 ok_real = abs(float(abs(W)) - 1.0) < 1e-15 and im_max < 1e-151184 # the brackets come from the twisted count itself; the validation arm's zeros enter only in the error reported1185 roots, levels, extra = secular_roots_scan(comb, 2.0, 40.0, const=0.0, theta_fn=theta_chi5_from0)1186 assert levels == [k + 0.5 for k in range(15)] and not extra, f"crossings {levels}, extra {extra}"1187 err = np.abs(roots - zeros)1188 # illustration: 20:E16 retired 2026-10-15 (T05)1189 check("C7c", "χ mod 5 complex: Z_χ real on the line; the twisted count reads 1 … 14 between the fifteen zeros below 40; "1190 "the secular condition (N = 10⁶ comb, brackets from the count alone, L never evaluated) recovers the fifteen heights to mean error 8.0e-5 (max 2.0e-4)",1191 ok_real and ok_count and abs(err.mean() - 8.0e-5) < 1.5e-5 and err.max() < 2.5e-4,1192 f"|W| = {float(abs(W)):.15f}, max|Im Z_χ| = {im_max:.1e}; {len(zeros)} zeros {zeros[0]:.4f} … {zeros[-1]:.4f}; count at midpoints "1193 + ", ".join(f"{c:.3f}" for c in cnt) + "; roots " + ", ".join(f"{r:.4f}" for r in roots[:3]) + f" … {roots[-1]:.4f}; mean error {err.mean():.2e}, max {err.max():.2e}",1194 kind="[approx]")11951196def block_C():1197 """Block C — the spectrum from the prime side with ζ never evaluated (the scale-spectrum, the secular condition raw and pole-corrected to depth 10⁸, the colleague and Jacobi matrices, the additive injection), the GUE statistics of the target spectrum, the χ-twisted comb for χ₋₄ and for the complex character mod 5; 12 checks ([approx], [chart])."""1198 _C_primeside(); _C_gue(); _C_chi()11991200# ======================================================================================================================1201# block D — d_classification.py (until 24 September 2026), its former docstring:1202# d_classification.py — block D: the classification and the Euler-product discriminator (00:D12)1203# =============================================================================================1204# Proposition turing: C(T) ≤ N_crit(T) ≤ N(T); RH below T is N = N_crit; C = N is Turing's practical1205# certificate; the shell reads N from the comb (raw count, uncertified) and C from the de-framing1206# (sign changes of the main sum). Definition screen: the classical hypothesis is read as the screen1207# value N − N_crit = 0. Numerical Observation dh: the value is discriminating — the Davenport–1208# Heilbronn function f (no Euler product) opens a phantom deficit of 2 at its first off-line pair1209# while its constituent L(s, χ) closes with deficit 0; on the comb side the secular condition1210# lands two "on-line" roots on the phantom height and the raw count near it does not settle with1211# depth; the ζ comb behaves the same way at ζ's own singularity to the right of the line, the pole.1212#1213# Paper-local predicates:1214#1215# D1 Prop. turing / Def. screen [approx] at T = 15, 30, 50.3 the raw comb count Ñ_N rounds to the exact1216# N = 1, 3, 10, and the sign changes of the horizon main sum Z_M on1217# (0, T) — the de-framing side — give C = 1, 3, 10 = N (Turing's1218# certificate read on the shell, both readings uncertified)1219# D2a Obs. dh [approx] Λ_f(½ + it) is real to 10⁻²⁰ (self-dual normalisation of the DH function)1220# D2b Obs. dh [approx] argument principle per band on [0, 87]: N_f = 45 strip zeros against C_f = 43 on-line sign1221# changes, the deficit 2 opening only in the band holding the off-line pair1222# s = 0.8085171825 + 85.6993484854 i (validated |f(s)| < 10⁻⁶); |Z_f| dips to 0.357 at1223# t = 85.71 without crossing; positive control L(s, χ): 45 = 45, deficit 0 in every band1224# D2c Obs. dh [approx] negative control: the secular condition on the Λ_f comb returns two "on-line" roots in1225# [84, 87], 85.63 / 85.76 (comb to 10⁵) and 85.65 / 85.75 (4×10⁵), at the phantom height —1226# the upward crossings of the levels 43.5 and 44.5 between the exact counts 43 and 45.1227# (Package finding: at 4×10⁵ the unsettled count also swings below 42.5 and above 45.5 in1228# [85.3, 86.1], adding two crossing pairs the paper's sentence does not mention; at 10⁵ the1229# excursion stays within (42.5, 45.5) and the two roots are the only crossings.)1230# D2d Obs. dh [approx] the raw DH comb count at t = 85.9 drifts monotonically with depth, 45.14 → 45.73 over1231# 5×10⁴ … 8×10⁵ (exact 45); at 85.3 it reads 42.66 … 42.73 (exact 43; the paper's two figures1232# are the endpoint depths, the interior depths dip to 42.64); away from the pair it is1233# stable: 28.005 ± 0.002 at 60.3 and 43.00 ± 0.01 at 84.01234# D2e Obs. dh [approx] the ζ comb at the pole: at t = 1 (exact 0) the raw count reads +0.09, −0.26, +0.36, −0.08,1235# −0.81, −1.42, −1.46 at depths 10⁴ … 8×10⁷; at t = 5 within ±0.02; at t = 10 within ±0.003;1236# at t = 15 within 2×10⁻³ of 1 at every depth; at t = 30 from 3.04 to 2.998; at 50.3 within1237# 10⁻² of 10 (depths to 8×10⁵)1238# D2f Prop. combformula, Obs. dh [approx] the pole-corrected ζ count N̂_N at t = 1 reads 0.0051, 0.0032, 0.0022, 0.0016, 0.0015,1239# 0.0013, 0.0012 at the same depths (exact 0), falling monotonically; at t = 5 and 10 within 4×10⁻⁴ from1240# 10⁶ on, falling in magnitude; at 15, 30, 50.3 the correction is below 10⁻³ and raw and corrected agree1241# D2g Rem. combsettle, Obs. dh [approx] the DH drift is the zero term of the off-line zero ρ₀: the count corrected by it and its constant reads1242# 44.960, 44.964, 44.967, 44.970, 44.972 at t = 85.9 (raw 45.14 … 45.73; exact 45) and 43.035 → 43.0221243# at 85.3 (raw 42.66 … 42.73; exact 43), both monotone toward the exact count1244# D3 Prop. combformula [approx] the identity's validation arm: the corrected tapered sum Σ_w − Π_N − log((s−1)/s) against1245# log ζ(½ + it) (mpmath) at t = 1, 5, 10, 15, 30: the real-part error falls with depth at every height1246# (0.0049 at t = 1, depth 8×10⁷) and the imaginary part is the corrected count of D2f1247# D4 Prop. turing (frame-exact inputs) [approx] on the shells p = 97, 1009, 4801 (primes ≡ 1 mod 4; ceilings T = 2πp ≈ 609, 6340,1248# 30166) the raw count from the frame-exact comb of depth ⌊√p⌋ = 9, 31, 69, evaluated midway between1249# consecutive zeros just below the ceiling (the test heights chosen with the validation arm: zeros from1250# Z(t) sign changes, the count from mpmath nzeros), rounds to the exact N(t) at every midpoint: maximum1251# deviation 0.079, 0.107, 0.164 (within 0.08, 0.11, 0.17); the 10⁶ comb within 0.0021252#1253# The "ζ never evaluated" discipline does not bind the control arms (Hurwitz zeta functions are1254# evaluated there). Figure: fig_dh.pdf. Master ledger: 00:D12.12551256mp.mp.dps = 201257I = mp.mpc(0, 1)1258CHI = {1: mp.mpc(1), 2: I, 3: -I, 4: mp.mpc(-1)} # complex character mod 5, χ(2) = i1259KAPPA = (mp.sqrt(10 - 2 * mp.sqrt(5)) - 2) / (mp.sqrt(5) - 1)1260TAU = sum(CHI[a] * mp.expjpi(mp.mpf(2 * a) / 5) for a in (1, 2, 3, 4))1261EPS_HALF = mp.sqrt(TAU / (I * mp.sqrt(5))) # half-phase of the root number of χ1262RHO_OFF = mp.mpc("0.8085171825", "85.6993484854") # Titchmarsh: the first off-line pair, with 1 − s̄1263T_HI = 87.01264DEPTHS_DH = [50_000, 100_000, 200_000, 400_000, 800_000]1265DEPTHS_Z = [10 ** 4, 10 ** 5, 10 ** 6, 10 ** 7, 2 * 10 ** 7, 4 * 10 ** 7, 8 * 10 ** 7] if not FAST else [10 ** 4, 10 ** 5, 10 ** 6, 10 ** 7]1266STATED_T1 = [0.09, -0.26, 0.36, -0.08, -0.81, -1.42, -1.46]12671268# ------------------------------------------------------------------ the validation arm: f and L(s, χ) via Hurwitz zeta1269def L_pair(s):1270 z = {a: mp.zeta(s, mp.mpf(a) / 5) for a in (1, 2, 3, 4)}1271 five = mp.power(5, -s)1272 return five * sum(CHI[a] * z[a] for a in (1, 2, 3, 4)), five * sum(mp.conj(CHI[a]) * z[a] for a in (1, 2, 3, 4))12731274def pair_values(s):1275 """(L(s,χ), f(s)) with f = (1 − iκ)/2 L(s,χ) + (1 + iκ)/2 L(s,χ̄)."""1276 Lchi, Lbar = L_pair(s)1277 return Lchi, (1 - I * KAPPA) / 2 * Lchi + (1 + I * KAPPA) / 2 * Lbar12781279def gamma_factor(s):1280 return mp.power(5 / mp.pi, (s + 1) / 2) * mp.gamma((s + 1) / 2)12811282def hardy_pair(t):1283 """Hardy-normalised on-line readouts (Z_χ(t), Z_f(t)), both real: |Z| = |value on the line|."""1284 s = mp.mpf(1) / 2 + I * t1285 g = gamma_factor(s); ph = g / abs(g)1286 Lchi, f = pair_values(s)1287 return ph * Lchi / EPS_HALF, ph * f12881289def theta_f(t):1290 """Twisted archimedean phase of the conductor-5 odd Γ-factor: Im log Γ(¾ + it/2) + (t/2) log(5/π)."""1291 t = np.asarray(t, float)1292 return np.imag(loggamma(0.75 + 0.5j * t)) + 0.5 * t * math.log(5 / math.pi)12931294def sign_changes(ts, vals):1295 pos = []1296 for k in range(1, len(ts)):1297 a, b = vals[k - 1], vals[k]1298 if a * b < 0:1299 pos.append(ts[k - 1] + (ts[k] - ts[k - 1]) * a / (a - b))1300 return np.array(pos)13011302def winding(fn_values, corners, step):1303 """Winding numbers of several functions around a rectangle, by continuous phase tracking with1304 adaptive refinement; fn_values(s) returns a tuple of complex values (one per function)."""1305 (s0, s1, s2, s3) = corners1306 sides = [(s0, s1), (s1, s2), (s2, s3), (s3, s0)]1307 total = None1308 for a, b in sides:1309 n = max(2, int(abs(b - a) / step) + 1)1310 pts = [a + (b - a) * k / (n - 1) for k in range(n)]1311 vals = [fn_values(s) for s in pts]1312 m = len(vals[0])1313 acc = np.zeros(m)1314 stack = list(zip(pts[:-1], pts[1:], vals[:-1], vals[1:], [0] * (n - 1)))1315 while stack:1316 pa, pb, va, vb, depth = stack.pop()1317 d = np.array([float(mp.arg(vb[i] / va[i])) for i in range(m)])1318 if np.max(np.abs(d)) > 0.8 and depth < 12:1319 pm = (pa + pb) / 2; vm = fn_values(pm)1320 stack.append((pa, pm, va, vm, depth + 1)); stack.append((pm, pb, vm, vb, depth + 1))1321 else:1322 acc += d1323 total = acc if total is None else total + acc1324 w = total / (2 * math.pi)1325 assert np.max(np.abs(w - np.round(w))) < 0.05, f"non-integer winding {w}"1326 return np.round(w).astype(int)13271328# ------------------------------------------------------------------ the generalized von Mangoldt comb of f1329def dh_lambda(M):1330 """Λ_f(n) for n ≤ M from −f'/f = Σ Λ_f(n) n^{−s}, by the divisor recursion a_n log n = Σ_{d|n} Λ_f(d) a_{n/d}, a_1 = 1."""1331 kappa = float(KAPPA)1332 n = np.arange(M + 1)1333 cvals = np.array([0, 1, 1j, -1j, -1], dtype=complex)[n % 5]1334 a = ((1 - 1j * kappa) / 2 * cvals + (1 + 1j * kappa) / 2 * np.conj(cvals)).real1335 logs = np.log(np.maximum(n, 1).astype(float))1336 Lf = np.zeros(M + 1); conv = np.zeros(M + 1)1337 for k in range(2, M + 1):1338 Lf[k] = a[k] * logs[k] - conv[k]1339 if Lf[k] != 0.0:1340 ms = np.arange(2, M // k + 1)1341 conv[k * ms] += Lf[k] * a[ms]1342 return Lf13431344def comb_support_from(Lf):1345 def lam_fn(N):1346 idx = np.nonzero(Lf[:N + 1])[0]1347 return idx.astype(float), Lf[idx]1348 return lam_fn13491350def count_at_depths(n_all, lam_all, t, depths, theta_fn, const):1351 """Raw comb count at height t for several depths from one support (taper log N per depth)."""1352 logn = np.log(n_all); out = []1353 for N in depths:1354 m = n_all <= N1355 w = 0.5 * (1 + np.cos(np.pi * logn[m] / math.log(N)))1356 terms = lam_all[m] / (np.sqrt(n_all[m]) * logn[m]) * w * np.sin(t * logn[m])1357 out.append(float(theta_fn(t)) / math.pi + const - math.fsum(terms) / math.pi)1358 return np.array(out)13591360def corrected_at_depths(n_all, lam_all, t, depths):1361 """Pole-corrected ζ count N̂_N(t) at several depths from one support (Prop. combformula)."""1362 raw = count_at_depths(n_all, lam_all, t, depths, theta, 1.0)1363 return raw - np.array([pole_count_term(N, t) for N in depths])13641365def zero_cut_term(N, t, rho0):1366 """The zero term of one fixed zero ρ₀ = β₀ + iγ₀ to the right of the line (½ < β₀ < 3/2) at s = ½ + it1367 (Remark combsettle): Π^{(ρ₀)}_N(s) = ∫_0^L W(u/L) (e^{(ρ₀−s)u} − e^{(ρ₀−1−s)u}) du/u, L = log N. The image of the1368 regulariser h_{ρ₀}(w) = log((w−ρ₀)/(w−ρ₀+1)) under the window inversion of Prop. combformula (i) is −Π^{(ρ₀)}_N,1369 the term the zero contributes to the smoothed sum. The routine evaluates this one integral at depth N and1370 asserts nothing about a limit N → ∞: the paper does not extend the identity of Prop. combformula to a series1371 with infinitely many zeros to the right of the line. The Davenport–Heilbronn correction1372 N̂_{f,N} = Ñ_{f,N} + Im[Π^{(ρ₀)}_N(s) + h_{ρ₀}(s)]/π (dh_corrected_at_depths, D2g) is validated over the depths1373 listed there."""1374 L = math.log(N); a = rho0 - mp.mpc(0.5, t); b = a - 11375 f = lambda u: 0.5 * (1 + mp.cos(mp.pi * u / L)) * (mp.exp(a * u) - mp.exp(b * u)) / u1376 return mp.quad(f, mp.linspace(0, L, 40))13771378def dh_corrected_at_depths(n_all, lam_all, t, depths):1379 """The DH comb count corrected by the zero term of ρ₀ and its constant (D2g; Remark combsettle):1380 N̂_{f,N} = Ñ_{f,N} + Im[Π^{(ρ₀)}_N(s) + h_{ρ₀}(s)]/π, with h_{ρ₀}(s) = −H0 below."""1381 s = mp.mpc(0.5, t); H0 = mp.log((s - RHO_OFF + 1) / (s - RHO_OFF))1382 raw = count_at_depths(n_all, lam_all, t, depths, theta_f, 0.0)1383 return raw + np.array([float(mp.im(zero_cut_term(N, t, RHO_OFF) - H0)) / math.pi for N in depths])13841385def Z_main(t):1386 """The horizon-length Riemann–Siegel main sum (Definition deframe), the de-framing side of Proposition turing."""1387 N = int(math.floor(math.sqrt(t / (2 * math.pi))))1388 th = float(theta(t))1389 return 2 * sum(math.cos(th - t * math.log(n)) / math.sqrt(n) for n in range(1, N + 1))13901391def block_D():1392 """Block D — Turing's count on the shell (T = 15, 30, 50.3 and the frame-exact combs of p = 97, 1009, 4801), the phantom of the Davenport–Heilbronn function on the line and on the comb side, the ζ comb at its pole raw and corrected, the smoothed explicit formula's validation arm; 10 checks ([approx])."""1393 mp.mp.dps = 201394 print("\n== block D: the screen reading and its Davenport–Heilbronn control (Prop. turing, Def. screen, Obs. dh) ==")1395 # ---------------- D1: the two readings on ζ1396 comb = Comb(10 ** 6)1397 Tq = np.array([15.0, 30.0, 50.3]); exact = [1, 3, 10]1398 raw = count_raw(comb, Tq)1399 tgrid = np.arange(0.5, 50.3, 0.02)1400 zm = np.array([Z_main(t) for t in tgrid])1401 Cm = [int(np.sum(zm[:-1][tgrid[:-1] < T] * zm[1:][tgrid[:-1] < T] < 0)) for T in Tq]1402 ok = all(int(round(r)) == e for r, e in zip(raw, exact)) and Cm == exact1403 # 20:F1 (p20047)1404 check("D1", "Ñ_N rounds to N = 1, 3, 10 at T = 15, 30, 50.3; sign changes of the main sum Z_M give C = 1, 3, 10 = N", ok,1405 "raw " + ", ".join(f"{r:.4f}" for r in raw) + f"; C_M = {Cm}", kind="[approx]")1406 # ---------------- D2a: reality of Λ_f on the line1407 with mp.workdps(30):1408 im = max(abs(mp.im(gamma_factor(mp.mpf(1) / 2 + I * t) * pair_values(mp.mpf(1) / 2 + I * t)[1])) for t in np.linspace(0.7, 86.3, 25))1409 # illustration: 20:F3 retired 2026-10-15 (T11)1410 check("D2a", "the completed Λ_f(½ + it) is real on the line to 10⁻²⁰", im < 1e-20, f"max |Im Λ_f| = {float(im):.1e} over 25 heights", kind="[approx]")1411 # ---------------- D2b: strip counts by the argument principle against on-line sign changes1412 with Timer("on-line grid (Hurwitz zeta)"):1413 dt = 0.02 if not FAST else 0.041414 ts = np.arange(0.1, T_HI + dt / 2, dt)1415 Zc = np.empty_like(ts); Zf = np.empty_like(ts)1416 for k, t in enumerate(ts):1417 zc, zf = hardy_pair(mp.mpf(float(t)))1418 assert abs(mp.im(zc)) < 1e-9 and abs(mp.im(zf)) < 1e-91419 Zc[k], Zf[k] = float(mp.re(zc)), float(mp.re(zf))1420 zeros_f, zeros_c = sign_changes(ts, Zf), sign_changes(ts, Zc)1421 off_err = abs(pair_values(RHO_OFF)[1])1422 band = (ts > 85.0) & (ts < 86.5); k_dip = np.argmin(np.abs(Zf[band])); t_dip, z_dip = ts[band][k_dip], abs(Zf[band][k_dip])1423 # band boundaries: multiples of 3, nudged away from every on-line zero of both functions1424 allz = np.concatenate([zeros_f, zeros_c])1425 bounds = [0.1]1426 for b in np.arange(3.0, T_HI + 0.1, 3.0):1427 while np.min(np.abs(allz - b)) < 0.15: b += 0.21428 bounds.append(min(b, T_HI))1429 with Timer(f"argument principle on {len(bounds) - 1} bands (Hurwitz zeta)"):1430 step = 0.05 if not FAST else 0.11431 Nf_band, Nc_band = [], []1432 for lo, hi in zip(bounds[:-1], bounds[1:]):1433 corners = (mp.mpc(-1.5, lo), mp.mpc(2.5, lo), mp.mpc(2.5, hi), mp.mpc(-1.5, hi))1434 wc, wf = winding(pair_values, corners, step)1435 Nc_band.append(wc); Nf_band.append(wf)1436 Cf_band = [int(np.sum((zeros_f > lo) & (zeros_f < hi))) for lo, hi in zip(bounds[:-1], bounds[1:])]1437 Cc_band = [int(np.sum((zeros_c > lo) & (zeros_c < hi))) for lo, hi in zip(bounds[:-1], bounds[1:])]1438 deficit_f = np.array(Nf_band) - np.array(Cf_band); deficit_c = np.array(Nc_band) - np.array(Cc_band)1439 kband = next(i for i, (lo, hi) in enumerate(zip(bounds[:-1], bounds[1:])) if lo < 85.6993 < hi)1440 ok = (sum(Nf_band) == 45 and len(zeros_f) == 43 and sum(Nc_band) == 45 and len(zeros_c) == 451441 and deficit_f[kband] == 2 and all(d == 0 for i, d in enumerate(deficit_f) if i != kband) and all(d == 0 for d in deficit_c)1442 and off_err < 1e-6 and abs(z_dip - 0.357) < 0.01 and abs(t_dip - 85.71) < 0.03)1443 check("D2b", "argument principle on [0,87]: f: 45 strip zeros vs 43 on-line, deficit 2 only in the band of the off-line pair; L(s,χ): 45 = 45, deficit 0 everywhere", ok,1444 f"f: N = {sum(Nf_band)}, C = {len(zeros_f)}, deficit by band {[int(d) for d in deficit_f]}; L(s,χ): N = {sum(Nc_band)}, C = {len(zeros_c)}; "1445 f"|f(s_off)| = {float(off_err):.1e}; |Z_f| dips to {z_dip:.3f} at t = {t_dip:.2f}", kind="[approx]")1446 # ---------------- D2c, D2d: the comb side of f1447 with Timer("Λ_f comb to 8e5 (divisor recursion)"):1448 Lf = dh_lambda(DEPTHS_DH[-1])1449 lam_fn = comb_support_from(Lf)1450 roots, extra, excursion = {}, {}, {}1451 for M in (10 ** 5, 4 * 10 ** 5):1452 c = Comb(M, lam_fn=lam_fn)1453 tt = np.linspace(84, 87, 601); v = count_raw(c, tt, const=0.0, theta_fn=theta_f)1454 cross = [] # every half-integer crossing (t, level, direction)1455 for i in range(len(tt) - 1):1456 lo, hi = v[i], v[i + 1]1457 for h in np.arange(math.floor(min(lo, hi)) + 0.5, max(lo, hi), 1.0):1458 if (lo - h) * (hi - h) < 0:1459 cross.append((brentq(lambda t: count_raw(c, t, const=0.0, theta_fn=theta_f) - h, tt[i], tt[i + 1]), h, np.sign(hi - lo)))1460 # the phantom roots: the upward crossings of the two levels 43.5, 44.5 between the exact counts 43 (below the pair) and 45 (above)1461 roots[M] = sorted(t for t, h, d in cross if h in (43.5, 44.5) and d > 0)1462 extra[M] = sorted((round(t, 2), h) for t, h, d in cross if h not in (43.5, 44.5))1463 excursion[M] = (v[(tt > 85) & (tt < 86.4)].min(), v[(tt > 85) & (tt < 86.4)].max())1464 ok = (len(roots[10 ** 5]) == 2 and len(roots[4 * 10 ** 5]) == 2 and len(extra[10 ** 5]) == 01465 and np.allclose(roots[10 ** 5], [85.63, 85.76], atol=0.015) and np.allclose(roots[4 * 10 ** 5], [85.65, 85.75], atol=0.015))1466 # illustration: 20:F4 retired 2026-10-15 (T11)1467 check("D2c", "the secular condition on the Λ_f comb returns two 'on-line' roots at the phantom height: 85.63/85.76 (10⁵), 85.65/85.75 (4×10⁵)", ok,1468 f"10⁵: {np.round(roots[10 ** 5], 3).tolist()} (no other crossing in [84,87]; count excursion {excursion[10 ** 5][0]:.2f}–{excursion[10 ** 5][1]:.2f}); "1469 f"4×10⁵: {np.round(roots[4 * 10 ** 5], 3).tolist()}, plus the swing of the unsettled count through the levels {sorted(set(float(h) for _, h in extra[4 * 10 ** 5]))} "1470 f"(excursion {excursion[4 * 10 ** 5][0]:.2f}–{excursion[4 * 10 ** 5][1]:.2f}; crossing pairs at {[t for t, _ in extra[4 * 10 ** 5]]})", kind="[approx]")1471 idx = np.nonzero(Lf)[0]; nf, lf = idx.astype(float), Lf[idx]1472 scan = {t: count_at_depths(nf, lf, t, DEPTHS_DH, theta_f, 0.0) for t in (85.9, 85.3, 60.3, 84.0)}1473 s859 = scan[85.9]1474 ok = (np.allclose(s859, [45.14, 45.27, 45.41, 45.57, 45.73], atol=0.015) and np.all(np.diff(s859) > 0)1475 and abs(scan[85.3][0] - 42.66) < 0.006 and abs(scan[85.3][-1] - 42.73) < 0.006 and scan[85.3].min() > 42.63 and scan[85.3].max() < 42.741476 and np.all(np.abs(scan[60.3] - 28.005) < 0.003) and np.all(np.abs(scan[84.0] - 43.0) < 0.011))1477 check("D2d", "raw DH comb count drifts monotonically at t = 85.9 (45.14 → 45.73) and reads 42.66 … 42.73 at 85.3 (interior depths dip to 42.64); stable at 60.3 (28.005 ± 0.002) and 84.0 (43.00 ± 0.01)", ok,1478 "; ".join(f"t={t}: " + ", ".join(f"{v:.3f}" for v in vals) for t, vals in scan.items()), kind="[approx]")1479 # ---------------- D2e: the ζ comb at the pole1480 with Timer(f"ζ comb support to {DEPTHS_Z[-1]:.0e}"):1481 n_all, lam_all = von_mangoldt_support(DEPTHS_Z[-1])1482 z1 = count_at_depths(n_all, lam_all, 1.0, DEPTHS_Z, theta, 1.0)1483 z5 = count_at_depths(n_all, lam_all, 5.0, DEPTHS_Z, theta, 1.0)1484 z10 = count_at_depths(n_all, lam_all, 10.0, DEPTHS_Z, theta, 1.0)1485 z15 = count_at_depths(n_all, lam_all, 15.0, DEPTHS_Z, theta, 1.0)1486 z30 = count_at_depths(n_all, lam_all, 30.0, DEPTHS_Z, theta, 1.0)1487 z503 = count_at_depths(n_all, lam_all, 50.3, [d for d in DEPTHS_Z if d <= 8 * 10 ** 5] + [8 * 10 ** 5], theta, 1.0)1488 k = len(DEPTHS_Z)1489 ok = (np.allclose(z1, STATED_T1[:k], atol=0.015) and np.all(np.abs(z5) < 0.021) and np.all(np.abs(z10) < 0.0031)1490 and np.all(np.abs(z15 - 1) < 2e-3) and abs(z30[0] - 3.04) < 0.006 and abs(z30[-1] - (2.998 if not FAST else z30[-1])) < 0.0021491 and np.all(np.abs(z503 - 10) < 1e-2))1492 # illustration: 20:F5 retired 2026-10-15 (T11)1493 check("D2e", "ζ comb at the pole: t = 1 reads +0.09, −0.26, +0.36, −0.08, −0.81, −1.42, −1.46 with depth (exact 0); t = 5, 10 within ±0.02, ±0.003; t = 15 within 2e-3 of 1; t = 30 3.04 → 2.998; t = 50.3 within 1e-2 of 10", ok,1494 "t=1: " + ", ".join(f"{v:+.2f}" for v in z1) + "; t=5: max|·| " + f"{np.max(np.abs(z5)):.3f}" + "; t=10: " + f"{np.max(np.abs(z10)):.4f}"1495 + "; t=15: max dev " + f"{np.max(np.abs(z15 - 1)):.1e}" + f"; t=30: {z30[0]:.3f} → {z30[-1]:.3f}; t=50.3: max dev {np.max(np.abs(z503 - 10)):.1e}", kind="[approx]")1496 # ---------------- D2f: the pole-corrected ζ count settles (Prop. combformula)1497 c1 = corrected_at_depths(n_all, lam_all, 1.0, DEPTHS_Z); c5 = corrected_at_depths(n_all, lam_all, 5.0, DEPTHS_Z); c10 = corrected_at_depths(n_all, lam_all, 10.0, DEPTHS_Z)1498 stated_c1 = [0.0051, 0.0032, 0.0022, 0.0016, 0.0015, 0.0013, 0.0012][:k]1499 i6 = DEPTHS_Z.index(10 ** 6)1500 ok = (np.allclose(c1, stated_c1, atol=0.0003) and np.all(np.diff(c1) < 0) and np.all(np.abs(c5[i6:]) < 4e-4) and np.all(np.abs(c10[i6:]) < 4e-4)1501 and np.all(np.diff(np.abs(c5)) < 0) and np.all(np.diff(np.abs(c10)) < 0))1502 # illustration: 20:F8 retired 2026-10-15 (T11)1503 check("D2f", "pole-corrected ζ count at t = 1: 0.0051, 0.0032, 0.0022, 0.0016, 0.0015, 0.0013, 0.0012 (exact 0), monotone; t = 5, 10 within 4e-4 from 10⁶ on, falling", ok,1504 "t=1: " + ", ".join(f"{v:+.4f}" for v in c1) + "; t=5: " + ", ".join(f"{v:+.5f}" for v in c5) + "; t=10: " + ", ".join(f"{v:+.5f}" for v in c10), kind="[approx]")1505 # ---------------- D2g: the DH drift is the zero term of the off-line zero (Remark combsettle)1506 with Timer("DH corrected counts"):1507 d859 = dh_corrected_at_depths(nf, lf, 85.9, DEPTHS_DH); d853 = dh_corrected_at_depths(nf, lf, 85.3, DEPTHS_DH)1508 ok = (np.allclose(d859, [44.960, 44.964, 44.967, 44.970, 44.972], atol=0.002) and np.all(np.diff(d859) > 0)1509 and np.allclose(d853, [43.035, 43.031, 43.028, 43.025, 43.022], atol=0.002) and np.all(np.diff(d853) < 0))1510 check("D2g", "DH count corrected by the zero term of ρ₀ and its constant: 44.960 → 44.972 at 85.9 (raw 45.14 → 45.73; exact 45), 43.035 → 43.022 at 85.3 (exact 43), monotone toward the exact count", ok,1511 "t=85.9: " + ", ".join(f"{v:.4f}" for v in d859) + "; t=85.3: " + ", ".join(f"{v:.4f}" for v in d853), kind="[approx]")1512 # ---------------- D3: the identity against log ζ (validation arm: mpmath)1513 logn_all = np.log(n_all); re_err = {}1514 for t, Nc in ((1.0, 0), (5.0, 0), (10.0, 0), (15.0, 1), (30.0, 3)):1515 sc = mp.mpc(0.5, t); lz_re = float(mp.log(abs(mp.zeta(sc))))1516 errs = []1517 for N in DEPTHS_Z:1518 m = n_all <= N; L = math.log(N); w = 0.5 * (1 + np.cos(np.pi * logn_all[m] / L))1519 S = np.sum(lam_all[m] * w / (np.sqrt(n_all[m]) * logn_all[m]) * np.exp(-1j * t * logn_all[m]))1520 errs.append(abs((S - pole_term(N, t) - complex(mp.log((sc - 1) / sc))).real - lz_re))1521 re_err[t] = np.array(errs)1522 ok = all(np.all(np.diff(e) < 0) for e in re_err.values()) and re_err[1.0][-1] < (0.0052 if not FAST else 0.0065)1523 check("D3", "the identity's validation arm: Re(Σ_w − Π_N − log((s−1)/s)) → log|ζ(½+it)| at t = 1, 5, 10, 15, 30, the error falling with depth (0.0049 at t = 1, 8×10⁷)", ok,1524 "; ".join(f"t={t}: " + " → ".join(f"{v:.4f}" for v in (e[0], e[len(e)//2], e[-1])) for t, e in re_err.items()), kind="[approx]")1525 # ---------------- D4: the count from the frame-exact comb on one shell (Proposition turing)1526 with Timer("frame-exact shells"):1527 shells = [97, 1009] + ([4801] if not FAST else []) # primes ≡ 1 (mod 4): admissible Subject shells1528 dev_fe, dev_deep, dev_smooth, npts = {}, {}, {}, {}1529 deep = Comb(10 ** 6)1530 for p_sh in shells:1531 Tc = 2 * math.pi * p_sh; depth = int(math.floor(math.sqrt(p_sh)))1532 grid = np.arange(Tc - 12.0, Tc + 1e-9, 0.1); zv = np.array([float(mp.siegelz(t)) for t in grid])1533 zs = np.array([brentq(lambda t: float(mp.siegelz(t)), grid[i], grid[i + 1], xtol=1e-7)1534 for i in range(len(grid) - 1) if zv[i] * zv[i + 1] < 0])1535 mids = 0.5 * (zs[1:] + zs[:-1]); nz = np.array([int(mp.nzeros(t)) for t in mids]); npts[p_sh] = len(mids)1536 v = count_raw(Comb(depth), mids); vd = count_raw(deep, mids)1537 vs = np.array([float(theta(t)) / math.pi + 1 for t in mids]) # the control: the smooth term alone1538 dev_fe[p_sh] = (float(np.max(np.abs(v - nz))), bool(np.all(np.round(v) == nz)))1539 dev_deep[p_sh] = (float(np.max(np.abs(vd - nz))), bool(np.all(np.round(vd) == nz)))1540 dev_smooth[p_sh] = (float(np.max(np.abs(vs - nz))), bool(np.all(np.round(vs) == nz)))1541 stated = {97: 0.079, 1009: 0.107, 4801: 0.164}1542 stated_smooth = {97: (0.268, True), 1009: (0.498, True), 4801: (0.523, False)}1543 ok = (all(r for _, r in dev_fe.values()) and all(r for _, r in dev_deep.values())1544 and all(abs(dev_fe[p_sh][0] - stated[p_sh]) < 0.02 for p_sh in shells) and max(d for d, _ in dev_deep.values()) < 0.0041545 and all(abs(dev_smooth[p_sh][0] - stated_smooth[p_sh][0]) < 0.02 and dev_smooth[p_sh][1] == stated_smooth[p_sh][1] for p_sh in shells))1546 check("D4", "frame-exact count on one shell: at p = 97, 1009, 4801 (T = 2πp) the depth-⌊√p⌋ comb, midway between consecutive zeros below the ceiling, rounds to the exact N(t) at every midpoint (max deviation 0.079, 0.107, 0.164 — within 0.08, 0.11, 0.17); the 10⁶ comb within 0.002; the control θ(t)/π + 1 alone deviates by 0.27, 0.498, 0.52 and fails to round at one midpoint of p = 4801", ok,1547 "; ".join(f"p={p_sh} (depth {int(math.sqrt(p_sh))}, {npts[p_sh]} midpoints): max dev {dev_fe[p_sh][0]:.3f} (10⁶ comb {dev_deep[p_sh][0]:.3f}; smooth term alone {dev_smooth[p_sh][0]:.3f}, rounds {dev_smooth[p_sh][1]})" for p_sh in shells), kind="[approx]")1548 # ---------------- figure1549 with Timer("window grid [84,87]"):1550 tw = np.arange(84.0, 87.0 + 1e-9, 0.005 if not FAST else 0.01)1551 Zcw = np.empty_like(tw); Zfw = np.empty_like(tw)1552 for j, t in enumerate(tw):1553 zc, zf = hardy_pair(mp.mpf(float(t))); Zcw[j], Zfw[j] = float(mp.re(zc)), float(mp.re(zf))1554 crimson, gray = "#c8102e", "0.40"1555 fig, ax = plt.subplots(1, 3, figsize=(16, 4.3), gridspec_kw={"width_ratios": [1.15, 1.0, 1.0]})1556 ax[0].axhline(0, color="0.65", lw=0.8)1557 ax[0].plot(tw, Zcw, color=gray, lw=1.6, label=r"$Z_\chi(t)$: constituent $L(s,\chi)$ (crosses)")1558 ax[0].plot(tw, Zfw, color=crimson, lw=2.0, label=r"$Z_f(t)$: Davenport–Heilbronn (dips, no crossing)")1559 ax[0].axvline(float(mp.im(RHO_OFF)), color=crimson, lw=1.0, ls=":")1560 ax[0].set_xlim(84, 87); ax[0].set_ylim(-7.2, 1.7); ax[0].set_xlabel("$t$"); ax[0].set_ylabel("rotated on-line value")1561 ax[0].set_title(f"D2b: the phantom on the line — dip {z_dip:.3f} at t = {t_dip:.2f}", fontsize=10); ax[0].legend(loc="lower left", fontsize=8)1562 # staircases: strip count from the band windings, on-line count from the sign changes1563 xs_on = np.concatenate([[0], np.sort(zeros_f), [T_HI]]); ys_on = np.concatenate([[0], np.arange(1, len(zeros_f) + 1), [len(zeros_f)]])1564 # strip count: on-line zeros plus the band excess placed at the off-line height1565 jumps = sorted([(z, 1) for z in zeros_f] + [(float(mp.im(RHO_OFF)), int(deficit_f[kband]))])1566 xs_s, ys_s, c = [0.0], [0], 01567 for z, j in jumps: xs_s.append(z); c += j; ys_s.append(c)1568 xs_s.append(T_HI); ys_s.append(c)1569 ax[1].step(xs_s, ys_s, where="post", color=crimson, lw=1.9, label=r"$N_f(t)$: strip count (argument principle by band)")1570 ax[1].step(xs_on, ys_on, where="post", color=gray, lw=1.4, ls="--", label=r"$N_{\rm crit}(t)$: on-line sign changes")1571 ax[1].axvline(float(mp.im(RHO_OFF)), color=crimson, lw=1.0, ls=":")1572 ax[1].set_xlim(0, T_HI); ax[1].set_ylim(0, 47.5); ax[1].set_xlabel("$t$"); ax[1].set_ylabel("cumulative zero count")1573 ax[1].set_title(f"deficit on [0,87]: $N_f={sum(Nf_band)}$, $C={len(zeros_f)}$; control $L(s,\\chi)$: {sum(Nc_band)}={len(zeros_c)}", fontsize=10)1574 ax[1].legend(loc="upper left", fontsize=8)1575 axins = ax[1].inset_axes([0.60, 0.08, 0.38, 0.40])1576 axins.step(xs_s, ys_s, where="post", color=crimson, lw=1.9); axins.step(xs_on, ys_on, where="post", color=gray, lw=1.4, ls="--")1577 axins.set_xlim(84, 87); axins.set_ylim(40, 46.5); axins.tick_params(labelsize=7); axins.set_title(r"$N_f-N_{\rm crit}=2$", fontsize=8)1578 # depth drift1579 ax[2].semilogx(DEPTHS_DH, scan[85.9], "o-", color=crimson, label="DH comb, t = 85.9 (exact 45)")1580 ax[2].semilogx(DEPTHS_DH, scan[84.0], "s-", color="#e59866", label="DH comb, t = 84.0 (exact 43)")1581 ax[2].axhline(45, color=crimson, lw=0.8, ls=":"); ax[2].axhline(43, color="#e59866", lw=0.8, ls=":")1582 ax2 = ax[2].twinx()1583 ax2.semilogx(DEPTHS_Z, z1, "^-", color="#10325f", label="ζ comb, t = 1 (exact 0)")1584 ax2.semilogx(DEPTHS_Z, c1, "^--", color="#1f6b4a", label="ζ comb, t = 1, pole-corrected")1585 ax2.semilogx(DEPTHS_Z, z15, "v-", color="#5b9bd5", label="ζ comb, t = 15 (exact 1)")1586 ax2.axhline(0, color="#10325f", lw=0.8, ls=":"); ax2.axhline(1, color="#5b9bd5", lw=0.8, ls=":")1587 ax[2].set_xlabel("comb depth $N$"); ax[2].set_ylabel("raw DH count", color=crimson); ax2.set_ylabel("raw ζ count", color="#10325f")1588 ax[2].set_title("D2d–D2f: the raw count near a singularity, and the pole-corrected count", fontsize=10)1589 h1, l1 = ax[2].get_legend_handles_labels(); h2, l2 = ax2.get_legend_handles_labels(); ax[2].legend(h1 + h2, l1 + l2, fontsize=7.5, loc="upper center", bbox_to_anchor=(0.5, -0.17), ncol=2)1590 plt.tight_layout(); plt.savefig(f"{FIGDIR}/fig_dh.pdf", bbox_inches="tight"); plt.savefig(f"{FIGDIR}/fig_dh.png", dpi=110, bbox_inches="tight"); plt.close()1591 print(f" wrote {FIGDIR}/fig_dh.pdf")1592 # the paper's Figure "the Euler-product discriminator": (a) the phantom on the line, (b) the deficit staircase with its inset;1593 # the depth panel above stays a diagnostic of the package1594 toff = float(mp.im(RHO_OFF))1595 def _dh(fig, ax):1596 a = ax[0]; a.axhline(0, color="0.75", lw=0.6)1597 a.plot(tw, Zcw, color=GRAY, lw=1.2, label=r"$Z_\chi(t)$, the constituent $L(s,\chi)$")1598 a.plot(tw, Zfw, color=CRIMSON, lw=1.5, label=r"$Z_f(t)$, Davenport–Heilbronn")1599 a.axvline(toff, color=CRIMSON, lw=0.8, ls=":")1600 a.text(toff + 0.06, 2.05, "off-line pair at $t=85.70$", fontsize=7.5, color=CRIMSON, ha="left", va="top")1601 a.set_xlim(84, 87); a.set_ylim(-7.2, 2.4); a.set_xlabel("$t$"); a.set_ylabel("rotated on-line value"); a.legend(loc="lower right"); panel_letter(a, "(a)")1602 b = ax[1]1603 b.step(xs_s, ys_s, where="post", color=CRIMSON, lw=1.4, label=r"$N_f(t)$, strip count")1604 b.step(xs_on, ys_on, where="post", color=GRAY, lw=1.0, ls="--", label=r"$N_{\mathrm{crit}}(t)$, on-line count")1605 b.axvline(toff, color=CRIMSON, lw=0.8, ls=":")1606 b.set_xlim(0, T_HI); b.set_ylim(0, 47.5); b.set_xlabel("$t$"); b.set_ylabel("cumulative zero count"); b.legend(loc="upper left"); panel_letter(b, "(b)")1607 axins = b.inset_axes([0.56, 0.08, 0.41, 0.35])1608 axins.step(xs_s, ys_s, where="post", color=CRIMSON, lw=1.4); axins.step(xs_on, ys_on, where="post", color=GRAY, lw=1.0, ls="--")1609 axins.fill_between([toff, T_HI], [43, 43], [45, 45], color=CRIMSON, alpha=0.12, step="post")1610 axins.set_xlim(84, 87); axins.set_ylim(41, 46.6); axins.set_xticks([84, 85, 86, 87]); axins.set_yticks([42, 44, 46]); axins.tick_params(labelsize=6.5, length=2)1611 axins.text(84.12, 45.15, r"$N_f-N_{\mathrm{crit}}=2$", fontsize=7, color=CRIMSON, ha="left", va="bottom")1612 paper_figure("fig_dh", _dh, 2.55)16131614 # ---------------- D6: the strip (Thm. strip, row D8) and the finite negative control (row D9); 13 October 20261615 # D6a the Davenport–Heilbronn zero off the line and its Klein-four orbit: f(ρ) = 0 at ρ = 0.808517182457 + 85.6993484854 i1616 # (A8), f(1 − ρ) = 0 and f(ρ̄) = 0 (the functional equation and the reality of the coefficients), 0 < Re ρ < 1, Re ρ ≠ ½;1617 # for ζ the strip is a theorem (Mathlib, Thm. strip): ζ(s) ≠ 0 for Re s ≥ 1, checked here at four points as a sanity reading1618 mp.mp.dps = 201619 rho = mp.mpc("0.808517182457", "85.6993484854")1620 vals = [abs(pair_values(s)[1]) for s in (rho, 1 - rho, mp.conj(rho))]1621 zeta_right = [abs(mp.zeta(s)) for s in (mp.mpf(1) + 20j, mp.mpf(2), mp.mpc(1, 85.7), mp.mpc(1.5, 1000))]1622 ok = all(v < mp.mpf(10) ** -9 for v in vals) and 0 < rho.real < 1 and abs(rho.real - mp.mpf(1) / 2) > 0.3 and all(z > 0.1 for z in zeta_right)1623 # 20:D8 (p20080)1624 check("D6a", "the strip: the DH zero ρ = 0.8085 + 85.699 i has 0 < Re ρ < 1 off the line, with f(1 − ρ) = f(ρ̄) = 0 (the Klein-four orbit); ζ ≠ 0 on four points with Re s ≥ 1", ok,1625 f"|f| at ρ, 1−ρ, ρ̄ = {[float(v) for v in vals]}; |ζ| on Re s ≥ 1: min {float(min(zeta_right)):.3f}", kind="[approx]")1626 # D6b the DH prime vector: Λ_f from the divisor recursion is supported outside the prime powers, and a(n) is not multiplicative1627 Lf = dh_lambda(40)1628 supp = [n for n in range(2, 41) if abs(Lf[n]) > 1e-12]1629 pp = {n for n in range(2, 41) if len(prime_factors(n)) == 1}1630 outside = [n for n in supp if n not in pp]1631 kappa = float(KAPPA); a = lambda n: [0, 1, kappa, -kappa, -1][n % 5]1632 ok = outside[:8] == [6, 12, 14, 18, 21, 24, 26, 28] and abs(a(6) - a(2) * a(3)) > 1.0 and abs(a(6) - 1) < 1e-151633 check("D6b", "the DH prime vector Λ_f from the divisor recursion is supported outside the prime powers (6, 12, 14, 18, 21, 24, 26, 28, …), and a(n) is not multiplicative: a(6) = 1 against a(2)a(3) = −κ²", ok,1634 f"support outside prime powers below 40: {outside}; a(6) = {a(6):.3f}, a(2)a(3) = {a(2) * a(3):.3f}", kind="[approx]")1635 # D6c the finite negative control (Thm. negative_control13, row D9): over F_13(√2), P(X) = ((X − 2)² − 2)((X + 1)² − 2) has exactly the1636 # roots 2 ± η and −1 ± η, traces 4, 4, 11, 11 against the line's 1 (Tr z = 1 ⟺ Re z = 2⁻¹ = 7), the root set closed under conjugation and1637 # under X ↦ 1 − X: finitude, conjugation and the reflection force no root onto the line1638 p13, nu13 = 13, 21639 E = [(a_, b_) for a_ in range(p13) for b_ in range(p13)]1640 add = lambda z, w: ((z[0] + w[0]) % p13, (z[1] + w[1]) % p13)1641 sub = lambda z, w: ((z[0] - w[0]) % p13, (z[1] - w[1]) % p13)1642 sq = lambda z: mul(z, z, nu13, p13)1643 P = lambda z: mul(sub(sq(sub(z, (2, 0))), (2, 0)), sub(sq(add(z, (1, 0))), (2, 0)), nu13, p13)1644 roots = sorted(z for z in E if P(z) == (0, 0))1645 expected = sorted([(2, 1), (2, 12), (12, 1), (12, 12)])1646 traces = [trace(z, p13) for z in roots]1647 line = [z for z in E if trace(z, p13) == 1]1648 conj_closed = all(((z[0], (-z[1]) % p13) in roots) for z in roots)1649 half_turn_closed = all((((1 - z[0]) % p13, (-z[1]) % p13) in roots) for z in roots) # 1 − z = (1 − a) − bη1650 ok = roots == expected and traces == [4, 4, 11, 11] and all(z[0] == 7 for z in line) and len(line) == 13 and conj_closed and half_turn_closed \1651 and all((x * x) % p13 != nu13 for x in range(1, p13))1652 # 20:D9 (p20083)1653 check("D6c", "the finite negative control over F_13(√2): ((X−2)²−2)((X+1)²−2) has exactly the roots 2 ± η, −1 ± η, traces 4, 4, 11, 11 against the line's 1, the root set closed under conjugation and under X ↦ 1 − X; none on the line", ok,1654 f"roots {roots}, traces {traces}, the line Re z = 7 ({len(line)} points)")1655 # ---------------- D7: the cells of the cycle (Def./Prop. cells, row D5); 13 October 20261656 # on Ω = 13 and 233 the cycle's Ω cells [2πy, 2π(y+1)): N_y the zeros of ζ in the cell (the argument principle, mpmath nzeros), C_y1657 # the sign changes of Z on the cell (the Riemann–Siegel formula with the first remainder term, a grid refined where a change is close);1658 # C_y ≤ N_y on every cell; every cell clean (C_y = N_y) and the totals agree, N(2πΩ) = Σ N_y = 21, 1038; about log y zeros per cell1659 def Z_rs(t):1660 t = np.asarray(t, float); a = np.sqrt(t / (2 * np.pi)); N = np.floor(a).astype(int); pfrac = a - N1661 th = np.imag(loggamma(0.25 + 0.5j * t)) - 0.5 * t * np.log(np.pi)1662 main = np.zeros_like(t)1663 for n in range(1, int(N.max()) + 1):1664 m = N >= n1665 main[m] += 2 * np.cos(th[m] - t[m] * np.log(n)) / np.sqrt(n)1666 psi = np.cos(2 * np.pi * (pfrac * pfrac - pfrac - 1 / 16)) / np.cos(2 * np.pi * pfrac)1667 return main + (-1) ** (N - 1) * a ** -0.5 * psi1668 for Om in (13, 233):1669 counts = [int(mp.nzeros(2 * mp.pi * (y + 1))) - int(mp.nzeros(2 * mp.pi * y)) for y in range(Om)]1670 changes = []1671 for y in range(Om):1672 lo, hi = 2 * math.pi * y, 2 * math.pi * (y + 1)1673 ts = np.linspace(lo, hi, 400 if y > 0 else 400); zs = Z_rs(ts) if y > 0 else None1674 if y == 0: # the first cell holds no zero (the first ordinate is 14.13 > 2π)1675 changes.append(0); continue1676 sc = int(np.sum(zs[:-1] * zs[1:] < 0))1677 if sc < counts[y]: # refine the grid where a close pair may hide1678 ts = np.linspace(lo, hi, 4000); zs = Z_rs(ts); sc = int(np.sum(zs[:-1] * zs[1:] < 0))1679 changes.append(sc)1680 total = int(mp.nzeros(2 * mp.pi * Om))1681 ok = all(c <= n for c, n in zip(changes, counts)) and changes == counts and sum(counts) == total and total == {13: 21, 233: 1038}[Om]1682 tail = counts[-10:]1683 # 20:D5 (p20029)1684 check(f"D7{'a' if Om == 13 else 'b'}", f"the cells of the cycle on Ω = {Om}: N_y from the argument principle, C_y from the sign changes of Z on each cell; C_y ≤ N_y and C_y = N_y on all {Om} cells, Σ N_y = N(2πΩ) = {total}", ok,1685 f"zeros per cell: max {max(counts)}, last ten {tail} against log y ≈ {math.log(Om):.2f}; empty cells {sum(1 for c in counts if c == 0)}", kind="[approx]")16861687# ======================================================================================================================1688# block E — e_resonance.py (until 24 September 2026), its former docstring:1689# e_resonance.py — block E: resonance, antipode, horizon-scale resolution, square-root cancellation1690# =================================================================================================1691# The units-chart side of the construction: Ramanujan sums c_q(n) as the standing waves of the1692# additive meridian, the Möbius-weighted resonance R_L(n) = Σ_{q≤L} μ(q)/φ(q) c_q(n), its limit the1693# von Mangoldt weight (Theorem resonance), the energy at the antipode (Obs. antipode), the resolving1694# band at the Subject horizon (Obs. horizon), and the 1/√p flatness of every single multiplicative-1695# chart mode against the prime indicator (Prop. flat, Obs. flat). Paper-local predicates:1696#1697# E1a Thm. resonance (b) EXACT Λ = μ ∗ log on the divisor lattice: −Σ_{d|n} μ(d) log d = log ℓ on prime powers1698# ℓ^k, 0 otherwise — verified symbolically (sympy) for n ≤ 300, in floating point1699# for n ≤ 50001700# E1b Thm. resonance (a),(c) [approx] Hardy: R_L(n) → (φ(n)/n) Λ(n) (Cesàro mean to L = 4000, prime powers ≤ 13:1701# within 0.02, and nearer that limit than Λ(n)); the intertwiner (n/φ(n)) R → Λ1702# E1c Def. resonance / Fig. emergence [approx] Ψ_L(N) = Σ_{n≤N} (n/φ(n)) R_L(n) → ψ(N): the maximal deviation on1703# N ≤ 30 falls monotonically over L = 15, 60, 2401704# E2a Obs. antipode EXACT per-mode energy μ²(q)/φ(q) on squarefree q: 1, ½, ¼, ⅙ at q = 2, 3, 5, 7, with its1705# unique maximum among the nonconstant modes at the antipode q = 2 (all 2 ≤ q ≤ 2000; q = 1 ties)1706# E2b Obs. antipode [approx] the additive-transform band energy of the prime indicator on Z/10007, Parseval-1707# normalised over the bins within three of each a/q, peaks at the antipode and1708# follows 1/φ(q): 1153 > 592 > 301 > 205 at q = 2, 3, 5, 7 (the paper's figures,1709# restated on 2026-09-13 from the unnormalised 1220 > 640 > 323 > 214)1710# E3 Obs. horizon [approx] the resolving threshold L*(H) — the least bandwidth at which R_L separates every prime1711# power in {2..H} from every non-prime-power — exists for H = 6..50 and tracks the1712# horizon, L* < 6H, far below the field scale H²1713# E4 Obs. flat [approx] the maximal correlation of the prime indicator with a multiplicative-chart mode χ_j1714# is a small multiple of the 1/√(p−2) floor of Prop. flat: 0.077 against 0.0315 at1715# p = 1009 (2.4×); 0.0087 against 0.0032 over the modes of p = 100049 (2.7×); by1716# contrast the units-chart resonance correlates 0.90 with the primes in its window1717#1718# Figures: fig_obstruction.pdf (E3, E4), fig_emergence_frc.pdf (E1c beside the classical reconstruction of ψ from zeros).17191720# ------------------------------------------------------------------ arithmetic tables1721def sieve_mu_phi(M):1722 mu = np.ones(M + 1, dtype=int); phi = np.arange(M + 1); comp = np.zeros(M + 1, dtype=bool); primes = []1723 mu[0] = 01724 for i in range(2, M + 1):1725 if not comp[i]:1726 primes.append(i); mu[i] = -1; phi[i] = i - 11727 for q in primes:1728 if i * q > M: break1729 comp[i * q] = True1730 if i % q == 0:1731 mu[i * q] = 0; phi[i * q] = phi[i] * q; break1732 mu[i * q] = -mu[i]; phi[i * q] = phi[i] * (q - 1)1733 return mu, phi17341735def divisors(n):1736 return [d for d in range(1, n + 1) if n % d == 0]17371738def ramanujan(q, n, mu):1739 """Kluyver: c_q(n) = Σ_{d | gcd(q,n)} d μ(q/d), an integer."""1740 g = gcd(q, n)1741 return sum(d * mu[q // d] for d in divisors(g))17421743def resonance_table(Qmax, Nmax, mu, phi):1744 """cum[L, n] = R_L(n) = Σ_{q≤L} μ(q)/φ(q) c_q(n), L = 0..Qmax, n = 0..Nmax."""1745 C = np.zeros((Qmax + 1, Nmax + 1))1746 for q in range(1, Qmax + 1):1747 for n in range(1, Nmax + 1):1748 C[q, n] = ramanujan(q, n, mu)1749 w = np.array([0.0] + [mu[q] / phi[q] for q in range(1, Qmax + 1)])1750 return np.cumsum(w[:, None] * C, axis=0)17511752def von_mangoldt_array(N):1753 n, lam = von_mangoldt_support(N)1754 out = np.zeros(N + 1); out[n.astype(int)] = lam1755 return out17561757def dlog_table(p):1758 g = int(sp.primitive_root(p)); lam = np.zeros(p, dtype=np.int64); x = 11759 for k in range(p - 1):1760 lam[x] = k; x = (x * g) % p1761 return lam17621763def chart_mode_correlations(p):1764 """|corr(1_Π, χ_j)| for every nontrivial multiplicative-chart mode χ_j(n) = ω^{j λ_p(n)} on F_p^×, by FFT1765 over the discrete logarithm; returns (max, the exact root-mean-square 1/√(p−2) of Prop. flat, the array)."""1766 lam = dlog_table(p)1767 ind = np.zeros(p); ind[sieve_primes(p - 1)] = 1.01768 u = np.zeros(p - 1); u[lam[1:]] = ind[1:] # u[m] = 1_Π(g^m)1769 u = u - u.mean()1770 F = np.fft.fft(u)1771 corr = np.abs(F[1:]) / (np.linalg.norm(u) * math.sqrt(p - 1))1772 return float(corr.max()), 1 / math.sqrt(p - 2), corr17731774def pearson(a, b):1775 a = a - a.mean(); b = b - b.mean()1776 return float(a @ b / (np.linalg.norm(a) * np.linalg.norm(b)))17771778def block_E():1779 """Block E — the units-chart resonance and the von Mangoldt weight, the staircase from the modes, the antipode, the resolving band at the horizon, the √p-flatness of the chart modes; 7 checks (EXACT, [approx])."""1780 mp.mp.dps = 201781 print("\n== block E: resonance, antipode, horizon, square-root cancellation ==")1782 Qmax, Nmax = 320, 601783 mu, phi = sieve_mu_phi(max(Qmax, 5000))1784 with Timer("resonance table"):1785 cum = resonance_table(Qmax, Nmax, mu, phi)1786 Lam = von_mangoldt_array(5000)1787 # ---------------- E1a: Λ = μ * log, exact1788 sym_ok = True1789 for n in range(2, 301):1790 expr = sp.expand_log(-sum(sp.Integer(int(mu[d])) * sp.log(d) for d in divisors(n) if d > 1), force=True)1791 f = sp.factorint(n)1792 target = sp.log(list(f)[0]) if len(f) == 1 else sp.Integer(0)1793 if sp.simplify(expr - target) != 0:1794 sym_ok = False; break1795 flt = np.array([-sum(mu[d] * math.log(d) for d in divisors(n) if d > 1) for n in range(2, 5001)])1796 flt_ok = np.max(np.abs(flt - Lam[2:5001])) < 1e-91797 # illustration: 20:C2 retired 2026-10-15 (T11)1798 check("E1a", "Λ = μ ∗ log: −Σ_{d|n} μ(d) log d = Λ(n), symbolic for n ≤ 300, floating point for n ≤ 5000", sym_ok and flt_ok,1799 f"symbolic {'ok' if sym_ok else 'FAIL'}; max float deviation {np.max(np.abs(flt - Lam[2:5001])):.1e}", kind="EXACT")1800 # ---------------- E1b: Hardy's limit and the intertwiner (Cesàro mean over L ≤ 4000, conditional convergence)1801 Qhi = 4000 if not FAST else 15001802 with Timer(f"Cesàro resonance to L = {Qhi}"):1803 cumhi = resonance_table(Qhi, 14, mu, phi)1804 pps = [2, 3, 4, 5, 7, 8, 9, 11, 13]; comp = [6, 10, 12, 14]1805 R = {n: float(cumhi[1:Qhi + 1, n].mean()) for n in pps + comp}1806 hardy = all(abs(R[n] - phi[n] / n * Lam[n]) < 0.02 and abs(R[n] - phi[n] / n * Lam[n]) < abs(R[n] - Lam[n]) for n in pps)1807 zero_ok = all(abs(R[n]) < 0.02 for n in comp)1808 inter = all(abs(n / phi[n] * R[n] - Lam[n]) < 0.03 for n in pps)1809 check("E1b", "Hardy: R_L(n) → (φ(n)/n) Λ(n) on prime powers, → 0 off them; intertwiner (n/φ(n)) R → Λ = log ℓ", hardy and zero_ok and inter,1810 f"R(2) = {R[2]:.4f} vs ½ log 2 = {0.5 * Lam[2]:.4f}; R(9) = {R[9]:.4f} vs ⅔ log 3 = {2 / 3 * Lam[9]:.4f}; (13/12) R(13) = {13 / 12 * R[13]:.4f} vs log 13 = {Lam[13]:.4f}; max |R| off prime powers {max(abs(R[n]) for n in comp):.3f}", kind="[approx]")1811 # ---------------- E1c: the staircase from the modes1812 ns = np.arange(2, 31); psi = np.cumsum(Lam[2:31])1813 def Psi(L): return np.cumsum([n / phi[n] * cum[L, n] for n in ns])1814 devs = [float(np.max(np.abs(Psi(L) - psi))) for L in (15, 60, 240)]1815 # illustration: 20:C3 retired 2026-10-15 (T05)1816 check("E1c", "Ψ_L(N) → ψ(N) on N ≤ 30: maximal deviation falls monotonically over L = 15, 60, 240", devs[0] > devs[1] > devs[2],1817 "max |Ψ_L − ψ| = " + ", ".join(f"{d:.3f}" for d in devs), kind="[approx]")1818 # ---------------- E2a: the per-mode energy, exact1819 E = {q: Fraction(int(mu[q] ** 2), int(phi[q])) for q in range(2, 2001)}1820 ok = (E[2], E[3], E[5], E[7]) == (Fraction(1), Fraction(1, 2), Fraction(1, 4), Fraction(1, 6)) and all(E[q] < 1 for q in E if q != 2)1821 # illustration: 20:B11 retired 2026-10-15 (T05)1822 check("E2a", "per-mode energy μ²(q)/φ(q) = 1, ½, ¼, ⅙ at q = 2, 3, 5, 7; unique maximum among the nonconstant modes at the antipode q = 2 (2 ≤ q ≤ 2000)", ok,1823 "energies " + ", ".join(f"q={q}: {E[q]}" for q in (2, 3, 5, 7, 11)), kind="EXACT")1824 # ---------------- E2b: the additive-transform band energy (reported against the paper's figures, not pinned)1825 p4 = 100071826 ind = np.zeros(p4); ind[sieve_primes(p4 - 1)] = 1.01827 F = np.fft.fft(ind); norm2 = float(ind @ ind)1828 band = {}1829 for q in (2, 3, 5, 7):1830 e = 0.01831 for a in range(1, q):1832 if gcd(a, q) == 1:1833 kc = a * p4 / q1834 e += sum(abs(F[k % p4]) ** 2 for k in range(int(math.floor(kc)) - 3, int(math.ceil(kc)) + 4))1835 band[q] = e / norm21836 law = [band[q] * phi[q] / band[2] for q in (3, 5, 7)]1837 stated = {2: 1153, 3: 592, 5: 301, 7: 205}1838 ok = band[2] > band[3] > band[5] > band[7] and all(abs(band[q] - stated[q]) < 2 for q in stated) and all(0.85 < r < 1.25 for r in law)1839 check("E2b", "additive-transform band energy of the prime indicator on Z/10007 (Parseval-normalised, bins within 3 of a/q): 1153 > 592 > 301 > 205, the law 1/φ(q)", ok,1840 " > ".join(f"q={q}: {band[q]:.0f}" for q in (2, 3, 5, 7)) + f"; φ(q)·E_q/E_2 = " + ", ".join(f"{r:.2f}" for r in law), kind="[approx]")1841 # ---------------- E3: the resolving threshold1842 Hs = list(range(6, 51, 2)); Ls = []1843 for H in Hs:1844 pp = [n for n in range(2, H + 1) if Lam[n] > 0]; npp = [n for n in range(2, H + 1) if Lam[n] == 0]1845 Lstar = next((L for L in range(2, Qmax + 1) if min(cum[L, n] for n in pp) > max(abs(cum[L, n]) for n in npp)), None)1846 Ls.append(Lstar)1847 ok = all(L is not None for L in Ls) and all(L < 6 * H for L, H in zip(Ls, Hs)) and all(L < H * H for L, H in zip(Ls, Hs) if H >= 8)1848 ratio = [L / H for L, H in zip(Ls, Hs)]1849 # illustration: 20:C4 retired 2026-10-15 (T05)1850 check("E3", "resolving threshold L*(H) exists for H = 6..50 and tracks the horizon: L* < 6H ≪ H²", ok,1851 f"L*/H in [{min(ratio):.2f}, {max(ratio):.2f}]; L*(12) = {Ls[Hs.index(12)]}, L*(20) = {Ls[Hs.index(20)]}, L*(50) = {Ls[-1]}", kind="[approx]")1852 # ---------------- E4: square-root cancellation of every chart mode1853 with Timer("chart-mode correlations"):1854 m1, f1, corr1 = chart_mode_correlations(1009)1855 p5 = int(sp.nextprime(100_000))1856 while p5 % 4 != 1: p5 = int(sp.nextprime(p5))1857 m5, f5, _ = chart_mode_correlations(p5)1858 cum200 = resonance_table(80, 200, mu, phi)1859 rn = np.arange(2, 201); res_corr = pearson(np.array([1.0 if sp.isprime(int(n)) else 0.0 for n in rn]), np.array([cum200[80, n] for n in rn]))1860 ok = (abs(m1 - 0.077) < 0.005 and abs(f1 - 0.0315) < 0.0005 and abs(m1 / f1 - 2.43) < 0.031861 and abs(m5 - 0.0087) < 0.0005 and abs(f5 - 0.0032) < 0.0001 and abs(m5 / f5 - 2.75) < 0.03 and res_corr > 0.85)1862 # illustration: 20:C6 retired 2026-10-15 (T05)1863 check("E4", "max chart-mode correlation with the primes is a small multiple of the 1/√(p−2) floor: 0.077 vs 0.0315 (p = 1009, 2.4×), 0.0087 vs 0.0032 (p = 100049, 2.7×); resonance correlates 0.90", ok,1864 f"p = 1009: max {m1:.4f}, floor {f1:.4f}, ratio {m1 / f1:.3f}; p = {p5}: max {m5:.4f} over {(p5 - 1) // 2} mode pairs, floor {f5:.4f}, ratio {m5 / f5:.3f}; corr(1_Π, R_80) on n ≤ 200 = {res_corr:.2f}", kind="[approx]")1865 # ---------------- figures1866 FIN, TEAL = "#3b34a8", "#1f9e8a"1867 fig, ax = plt.subplots(1, 2, figsize=(12, 4.0))1868 ax[0].plot(np.arange(1, len(corr1) // 2 + 1), corr1[:len(corr1) // 2], color=TEAL, lw=0.6)1869 ax[0].axhline(f1, color="#888", ls="--", lw=1, label=r"floor $1/\sqrt{p-2}$")1870 ax[0].axhline(res_corr, color=FIN, lw=1.8, label=f"units-chart resonance ({res_corr:.2f})")1871 ax[0].set_ylim(0, 1); ax[0].set_xlabel("chart mode index $j$"); ax[0].set_ylabel(r"$|{\rm corr}(1_\Pi,\chi_j)|$")1872 ax[0].set_title(f"E4: chart modes are blind to primality (p = 1009, max {m1:.3f})", fontsize=10); ax[0].legend(fontsize=8, loc="center right")1873 ax[1].plot(Hs, Ls, "o-", color=FIN, lw=1.4, ms=4, label=r"resolving scale $L^*(H)$")1874 ax[1].plot(Hs, Hs, "--", color=TEAL, lw=1.2, label=r"horizon $H=\lfloor\sqrt{p}\rfloor$")1875 ax[1].plot(Hs, [h * h for h in Hs], ":", color="#b00", lw=1.5, label=r"field $p\sim H^2$")1876 ax[1].set_yscale("log"); ax[1].set_xlabel("window horizon $H$"); ax[1].set_ylabel("frequency scale")1877 ax[1].set_title("E3: the resolving band sits at the horizon", fontsize=10); ax[1].legend(fontsize=8, loc="upper left")1878 plt.tight_layout(); plt.savefig(f"{FIGDIR}/fig_obstruction.pdf", bbox_inches="tight"); plt.savefig(f"{FIGDIR}/fig_obstruction.png", dpi=110, bbox_inches="tight"); plt.close()1879 print(f" wrote {FIGDIR}/fig_obstruction.pdf")1880 # emergence: classical zeros (left) against units-chart modes (right)1881 g = riemann_zeros(80)1882 X0, X1 = 2, 301883 xs = np.linspace(X0 + 0.02, X1, 1400)1884 rho = [mp.mpc(0.5, float(gg)) for gg in g]1885 def psi_M(x, M):1886 s = sum(2 * mp.re(mp.power(x, r) / r) for r in rho[:M])1887 return float(x - mp.log(2 * mp.pi) - 0.5 * mp.log(1 - x ** -2) - s)1888 CL = ["#f0a8a8", "#df6b6b", "#c0392b"]; FR = ["#a8c4ec", "#5f8fd6", "#1f5fbf"]; GRIDC = "#d9d9d9"1889 fig, ax = plt.subplots(2, 2, figsize=(13.5, 8.2))1890 A = ax[0, 0]1891 for i, m in enumerate([1, 2, 3]):1892 A.plot(xs, [float(-2 * mp.re(mp.power(x, rho[m - 1]) / rho[m - 1])) for x in xs], color=CL[i], lw=1.4, label=fr"$\rho_{m}=\frac{{1}}{{2}}+i\,{float(g[m - 1]):.2f}$")1893 A.axhline(0, color=GRIDC, lw=0.8); A.set_title("classical spectral modes: the zeros"); A.set_xlabel("$x$"); A.set_ylabel(r"$-2\,{\rm Re}\,x^{\rho_m}/\rho_m$"); A.legend(fontsize=8); A.set_xlim(X0, X1)1894 B = ax[0, 1]1895 for i, q in enumerate([2, 3, 5]):1896 B.step(ns, [mu[q] / phi[q] * ramanujan(q, int(n), mu) for n in ns], where="mid", color=FR[i], lw=1.5, label=fr"$q={q}$ (energy {E[q]})")1897 B.axhline(0, color=GRIDC, lw=0.8); B.set_title(r"units-chart modes $\frac{\mu(q)}{\varphi(q)}c_q$ for $q=2,3,5$ (decreasing energy)"); B.set_xlabel("$n$"); B.legend(fontsize=8); B.set_xlim(X0, X1)1898 C = ax[1, 0]1899 C.step(np.concatenate([[X0], ns]), np.concatenate([[0], psi]), where="post", color="#111", lw=2.0, label=r"$\psi(x)$")1900 for i, M in enumerate([5, 20, 80]):1901 C.plot(xs, [psi_M(x, M) for x in xs], color=CL[i], lw=1.3, label=f"{M} zeros")1902 C.set_title(r"classical: $\psi(x)$ from $M$ zeros"); C.set_xlabel("$x$"); C.legend(fontsize=8, loc="upper left"); C.set_xlim(X0, X1); C.set_ylim(0, X1 + 4)1903 D = ax[1, 1]1904 D.step(np.concatenate([[X0], ns]), np.concatenate([[0], psi]), where="post", color="#111", lw=2.0, label=r"$\psi(N)$")1905 for i, L in enumerate([15, 60, 240]):1906 D.plot(ns, Psi(L), color=FR[i], lw=1.4, label=f"bandwidth $L={L}$ (max dev {devs[i]:.2f})")1907 D.set_title(r"E1c: $\Psi_L(N)=\sum_{n\leq N}\frac{n}{\varphi(n)}R_L(n)$ from $L$ units-chart modes"); D.set_xlabel("$N$"); D.legend(fontsize=8, loc="upper left"); D.set_xlim(X0, X1); D.set_ylim(0, X1 + 4)1908 plt.tight_layout(); plt.savefig(f"{FIGDIR}/fig_emergence_frc.pdf", bbox_inches="tight"); plt.savefig(f"{FIGDIR}/fig_emergence_frc.png", dpi=110, bbox_inches="tight"); plt.close()1909 print(f" wrote {FIGDIR}/fig_emergence_frc.pdf")191019111912# ======================================================================================================================1913# block T — Turing's certificate over the laboratory Carriers' cycles (15 October 2026, T03 of reports/blueprint-20261013.md)1914# =======================================================================================1915# For a Carrier Ω the cycle of heights is [0, 2πΩ). The certificate decides N(2πΩ) = C(2πΩ): every zero of ζ with ordinate1916# below 2πΩ is simple and on the line, cell by cell (Thm. hosted, row 20:F2). Method (Turing 1953; Lehman 1970; Brent 1979;1917# the constants of Trudgian 2011 for t ≥ 168π):1918# * Z(t) by the Riemann–Siegel formula with the first correction term and Gabcke's remainder bound |R_1| ≤ 0.053 τ^{-5/4}1919# (τ = t/2π, t ≥ 200), in double precision with an explicit rounding bound; a sign is certified only where |Z| exceeds1920# both bounds;1921# * Gram points g_n (θ(g_n) = nπ) by Newton from the Lambert-W initial value; g_n is good when (−1)^n Z(g_n) > 0 is1922# certified; a Gram block is the span between consecutive good Gram points; Rosser's rule for a block of k Gram1923# intervals: at least k certified sign changes of Z in it (a grid of m points per interval, refined where fewer are found);1924# * Turing's theorem at a good Gram point g_n: if the blocks after g_n covering h Gram intervals satisfy Rosser's rule and1925# g_{n+h} − g_n > 2.067 + 0.059 log(g_{n+h}/2π), then S(g_n) ≤ 0; with the same on the blocks before, S(g_n) = 0, that is1926# N(g_n) = n + 1; the certified sign changes below g_n then equal N(g_n) or fall short, and equality gives C = N: every1927# zero below g_n simple and on the line. The rule may fail in the interior (Lehman 1970; Brent 1979: the first exception1928# at g_13999825): a block short at the finest grid keeps its brackets, and its union with two blocks on each side is read1929# at the finest grid, the missing zeros lying in a neighbour; the count is pinned by the two ends, not by the interior;1930# * below the first good Gram point above 168π the count and the sign changes are the host's at twenty-five digits (mpmath).1931# The cells: N_y is the number of certified sign-change brackets in [2πy, 2π(y+1)); a bracket straddling a cell boundary is1932# bisected with certified signs until it lies in one cell.1933# T1 Thm. hosted [approx] Ω = 233: N(2πΩ) = 1038, every cell clean (≈ 1.5 min, the low range by mpmath)1934# T2 Thm. hosted [approx] Ω = 30089: N(2πΩ) = 280186, every cell clean (≈ 1.5 min)1935# T3 Thm. hosted [approx] Ω = 2408561: N(2πΩ) = 32984135, every cell clean — the embedded record CERT_2408561 of this1936# script's run of 15 October 2026 (RH_CERT_FULL=1 recomputes it, about an hour on one core)19371938TWO_PI = 2 * math.pi1939EPS = 2.0 ** -52194019411942def cert_gram_points(n):1943 """g_n with θ(g_n) = nπ, by Newton from the Lambert-W initial value (vectorised)."""1944 n = np.asarray(n, float)1945 x = (n + 7 / 8) / math.e1946 t = TWO_PI * (n + 7 / 8) / np.real(lambertw(x))1947 for _ in range(6):1948 t = t - (theta(t) - n * math.pi) / (0.5 * np.log(t / TWO_PI))1949 return t195019511952# Ψ(p) = cos(2π(p² − p − 1/16)) / cos(2πp) and its third derivative, with the removable singularities at p = 1/4, 3/41953_p = sp.symbols("p")1954_Psi = sp.cos(2 * sp.pi * (_p ** 2 - _p - sp.Rational(1, 16))) / sp.cos(2 * sp.pi * _p)1955_Psi3 = sp.diff(_Psi, _p, 3)1956_f0 = sp.lambdify(_p, _Psi, "numpy")1957_f3 = sp.lambdify(_p, _Psi3, "numpy")1958_ser = {}1959for _c in (sp.Rational(1, 4), sp.Rational(3, 4)):1960 _s0 = sp.series(_Psi, _p, _c, 10).removeO()1961 _s3 = sp.diff(_s0, _p, 3)1962 _ser[float(_c)] = (sp.lambdify(_p, _s0, "numpy"), sp.lambdify(_p, _s3, "numpy"))196319641965def cert_psi_c0_c1(p):1966 p = np.asarray(p, float)1967 c0 = np.empty_like(p); c3 = np.empty_like(p)1968 near = np.zeros(p.shape, bool)1969 for c, (f0s, f3s) in _ser.items():1970 m = np.abs(p - c) < 0.021971 c0[m] = f0s(p[m]); c3[m] = f3s(p[m]); near |= m1972 c0[~near] = _f0(p[~near]); c3[~near] = _f3(p[~near])1973 return c0, -c3 / (96 * math.pi ** 2)197419751976def cert_Z_certified(t):1977 """Z(t) with the first correction term, and the bound |error| ≤ Gabcke's R_1 + the rounding bound. t ≥ 200."""1978 t = np.asarray(t, float)1979 tau = t / TWO_PI1980 a = np.sqrt(tau)1981 N = np.floor(a).astype(np.int64)1982 p = a - N1983 th = theta(t)1984 main = np.zeros_like(t)1985 fp = np.zeros_like(t)1986 Nmax = int(N.max())1987 for n in range(1, Nmax + 1):1988 m = N >= n1989 if n == 1:1990 main += 2 * np.cos(th)1991 fp += 2 * (np.abs(th) * 4 * EPS + EPS)1992 else:1993 ln = math.log(n); inv = 1 / math.sqrt(n)1994 main[m] += 2 * np.cos(th[m] - t[m] * ln) * inv1995 fp[m] += 2 * ((np.abs(th[m]) + t[m] * ln) * 4 * EPS + EPS) * inv1996 c0, c1 = cert_psi_c0_c1(p)1997 corr = (-1.0) ** (N - 1) * tau ** -0.25 * (c0 + c1 * tau ** -0.5)1998 bound = 0.053 * tau ** -1.25 + fp + 1e-13 * (np.abs(main) + 1)1999 return main + corr, bound200020012002def cert_certified_signs(t):2003 z, b = cert_Z_certified(t)2004 s = np.zeros(len(t), np.int8)2005 s[z > b] = 1; s[z < -b] = -12006 return z, s200720082009def cert_brackets(ts, s):2010 """Brackets (t_i, t_j) of consecutive certified points with opposite signs; uncertified points skipped."""2011 idx = np.nonzero(s)[0]2012 if len(idx) < 2: return np.zeros((0, 2))2013 sc = s[idx]2014 ch = np.nonzero(sc[:-1] * sc[1:] < 0)[0]2015 return np.stack([ts[idx[ch]], ts[idx[ch + 1]]], axis=1)201620172018def certify(Om, log=print, grid=4, max_grid=4096, chunk=20000, window=2, cache=None):2019 """Turing's certificate over the cycle of the Carrier Ω: the record and the per-cell counts (block T)."""2020 T = TWO_PI * Om2021 t0 = time.time()2022 # --- the low range by the host at high precision (mpmath): below the first good Gram point above 168π2023 import mpmath as mp2024 mp.mp.dps = 252025 n_a = int(math.ceil((theta(168 * math.pi + 2) / math.pi)))2026 g = cert_gram_points(np.arange(n_a, n_a + 40))2027 z, s = cert_certified_signs(g)2028 good = s * ((-1) ** (np.arange(n_a, n_a + 40) % 2)) > 02029 i = int(np.argmax(good)); n_a += i; g_a = float(g[i])2030 N_low = int(mp.nzeros(g_a))2031 assert N_low == n_a + 1, (N_low, n_a + 1)2032 low_grid = np.linspace(14.0, g_a, 40 * N_low)2033 low_z = np.array([float(mp.siegelz(x)) for x in low_grid])2034 low_s = np.sign(low_z).astype(np.int8)2035 low_br = cert_brackets(low_grid, low_s)2036 assert len(low_br) == N_low, (len(low_br), N_low)2037 log(f" low range: g_{n_a} = {g_a:.3f}, N = {N_low} zeros, {len(low_br)} sign changes of Z (mpmath, 25 digits) [{time.time()-t0:.0f} s]")2038 # --- Gram points and their signs up to beyond T2039 n_b_guess = int(theta(T) / math.pi) + 602040 ns = np.arange(n_a, n_b_guess + 1)2041 if cache and os.path.exists(cache):2042 z = np.load(cache); gs, ss = z["gs"], z["ss"]; assert len(gs) == len(ns)2043 log(f" Gram points and signs read from {cache}")2044 else:2045 gs = np.empty(len(ns)); ss = np.zeros(len(ns), np.int8)2046 for c in range(0, len(ns), chunk):2047 gg = cert_gram_points(ns[c:c + chunk]); gs[c:c + chunk] = gg2048 ss[c:c + chunk] = cert_certified_signs(gg)[1]2049 if cache: np.savez(cache, gs=gs, ss=ss)2050 par = np.where(ns % 2 == 0, 1, -1).astype(np.int8)2051 good = (ss * par) > 02052 log(f" Gram points g_{n_a}..g_{ns[-1]}: {len(ns)}, good {int(good.sum())} ({100*good.mean():.1f} %), uncertified signs {int((ss == 0).sum())} [{time.time()-t0:.0f} s]")2053 gi = np.nonzero(good)[0]2054 # --- the blocks: between consecutive good Gram points; Rosser's rule by certified sign changes on a grid.2055 # Memory-bounded (15 October, second version): the blocks are arrays, the grids are built per batch of blocks of one2056 # length, the brackets are collected per batch; no per-block Python objects.2057 nblocks = len(gi) - 1; blen = np.diff(gi)2058 log(f" blocks: {nblocks}, longest {int(blen.max())} Gram intervals")2059 br_parts = [low_br]2060 cert_brackets_fn = cert_brackets2061 # blocks of length 1 with the Gram points' own certified signs: one change each (the two signs differ by goodness)2062 one = gi[np.nonzero(blen == 1)[0]]2063 br_parts.append(np.column_stack([gs[one], gs[one + 1]]))2064 del one2065 # longer blocks: grids, refined where fewer than k changes are found2066 pend = np.nonzero(blen > 1)[0] # block indices2067 viol = np.zeros(nblocks, bool)2068 m = grid2069 def signs_of(flat):2070 s = np.zeros(len(flat), np.int8)2071 for c in range(0, len(flat), chunk):2072 s[c:c + chunk] = cert_certified_signs(flat[c:c + chunk])[1]2073 return s2074 while len(pend) and m <= max_grid:2075 nxt = []2076 for k in np.unique(blen[pend]):2077 sel = pend[blen[pend] == k]; L = int(k) * m + 12078 nb_batch = max(1, min(len(sel), 2 * chunk // L))2079 for c in range(0, len(sel), nb_batch):2080 ids = sel[c:c + nb_batch]; nb = len(ids)2081 A = gs[gi[ids]]; Bv = gs[gi[ids] + k]2082 ts = A[:, None] + (Bv - A)[:, None] * np.linspace(0.0, 1.0, L)[None, :]2083 s = signs_of(ts.ravel()).reshape(nb, L)2084 last = np.maximum.accumulate(np.where(s != 0, np.arange(L)[None, :], 0), axis=1) # index of the last certified point at or before i2085 prev = last[:, :-1]; cur = s[:, 1:]2086 rows = np.arange(nb)[:, None]2087 change = (cur != 0) & (s[rows, prev] * cur < 0)2088 cnt = change.sum(axis=1)2089 ok = cnt >= k2090 keep = ok | (m * 4 > max_grid) # at the last grid the short blocks keep their brackets too (the exceptions)2091 if keep.any():2092 r, i = np.nonzero(change[keep])2093 tso = ts[keep]; pro = prev[keep]2094 br_parts.append(np.column_stack([tso[r, pro[r, i]], tso[r, i + 1]]))2095 if (~ok).any():2096 nxt.append(ids[~ok])2097 nxt = np.concatenate(nxt) if nxt else np.zeros(0, np.int64)2098 log(f" grid {m} per Gram interval: {len(pend)} blocks, {len(nxt)} short of Rosser's rule [{time.time()-t0:.0f} s]")2099 pend = nxt; m *= 42100 viol[pend] = True2101 violations = pend2102 # the exceptions to Rosser's rule (Lehman 1970, Brent 1979): a block short of its count at the finest grid has its missing2103 # zeros in a neighbouring block; the union of the block with `window` blocks on each side is read at the finest grid, its2104 # brackets replacing the coarse ones, and the union is resolved when its certified sign changes reach its Gram intervals.2105 # The count below the Turing point does not depend on the rule holding in the interior: it is pinned by the host's count2106 # at g_a and Turing's theorem at g_n; the rule drives the refinement, and the unions name where the grid looked harder.2107 unions = []2108 if len(violations):2109 log(f" blocks short of Rosser's rule at grid {max_grid}: {len(violations)}, Gram indices {[int(ns[gi[j]]) for j in violations[:10]]}")2110 spans = sorted((max(0, j - window), min(nblocks - 1, j + window)) for j in violations)2111 merged = []2112 for a, b in spans:2113 if merged and a <= merged[-1][1] + 1: merged[-1] = (merged[-1][0], max(merged[-1][1], b))2114 else: merged.append((a, b))2115 fine_parts = []2116 for (ja, jb) in merged:2117 a = gi[ja]; b = gi[jb + 1]; k = int(b - a); L = k * max_grid + 12118 ts = np.linspace(gs[a], gs[b], L); s = signs_of(ts)2119 br = cert_brackets_fn(ts, s)2120 resolved = len(br) >= k; host_pts = 02121 if not resolved and (s == 0).any():2122 # a pair of zeros closer than the double-precision certificate resolves (|Z| between them below the bound): the2123 # uncertified points of the union are read by the host at thirty digits, as the low range is, and the signs2124 # of the union are completed with them2125 mp.mp.dps = 302126 unc = np.nonzero(s == 0)[0]; host_pts = len(unc)2127 for i_ in unc:2128 zv = mp.siegelz(mp.mpf(float(ts[i_])))2129 s[i_] = 1 if zv > 0 else (-1 if zv < 0 else 0)2130 br = cert_brackets_fn(ts, s)2131 resolved = len(br) >= k2132 short = [int(ns[gi[j]]) for j in violations if ja <= j <= jb] # the Gram indices of the short blocks in the union2133 unions.append(dict(blocks=[int(ja), int(jb)], gram=[int(ns[a]), int(ns[b])], short=short, intervals=k, changes=int(len(br)), host_points=int(host_pts), resolved=bool(resolved)))2134 fine_parts.append(br)2135 log(f" union of blocks {ja}..{jb} (g_{ns[a]}..g_{ns[b]}, {k} Gram intervals): {len(br)} certified sign changes at grid {max_grid}"2136 + (f", {host_pts} uncertified points read by the host at thirty digits" if host_pts else "") + f" — {'resolved' if resolved else 'UNRESOLVED'}")2137 # the coarse brackets inside the unions are replaced by the fine ones2138 lo = np.array([gs[gi[ja]] for ja, jb in merged]); hi = np.array([gs[gi[jb + 1]] for ja, jb in merged])2139 def outside(br):2140 mid = 0.5 * (br[:, 0] + br[:, 1])2141 inside = np.zeros(len(br), bool)2142 for l_, h_ in zip(lo, hi): inside |= (mid > l_) & (mid < h_)2143 return br[~inside]2144 br_parts = [outside(br) for br in br_parts if len(br)] + fine_parts2145 # --- Turing's theorem at the end point: the last good Gram point g_n ≥ T with the blocks after it satisfying the bound2146 B = lambda t: 2.067 + 0.059 * math.log(t / TWO_PI)2147 n_end = None2148 for j in range(len(gi) - 1, -1, -1):2149 a = gi[j]2150 if gs[a] < T: break2151 ok = False; h = 02152 for jj in range(j, len(gi) - 1):2153 if viol[jj]: break2154 h = gi[jj + 1] - a2155 if gs[a + h] - gs[a] > B(gs[a + h]): ok = True; break2156 if ok: n_end = (j, a, h)2157 assert n_end is not None, "no certifiable end point in range"2158 j, a, h = n_end2159 n_cert = int(ns[a]); g_cert = float(gs[a])2160 log(f" Turing at g_{n_cert} = {g_cert:.3f} (≥ T = {T:.3f}): the blocks after it span {int(h)} Gram intervals, {gs[a+h]-g_cert:.3f} > {B(gs[a+h]):.3f}; N(g_{n_cert}) = {n_cert + 1}")2161 # the blocks before the point over the bound's span satisfy Rosser's rule (S(g_n) ≥ 0; the symmetric half of Turing's theorem)2162 before_ok = False; hb = 02163 for jj in range(j - 1, -1, -1):2164 if viol[jj]: break2165 hb = a - gi[jj]2166 if gs[a] - gs[a - hb] > B(gs[a]): before_ok = True; break2167 log(f" the blocks before g_{n_cert} span {int(hb)} Gram intervals, {gs[a]-gs[a-hb]:.3f} {'>' if before_ok else '<='} {B(gs[a]):.3f}: S(g_{n_cert}) {'≥ 0' if before_ok else 'not bounded below'}")2168 # --- the count2169 allbr = np.concatenate([br for br in br_parts if len(br)]); del br_parts2170 allbr = allbr[allbr[:, 1] <= g_cert + 1e-9]2171 C_cert = len(allbr)2172 log(f" certified sign changes below g_{n_cert}: {C_cert} against N = {n_cert + 1}: {'EQUAL' if C_cert == n_cert + 1 else 'DEFICIT ' + str(n_cert + 1 - C_cert)}")2173 ok_total = (C_cert == n_cert + 1) and before_ok and all(u["resolved"] for u in unions)2174 # --- the zeros below T, and the cells2175 work = allbr[allbr[:, 0] < T]; del allbr # brackets below T and those straddling T2176 def cell(x): return np.floor(x / TWO_PI).astype(np.int64)2177 # bisect brackets that straddle T or a cell boundary (certified signs), until each lies in one cell2178 need = cell(work[:, 0]) != cell(work[:, 1])2179 it = 02180 while need.any() and it < 60:2181 w = work[need]2182 mid = 0.5 * (w[:, 0] + w[:, 1])2183 smid = signs_of(mid); sl = signs_of(w[:, 0])2184 left = (smid * sl) < 0 # the change is in the left half2185 unc = smid == 02186 w2 = w.copy()2187 w2[left, 1] = mid[left]; w2[~left & ~unc, 0] = mid[~left & ~unc]2188 w2[unc, 0] = w[unc, 0] + 0.37 * (w[unc, 1] - w[unc, 0]) # an uncertified midpoint: shrink from an off-centre point2189 work[need] = w22190 need = cell(work[:, 0]) != cell(work[:, 1]); it += 12191 assert not need.any(), "a bracket could not be confined to one cell"2192 work = work[work[:, 1] <= T]2193 N_T = len(work)2194 cells = np.bincount(cell(work[:, 0]), minlength=Om)[:Om]2195 rec = dict(Omega=Om, T=T, n_low=n_a, g_low=g_a, N_low=N_low, n_end=n_cert, g_end=g_cert, N_end=n_cert + 1, C_end=C_cert,2196 gram_points=len(ns), good_fraction=float(good.mean()), uncertified=int((ss == 0).sum()), blocks=nblocks,2197 longest_block=int(blen.max()), rosser_violations=len(violations), rosser_unions=unions, turing_bound_intervals=int(h),2198 turing_before_intervals=int(hb),2199 N_T=N_T, cells_max=int(cells.max()), cells_empty=int((cells == 0).sum()), cells_last10=cells[-10:].tolist(),2200 cells_sum=int(cells.sum()), certified=bool(ok_total and N_T == cells.sum()), seconds=round(time.time() - t0, 1))2201 log(f" N(2πΩ) = {N_T}, every zero simple and on the line: {rec['certified']}; cells: max {rec['cells_max']}, empty {rec['cells_empty']}, last ten {rec['cells_last10']} [{rec['seconds']} s]")2202 return rec, cells220322042205# The record of the run of certify(2408561) of 15 October 2026 (the cloud container, one core, the memory-bounded block phase):2206# the sdist ships the script alone, so T3 reads this record unless RH_CERT_FULL=1 recomputes it (about an hour).2207CERT_2408561 = {2208 "Omega": 2408561,2209 "T": 15133435.08664577,2210 "n_low": 290,2211 "g_low": 530.5318655440009,2212 "N_low": 291,2213 "n_end": 32984135,2214 "g_end": 15133435.125860486,2215 "N_end": 32984136,2216 "C_end": 32984136,2217 "gram_points": 32983905,2218 "good_fraction": 0.8464814581536055,2219 "uncertified": 231,2220 "blocks": 27920263,2221 "longest_block": 7,2222 "rosser_violations": 4,2223 "rosser_unions": [2224 {2225 "blocks": [2226 11921435,2227 119214392228 ],2229 "gram": [2230 13999520,2231 139995292232 ],2233 "intervals": 9,2234 "changes": 9,2235 "host_points": 0,2236 "resolved": True,2237 "short": [2238 139995252239 ]2240 },2241 {2242 "blocks": [2243 20663429,2244 206634332245 ],2246 "gram": [2247 24360728,2248 243607362249 ],2250 "intervals": 8,2251 "changes": 8,2252 "host_points": 22,2253 "resolved": True,2254 "short": [2255 243607302256 ]2257 },2258 {2259 "blocks": [2260 26069384,2261 260693882262 ],2263 "gram": [2264 30783327,2265 307833332266 ],2267 "intervals": 6,2268 "changes": 6,2269 "host_points": 0,2270 "resolved": True,2271 "short": [2272 307833292273 ]2274 },2275 {2276 "blocks": [2277 26193457,2278 261934612279 ],2280 "gram": [2281 30930925,2282 309309322283 ],2284 "intervals": 7,2285 "changes": 7,2286 "host_points": 0,2287 "resolved": True,2288 "short": [2289 309309272290 ]2291 }2292 ],2293 "turing_bound_intervals": 7,2294 "turing_before_intervals": 7,2295 "N_T": 32984135,2296 "cells_max": 17,2297 "cells_empty": 2,2298 "cells_last10": [2299 15,2300 15,2301 14,2302 15,2303 14,2304 15,2305 14,2306 16,2307 14,2308 142309 ],2310 "cells_sum": 32984135,2311 "certified": True,2312 "seconds": 1916.62313}231423152316def block_T():2317 """Block T — Turing's certificate over the cycles of the laboratory Carriers 233, 30089 and 2408561 (Thm. hosted, row F2)."""2318 print("\n== block T: Turing's certificate over the laboratory cycles (Thm. hosted, row F2) ==")2319 here = os.path.dirname(os.path.abspath(__file__))2320 for pid, Om in (("T1", 233), ("T2", 30089), ("T3", 2408561)):2321 rec_path = os.path.join(here, f"turing_{Om}.json")2322 if Om == 2408561 and not os.environ.get("RH_CERT_FULL"):2323 rec = CERT_2408561; source = "the recorded run of 15 October 2026 (RH_CERT_FULL=1 recomputes it)"2324 else:2325 rec, cells = certify(Om, log=lambda s: print(" " + s))2326 try: json.dump(rec, open(rec_path, "w"), indent=1)2327 except OSError: pass2328 source = "computed now"2329 expected = {233: 1038, 30089: 280186, 2408561: 32984135}[Om]2330 unions = rec.get("rosser_unions", [])2331 ok = rec["certified"] and rec["N_T"] == expected and rec["C_end"] == rec["N_end"] and all(u["resolved"] for u in unions)2332 # 20:F2 (p20081)2333 check(pid, f"Turing's certificate over the cycle of Ω = {Om}: N(g_{rec['n_end']}) = {rec['N_end']} by Turing's theorem, {rec['C_end']} certified sign changes; N(2πΩ) = {rec['N_T']}, every cell clean", ok,2334 f"{source}; Gram points {rec['gram_points']}, good {100*rec['good_fraction']:.1f} %, uncertified {rec['uncertified']}, blocks {rec['blocks']}, longest {rec['longest_block']}, exceptions to Rosser's rule {rec['rosser_violations']} ({sum(u['resolved'] for u in unions)} unions resolved at the finest grid), Turing span {rec['turing_bound_intervals']} after / {rec.get('turing_before_intervals', '?')} before; cells max {rec['cells_max']}, empty {rec['cells_empty']}, last ten {rec['cells_last10']}; {rec['seconds']} s",2335 kind="[approx]")23362337if __name__ == "__main__":2338 want = [a.upper() for a in sys.argv[1:]] or sorted(BLOCK)2339 bad = [b for b in want if b not in BLOCK]2340 if bad: sys.exit(f"no block {', '.join(bad)}: the blocks are {', '.join(sorted(BLOCK))}")2341 t0 = time.time()2342 for b in want:2343 t = time.time(); _run_block(b); print(f" [block {b}: {time.time() - t:.0f} s]")2344 ok = summary(write=(want == sorted(BLOCK)))2345 kinds = {}2346 for r in RESULTS: kinds[r["kind"]] = kinds.get(r["kind"], 0) + 12347 print("by kind: " + ", ".join(f"{k} {v}" for k, v in sorted(kinds.items())) + f"; {len(RESULTS)} checks in {time.time() - t0:.0f} s" + ("; results.json written" if want == sorted(BLOCK) else ""))2348 sys.exit(0 if ok else 1)2349