src/13-epi — epi.py

1"""2epi.py — the validation package of "Finite Field Realisation of the Classical Constants π and e" (Akhtman, 2026; Preprints3doi 10.20944/preprints202607.1141.v1), the paper 13-epi of the FRC corpus (finite-ring-space/src/13-epi); one script since424 September 2026 (the registry epicommon.py and the six suites validate_e, validate_pi, validate_pi2, validate_towers,5kurepa_wall, frame_invariants merged; the two C passes stay beside it as kurepa_wall.c and frame_invariants.c, compiled on the6fly — a pure-python pass to a smaller bound runs without a compiler).7========================================================================================================================89One script, seven blocks, forty-seven family checks of exact micro-checks, standard library only. Every check is exact: integers,10residues, exact rationals (fractions.Fraction); the external targets e and π enter only as certified rational brackets of11the paper's own chains (the subfactorial chain for e, the Machin chain for π), and the binary64 constants enter as the12objects of study of the readout theorems (compared exactly with correctly rounded framed rationals). No floating-point13reference value of either constant decides any check; the few floating-point figures printed (a Pearson correlation, a14median, the Gauss sums) are [approx]-tagged diagnostics of the paper and decide nothing.1516A check of the registry is a *family* of micro-checks — one labelled claim of a block (e.R1, pi.W1, pi2.B2, tow.H,17kur.K1, …) — and names the predicate(s) of the paper's ledger it witnesses (LEDGER; predicates cited as 13:XN); the18deciding family's `# 13:XN (<key>)` marker is what the ledger's source column links (PREDICATES;19finitering.space/src/13-epi/#<key>). Master-ledger predicates reached through the paper predicates: 00:B820(13:J1, 13:J3), 00:B9 (13:J4, 13:J5), 00:C13 (13:J3), 00:C14 (13:B2, 13:H2), 00:B1 (13:B4), 00:C24 (13:G7 with A8, the21third-order wall) and 00:C26 (13:E5, the four chains as E- and G-partial sums); 13:Z1 and 13:Z2 are the paper's22Ω-hard certificates (block Z, the horizon: read over the shells p < √Ω of the totality, decided on the laboratory23Carriers over their own shells by block lab; 12 October 2026).2425    python3 epi.py                every block, results.json written; exit 1 if a family check fails (≈ 40 s, ≈ 30 s of it the C pass over 22 043 primes)26    python3 epi.py pi tow         the blocks named (e, pi, pi2, tow, kur, frm, lab); no results.json27    from frc_13_epi import predicate; predicate("13:F2")    one predicate: its blocks run once per session28"""29import os, json, sys, re, shutil, subprocess, time, math30from collections import OrderedDict31from fractions import Fraction32from fractions import Fraction as F33from fractions import Fraction as Fr34from decimal import Decimal, getcontext35from math import comb, gcd36import math as _m3738SCRIPT = os.path.splitext(os.path.basename(__file__))[0]        # "epi": the one script, the name results.json and the site pages carry39HERE = os.path.dirname(os.path.abspath(__file__))              # the C sources beside the script4041RESULTS = []42MICRO = []             # (script tag, family, label, ok)4344LEDGER = {45    # validate_e.py46    "e.R1": "13:D3",47    "e.R2": "13:D4, 13:I1",48    "e.R3": "13:D4",49    "e.R4": "13:D5, 13:E5",50    "e.P1": "13:C2",51    "e.T1": "13:D2",52    "e.W1": "13:F1, 13:F2",53    "e.W2": "13:F3",54    "e.W3": "13:F4",55    "e.W4": "13:F6",56    "e.W5": "13:F8, 13:Z1",57    "e.S1": "13:F5, 13:I1, 13:Z2",58    "e.A1": "13:H3, 13:H4, 13:I1",59    "e.N1": "13:H5, 13:I1, 13:Z2",60    # validate_pi.py61    "pi.R1": "13:E2, 13:C5",62    "pi.R2": "13:G8",63    "pi.R3": "13:E4, 13:E5, 13:I2",64    "pi.R4": "13:E3",65    "pi.W1": "13:G1, 13:G2, 13:I2",66    "pi.W2": "13:G2",67    "pi.S1": "13:G4",68    "pi.S2": "13:G4, 13:I2, 13:J6",69    "pi.S3": "13:G4, 13:G10",70    "pi.L1": "13:G5, 13:G8",71    "pi.V1": "13:G6",72    # validate_pi2.py73    "pi2.A1": "13:G3",74    "pi2.A2": "13:G3, 13:I2",75    "pi2.B1": "13:G7",76    "pi2.B2": "13:G7",77    "pi2.B3": "13:G7",78    "pi2.C1": "13:G6",79    "pi2.C2": "13:G6",80    "pi2.C3": "13:G6",81    "pi2.C4": "13:G6",82    "pi2.D1": "13:G7",83    "pi2.D2": "13:G5",84    # validate_towers.py85    "tow.E": "13:F7, 13:B3",86    "tow.P": "13:G9, 13:B3",87    "tow.C": "13:C3",88    "tow.O": "13:B2",89    "tow.H": "13:J2, 13:J3, 13:J5",90    "tow.W": "13:J1, 13:J2, 13:G2, 13:H2, 13:I4",91    "tow.D": "13:I3",92    # kurepa_wall.c93    "kur.K1": "13:F2, 13:J6, 13:Z1",94    # frame_invariants.c95    "frm.I1": "13:I5, 13:Z2",96    # the laboratory Carriers' shells (12 October 2026)97    "lab.Z1": "13:Z1, 13:Z2",98    "lab.Z2": "13:Z2",99}100101LABELS = {102    "e.R1": "superexponential enclosure of the derangement chain ε_n = n!/!n: |ε_n − e| < 8/(n+1)! and the alternating enclosure ε_even < e < ε_odd for 2 ≤ n ≤ 60, against the chain's own certified bracket",103    "e.R2": "binary64 readout: the IEEE-754 double nearest e is the correctly rounded ε_18 = 18!/!18 = 6402373705728000/2355301661033953, and n = 18 is minimal (ε_17 fails readout equality); the plateau thresholds below 200 bits, exact: ro_k(ε_n) within one ulp of ro_k(e) from n0(k), equal from n1(k) ≤ n0(k) + 1, with n1 = n0 + 1 exactly at k = 6, 21, 28, 113, 161, 177, 194, 198 (at k = 6, ε_5 = 30/11 reads out as 44/16, e as 43/16)",104    "e.R3": "100-digit determination: n = 70 is the least n with 8/(n+1)! < 10⁻¹⁰⁰, and ε_70 agrees with the certified value of e in its first 101 significant characters",105    "e.P1": "the derangement selector: the rencontres polynomial P_n(X) = Σ_{σ∈S_n} X^fix(σ) equals n! Σ_{k≤n} (X−1)^k/k! — the fixed-point census of S_n enumerated for n ≤ 8 (40 320 permutations at n = 8) against C(n,k)·D_{n−k} and against the expanded right side, the polynomial identity in exact integer coefficients for n ≤ 60 (P_3 = X³ + 3X + 2, P_4 = X⁴ + 6X² + 8X + 9), and the normalised kernel Z_n = P_n/n! with Z_n(1) = 1 and Z_n(0) = D_n/n!, so the unit registration ratio R_n = n!/D_n is the chain of e.R1",106    "e.T1": "the tail form of the derangement numbers: D_n = n!/e + (−1)^n δ_n with δ_n strictly in (1/(n+2), 1/(n+1)), so D_n is the nearest integer to n!/e, for n = 2 … 60 — decided with e bracketed by the chain itself at n = 109, 110 (width below 10⁻¹⁷⁰), each bound holding for the whole bracket",107    "e.R4": "feasibility selection: the exact step counts of the derangement chain for 2⁻ᵏ precision, k = 24, 53, 333, are n = 11, 18, 70 (the paper's table; the compound-chain column is the chart-side estimate e·2^{k−1})",108    "e.W1": "the wall of e on p ∈ {5, 13, 29, 101, 257, 1009, 10007}: Wilson reflection 1/k! ≡ −(−1)^k (p−1−k)! for every k < p; the wall identity !(p−1) ≡ K(p) with K(p) = Σ_{k<p} k! the left factorial; K(p) ≢ 0, so the terminal residue [ε_{p−1}] = −K(p)⁻¹ exists on every sample shell",109    "e.W2": "antiperiodicity !(n+p) ≡ −!n (mod p) over a full period n < p on the seven sample shells; the blind set Z₀(p) = {n < p : !n ≡ 0} always contains 1 and never p−1 there",110    "e.W3": "series duals at the wall for p ∈ {7, 13, 29, 101, 257, 1009}: Σ(−1)^k/k! ≡ −K(p), Σ 1/k! ≡ −A(p); the group law breaks, K(p)·A(p) ≡ 3, 0, 21, 28, 132, 258, never 1; the asymmetry at p = 13: A(13) ≡ 0, K(13) ≡ 10",111    "e.W4": "the p = 13 picture: the residue line [ε_n]₁₃, n = 2 … 12, is 2, 3, 7, 11, 1, 6, blind, 2, 8, 10, 9 (blind at n = 8, !8 = 13·1141), terminating on the wall invariant 9 = −K(13)⁻¹, the reading of 6/5 and of 11/7",112    "e.W5": "the tail identity j!·!(p−1−j) ≡ (−1)^j (K(p) − K(j)) for every j < p on the seven sample shells (j = 0 the wall identity), and the collision reading of the blind set: n ∈ Z₀(p) exactly when the partial sum K(p−1−n) collides with K(p); K(p−2) ≡ K(p) is the forced collision n = 1, and Kurepa's hypothesis at p is K(0) = 0 ≢ K(p)",113    "e.S1": "blind-scale statistics over the 501 primes in (1000, 5000): |Z₀(p)| distributes as 1¹⁶⁹ 2²⁰⁹ 3⁹² 4²² 5⁶ 6² 7¹, mean 1.996 (exact counts; the Poisson(1) reading of the excess is the paper's [approx] model, not a predicate)",114    "e.A1": "radian calibration over p < 10⁶: 39 175 shells p ≡ 1 (mod 4); exactly 260 of them satisfy |θ − 1| < 0.01 for one chirality, θ = 2π λ(i)/(p−1), decided with π as a Machin bracket of width < 10⁻¹⁰⁰; the best six (p, λ(i)) are (778933, 123966), (497773, 79227), (962237, 153133), (261017, 41537), (345181, 54945), (891997, 141986) in that order",115    "e.N1": "the null experiment over the 210 primes p ≡ 1 (mod 4), 5 < p < 3000: the index condition gcd(s, p−1) = gcd(λ(±i), p−1) linking the exponential unit to the wall residue −K(p)⁻¹ — the primitive root and the orientation of the quarter-turn varied independently — is solvable for exactly 94 primes against 96.0 ± 6.6 expected under independence; with the orientation tied to the generator (i = g^(−κ), Section 2) for exactly 54 against 57.9 ± 6.2; the Pearson correlation of the two angular addresses (≈ −0.10) is printed as the paper's [approx] reading",116    "pi.R1": "the Wallis pair v_n = 2·16ⁿ/((2n+1)C(2n,n)²), w_n = 16ⁿ/(n C(2n,n)²): v_n < π < w_n for n ≤ 300 against Machin brackets of width < 10⁻⁴²⁰, the exact width identity w_n − v_n = w_n/(2n+1), and strict monotonicity",117    "pi.R2": "the p = 13 picture: the residue line [w_n]₁₃, n = 1 … 6, is 4, 5, 10, 9, 6, 11",118    "pi.R3": "the Machin chain M_N: the IEEE-754 double nearest π is the correctly rounded M_10 and N = 10 is minimal; |M_N − π| < 17/((2N+3)5^{2N+3}) for 2 ≤ N ≤ 39; 100 decimal digits fixed at N = 71; the alternating pairing of the two arctangent series gives a certified two-sided enclosure at (12, 4) terms",119    "pi.R4": "the arcsin chain s_n = 3 Σ_{k≤n} C(2k,k)/((2k+1)16^k): the tail bound |π − s_n| < 4⁻ⁿ/(2n+3) for 1 ≤ n ≤ 120, against the Machin bracket",120    "pi.W1": "the legibility window on p ∈ {5, 13, 29, 101, 257, 1009, 10007}: p ∤ C(2n,n) for 1 ≤ n ≤ m = (p−1)/2 and p | C(2n,n) for m < n < p; the half-wall terminus [w_m] ≡ −2 = m⁻¹ = 4π_A; the lower chain blind exactly at the wall (2m+1 = p)",121    "pi.W2": "the wall residue −2 is universal: [w_m] ≡ −2 also for p ≡ 3 (mod 4), p ∈ {7, 23, 1031}",122    "pi.S1": "second order at the half wall for p ∈ {13, 29, 37, 41, 53, 101}: Morley's congruence (−1)^m C(p−1,m) ≡ 4^{p−1} (mod p³) and the formula w_m ≡ −2 + 2p(q_p(4) − 1) (mod p²)",123    "pi.S2": "the π-Wieferich search: 4^{p−1} ≡ 1 + p (mod p²) holds for exactly p ∈ {5, 45827} among the primes below 10⁶, and the recorded search of the standalone C programs (results_wieferich.txt) tiles [3, 2³⁵) with the same two hits (the normalised quotients' moments over p < 2·10⁴ are printed as the paper's [approx] reading)",124    "pi.S3": "the π-Wieferich condition in three forms on p ∈ {5, 13, 29, 37, 41, 53, 101, 45827}: q_p(4) = 2q_p(2) + p·q_p(2)² exactly, 4^{p−1} ≡ 1 + p (mod p²) iff q_p(4) ≡ 1 iff 2q_p(2) ≡ 1 (mod p); Eisenstein's 2q_p(2) ≡ −H_{(p−1)/2}; and the OEIS A355959 form (p+2)^{p−1} ≡ 1 (mod p²) iff 2q_p(2) ≡ 1, holding at 5 and 45827 only",125    "pi.L1": "Lucas revivals: [v_p] ≡ 8 for p ∈ {5, 7, 11, 13, 29, 37}; on p = 13, C(26,13) ≡ 2 = C(2,1)·C(0,0), C(28,14) ≡ 4 = C(2,1)² digitwise, and the two-digit revival scales number exactly (2κ)² = 36",126    "pi.V1": "first-order wall vanishing: σ_p = Σ_{k<m} C(2k,k)/((2k+1)16^k) ≡ 0 (mod p) for all 428 odd primes 5 ≤ p < 3000, with σ_3 = 1 the sole exception; the rational numerator of σ_p is divisible by p for 5 ≤ p < 200",127    "pi2.A1": "Gauss's central-binomial congruence C(2κ,κ) ≡ 2a* (mod p), a* ≡ 1 (mod 4) the odd part of p = a² + b², for all 211 primes p ≡ 1 (mod 4) below 3000",128    "pi2.A2": "the quarter-wall two-squares invariant [w_κ] = κ⁻¹(2a)⁻² = −(a²)⁻¹ = (b²)⁻¹ (mod p) for the same 211 primes; on p = 13 = 3² + 2², [w_3] = 10 = −(9)⁻¹",129    "pi2.B1": "Sun's supercongruence σ_p ≡ 0 (mod p²) re-verified in exact rational arithmetic for all primes 5 ≤ p < 300",130    "pi2.B2": "the third-order Bernoulli law σ_p/p² ≡ (−1)^{(p+1)/2} B_{p−3}/36 (mod p) for all sixty primes 5 ≤ p < 300, B_{p−3} from the exact Bernoulli recurrence",131    "pi2.B3": "the third-order Bernoulli law for all 1227 primes 5 ≤ p < 10⁴, σ_p modulo p³ by the recursion C(2k,k) = C(2k−2,k−1)·2(2k−1)/k (every divisor a p-unit) and B_{p−3} modulo p by the power sum Σ_{j<p} j^{p−3} ≡ p B_{p−3} (mod p²); the power-sum Bernoulli values agree with the exact recurrence on the sixty primes below 300",132    "pi2.C1": "the binomial transfer C(2k,k) ≡ (−4)^k C(m,k) (mod p) for k < m on p ∈ {11, 13, 29, 37}",133    "pi2.C2": "Lerch's congruence H_{(p−1)/2} ≡ −2 q_p(2) (mod p) for all primes 5 ≤ p < 500",134    "pi2.C3": "the Wallis evaluation L_m = Σ_j C(m,j)(−1)^j/(2j+1) = 4^m (m!)²/(2m+1)! over Q for m ≤ 60",135    "pi2.C4": "the key identity A_m + B_m = 2L_m over Q for m ≤ 60, with A = 2(F(½) − F(0)), B = 2(F(1) − F(½)), L = F(1) − F(0) for the formal antiderivative F of the expanded (1 − s²)^m — coefficient bookkeeping only",136    "pi2.D1": "the blind-range Euler congruence Σ_{m<k<p} C(2k,k)/((2k+1)16^k) ≡ p E_{p−3}/3 (mod p²) for all primes 5 ≤ p < 80",137    "pi2.D2": "the first revival to second order, v_p ≡ 8 + 16p(2q_p(2) − 1) (mod p²), for p ∈ {5, 7, 11, 13, 29, 37}",138    "tow.E": "the fixed-shell tower of e on p = 13 and 29, grades m ≤ 3: !(mp) ≡ (−1)^m, the shell reading [q_m] = e_p exactly, and q_m inside the certified window around e",139    "tow.P": "the fixed-shell tower of π on p = 13 and 29, n = p, p²: C(2n,n) ≡ 2, 16ⁿ ≡ 16, δ_p = ⟨−34⟩, the shell reading [π̂_{p,r}] = π_A exactly, and the r = 1 member inside the certified brackets",140    "tow.C": "the Cayley quarter-turn map C(x) = (1 + ix)/(1 − ix) on F_13 and F_29: C(0) = 1, C(1) = i, and the composition law C(x)C(y) = C((x+y)/(1−xy)) on every pair where the denominators are units",141    "tow.O": "orientation transport on F_13 and F_29, exhaustively over the units u of the exponent cycle: g' = g^u re-derives i' = i when u ≡ 1 (mod 4) and i' = −i when u ≡ 3 (mod 4)",142    "tow.H": "the height run over all 500 shells p ≡ 1 (mod 4), p ≤ 8009, in the smallest-primitive-root frame with H minimised over the two chiralities: π_A = [−1/2] pinned at height 2 on every shell; H(e_p) ≤ 2√p on every shell; exactly 68 shells with H ≤ 10; height 2 on exactly {13, 1933, 4177, 5857} (median H/√p ≈ 0.525 printed as the diagnostic)",143    "tow.D": "the duality on the shells p = 13 and 29: π angular-exact (g^{2κ} = −1) and wall-trivial (the terminus [w_m] = −2, the quarter invariant [w_κ] = −(a²)⁻¹, no blind scale below the wall), e angular-approximate (no exact calibration, decided on the Machin bracket) and wall-wild (the terminus K(p) ≠ 0 open at every shell, the blind set nonempty, the group law broken: K(p)A(p) ≠ 1)",144    "tow.W": "the wrap-free window and the calibration pin on p ∈ {13, 29, 101, 997, 8009}: ⌊√p⌋² < p and 2⌊√p⌋ < p; the pinning relation 2π_A + 1 ≡ 0; the half-turn tautology g^{2κ} = −1 for the smallest primitive root; the quarter-turn pin i² + 1 ≡ 0",145    "kur.K1": "Kurepa's wall in the derangement form for all 22 043 odd primes p < 2.5·10⁵ (kurepa_wall.c, one O(p) pass per prime in 128-bit modular arithmetic): !(p−1) ≡ K(p) (mod p) and K(p) ≢ 0 (mod p), zero failures of either condition",146    "frm.I1": "the three frame invariants of the two constants on the 4 783 shells p ≡ 1 (mod 4) below 10⁵ (frame_invariants.c): every triple (K(p) mod p, q_p(4) mod p, a + b i) exact — p = 13: (10, 6, −3+2i), p = 1093 the one shell with q_p(4) = 0 — the 4×4×4 contingency table of (K/p, q_p(4)/p, φ/π) and its chi-square 59.6 on 63 df exact, the means 0.502, 0.499 and variances 0.0828, 0.0831 by certified rational bounds; the independence reading is the paper's [approx] statement",147    "lab.Z1": "Kurepa's certificate on the laboratory Carriers Ω = 233, 30089, 2 408 561, over their own shells (the odd primes p with p² < Ω: 5, 39, 243 of them, the largest 13, 173, 1549): K(p) ≢ 0 (mod p) on every shell, Z = 0 on each Carrier, the walks costing Σ(p−1) = 34, 3 046, 172 452 steps; the Carriers' own residues as shells of the physical totality K(233) ≡ 69, K(30089) ≡ 5260, K(2408561) ≡ 435994; on the same shells the two range laws of Z2 (c): the blind-set sizes |Z₀(p)| (1⁸⁷ 2⁹³ 3⁴³ 4¹⁵ 5³ 6¹ 7¹ on the 243, mean 490/243 = 2.016 against Poisson(1)'s 2) and the index condition of H5 (solvable on 48 of the 116 shells p ≡ 1 (mod 4), p > 5, against 52.6 ± 4.9 under independence [approx])",148    "lab.Z2": "the calibration minimum over the laboratory Carriers' shells p ≡ 1 (mod 4) — 2, 18, 117 of them — decided exactly with π as a Machin bracket: δ = |2πλ/(p−1) − 1| is least at (p, λ) = (13, 5), (37, 6), (397, 63), δ = 1.618, 4.720·10⁻², 4.023·10⁻⁴, against H4's nearest-hit scale π/(2N) = 0.785, 0.0873, 0.0134 [approx]; the three frame invariants (K(p)/p, q_p(4)/p, φ/π) exact on the 117 shells of 2 408 561 (p = 397: 152, 383, −19+6i), their pairwise correlations +0.024, +0.067, +0.144 within 2σ = 0.185 of zero [approx]",149}150151BLOCK = {"e": "the derangement chain of e: the wall, the blind statistics, the calibration scan, the null experiment", "pi": "the chains of π: Wallis, Machin, arcsin; the window, the terminus, Morley, the Wieferich search, the revivals",      # check-id prefix -> the block (the function block_<name> below)152         "pi2": "the quarter invariant, Sun's supercongruence and the Bernoulli law, the proof ingredients", "tow": "the fixed-shell towers, the Cayley map, orientation transport, the height run, the pins",153         "kur": "the Kurepa wall over the 22 043 odd primes below 2.5e5 (a C pass)", "frm": "the frame invariants of the 4 783 shells p ≡ 1 (mod 4) below 1e5 (a C pass)",154         "lab": "the laboratory Carriers' own shells (p² < Ω for Ω = 233, 30089, 2 408 561): the certificates of block Z decided where they are decided"}155_FAM = [None, None]156157# the deciding family of each witnessed predicate: the one whose checks decide the predicate's statement (the other families158# that touch it are corroboration, listed by predicate() from the records)159PREDICATES = {160    "13:B2": "tow.O", "13:B3": "tow.E", "13:C2": "e.P1", "13:C3": "tow.C", "13:C5": "pi.R1", "13:D2": "e.T1", "13:D3": "e.R1", "13:D4": "e.R2",161    "13:D5": "e.R4", "13:E2": "pi.R1", "13:E3": "pi.R4", "13:E4": "pi.R3", "13:E5": "pi.R3", "13:F1": "e.W1",162    "13:F2": "e.W1", "13:F3": "e.W2", "13:F4": "e.W3", "13:F5": "e.S1", "13:F6": "e.W4", "13:F7": "tow.E",163    "13:F8": "e.W5", "13:G1": "pi.W1", "13:G2": "pi.W1", "13:G3": "pi2.A1", "13:G4": "pi.S1", "13:G5": "pi.L1",164    "13:G6": "pi.V1", "13:G7": "pi2.B2", "13:G8": "pi.R2", "13:G9": "tow.P", "13:G10": "pi.S3", "13:H2": "tow.W", "13:H3": "e.A1", "13:H4": "e.A1", "13:H5": "e.N1",165    "13:I1": "e.R2", "13:I2": "pi.R3", "13:I3": "tow.D", "13:I4": "tow.W", "13:I5": "frm.I1", "13:J1": "tow.W", "13:J2": "tow.H",166    "13:J3": "tow.H", "13:J5": "tow.H", "13:J6": "kur.K1", "13:Z1": "lab.Z1", "13:Z2": "lab.Z2",167}168_FAM = [None, None]169_RAN = set()                                            # blocks already run in this session (predicate() runs each once)170171def markers():172    """predicate label -> (script file, line) of its marker `# <paper>:<label> (<key>)`: the line of the check that decides it."""173    import re174    out = {}175    for i, line in enumerate(open(os.path.abspath(__file__), encoding="utf-8"), 1):176        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)177        if m:178            for lab in re.findall(r"\d+:[A-Z]+\d+[a-z]?", m.group(1)): out.setdefault(lab, (SCRIPT + ".py", i))179    return out180181def _run_block(name):182    if name not in _RAN:183        globals()[f"block_{name}"](); _RAN.add(name)184185def predicate(label, lines=14):186    """Verify one ledger predicate: run the blocks of the families that cite it (each once per session), print the deciding187    check's source (from its marker) and every family record that cites the predicate, and return True iff all pass."""188    fam = PREDICATES.get(label)189    if fam is None:190        print(f"{label}: no python witness (see the predicate's Lean witness or its source)"); return None191    citing = {f.split(".")[0] for f, rows in LEDGER.items() if label in [t.strip() for t in rows.split(",")]}192    for b in [b for b in BLOCK if b in {fam.split(".")[0]} | citing]: _run_block(b)     # the deciding block and every block whose families cite the predicate193    mk = markers().get(label)194    if mk:195        src = open(os.path.abspath(__file__), encoding="utf-8").read().split("\n")196        print(f"— {mk[0]}:{mk[1]} (the check that decides {label}; family {fam})")197        for j in range(mk[1] - 1, min(mk[1] - 1 + lines, len(src))): print(f"{j + 1:5d}  {src[j]}")198    else:199        print(f"— {SCRIPT}.py, block {fam.split('.')[0]} (the whole block witnesses {label})")200    recs = [r for r in RESULTS if label in [t.strip() for t in r["rows"].split(",")]]201    ok = all(r["ok"] for r in recs)202    for r in recs:203        role = "(deciding)" if r["id"] == fam else "(corroborating)"204        print(f"  [{'PASS' if r['ok'] else 'FAIL'}] {r['id']:8s} {role:16s} {r['detail']}")205    print(f"{label}: {'VERIFIED' if ok and recs else 'FAILED'} — {len(recs)} family record(s)")206    return ok207208def verify_all():209    """Run every block (those already run in this session are not re-run) and print the summary; True iff every family check passed."""210    for b in BLOCK: _run_block(b)                          # in the declared order211    return summary(write=False)212213def family(tag, fam):214    _FAM[0], _FAM[1] = tag, fam215216def chk(label, ok, tag=None, fam=None):217    """A micro-check, recorded under the current (tag, family); never exits."""218    MICRO.append((tag or _FAM[0], fam or _FAM[1], label, bool(ok)))219    if not ok:220        print(f"    FAIL: {label}")221222def flush(tag, order=None, details=None):223    fams = OrderedDict()224    for t, fam, label, ok in MICRO:225        if t == tag:226            fams.setdefault(fam, []).append(ok)227    for fam in (order or fams):228        oks = fams.get(fam, [])229        d = f"{len(oks)} exact micro-checks"230        if details and fam in details:231            d += "; " + details[fam]232        check(f"{tag}.{fam}", ok=bool(oks) and all(oks), detail=d)233234def check(pid, ok, detail="", kind="EXACT", label=None):235    ok = bool(ok)236    rows = LEDGER.get(pid, "")237    label = label or LABELS.get(pid, pid)238    RESULTS.append({"id": pid, "rows": rows, "block": pid.split(".")[0], "script": SCRIPT, "label": label, "ok": ok, "detail": detail, "kind": kind})239    print(f"  [{'PASS' if ok else 'FAIL'}] {pid:7s} {kind:5s} [{rows}] {label[:110]}" + (f"  --  {detail}" if detail else ""))240    return ok241242def summary(write=True):243    n_ok = sum(r["ok"] for r in RESULTS)244    print(f"\nSUMMARY: {n_ok}/{len(RESULTS)} checks passed ({len(MICRO)} exact micro-checks)" + ("" if n_ok == len(RESULTS) else "  <-- FAILURES"))245    if write:246        with open("results.json", "w") as f:247            json.dump(RESULTS, f, indent=1, ensure_ascii=False)248    return n_ok == len(RESULTS)249250# ------------------------------------------------------------------------------------------------------------251# block e — validate_e.py: FRC emergence of e: verification suite.252# Sections:253#   [1] Metric emergence of e from the derangement chain n!/D_n (exact rationals).254#   [2] Frame-internal arithmetic: antiperiodicity, blind scales, wall identity, series duals.255#   [3] Angular unit e_t = g^{i_t}: two-valued angle, radian calibration search over p < 10^6.256#   [4] Null experiment: pointwise identification g^{i_t} = k(e)-residue solvability statistics.257#258# Package form (2026-09): the memorandum script as written, its PASS/FAIL lines and the paper's stated figures turned259# into registry predicates (families e.R1-R4, e.W1-W5, e.S1, e.A1, e.N1). The radian-calibration scan of [3] now decides260# |theta - 1| < 0.01 and the ordering of the best shells in exact rational arithmetic with pi as a Machin bracket261# (width < 1e-100), as the paper states; the binary64 constant math.e is the object of study of the readout theorem262# and is compared exactly with the correctly rounded framed rationals. Every check is exact; the Pearson correlation263# of [4] is printed as the paper's [approx] reading and decides nothing.264OUT_e = []265def log_e(s=""):266    print(s); OUT_e.append(str(s))267268def pi_bracket(N=80):269    """Certified Machin bracket L < pi < U by alternating pairing; width < 1e-100 at N = 80."""270    def S(x, N):271        run, part = Fraction(0), []272        for k in range(N + 2):273            run += Fraction((-1)**k, (2*k + 1) * x**(2*k + 1)); part.append(run)274        return min(part[N], part[N+1]), max(part[N], part[N+1])275    l5, u5 = S(5, N); l239, u239 = S(239, N)276    return 16*l5 - 4*u239, 16*u5 - 4*l239277278# ---------- helpers ----------279def sieve_e(N):280    comp = bytearray(N + 1)281    for i in range(2, int(N**0.5) + 1):282        if not comp[i]:283            comp[i*i::i] = b'\x01' * len(comp[i*i::i])284    return [i for i in range(2, N + 1) if not comp[i]]285286def derangements_exact(N):287    D = [1, 0]288    for n in range(2, N + 1):289        D.append(n * D[n-1] + (-1)**n)290    return D291292def e_fraction(K=260):293    s, f = Fraction(0), 1294    for k in range(K):295        if k: f *= k296        s += Fraction(1, f)297    return s   # |s - e| < 3/K! < 10^-500 for K=260298299def frac_to_dec(fr, prec=140):300    getcontext().prec = prec301    return Decimal(fr.numerator) / Decimal(fr.denominator)302303_DIAG = [""]304305def block_e():306    """Block e — the derangement chain of e — enclosure, readouts, feasibility, the Kurepa wall, antiperiodicity, duals, the p = 13 line, blind statistics, radian calibration, the null experiment (e.R1 … e.N1)."""307    global E, D, fact, NMAX308    # =========================================================309    log_e("="*72); log_e("[1] METRIC EMERGENCE: the chain eps_n = n!/D_n")310    log_e("="*72)311    E = e_fraction(300)312    NMAX = 110313    D = derangements_exact(NMAX)314    fact = [1]315    for n in range(1, NMAX + 1): fact.append(fact[-1] * n)316317    # (a) error bound |n!/D_n - e| < 8/(n+1)!  and alternating enclosure318    bound_ok, encl_ok = True, True319    for n in range(2, 61):320        err = Fraction(fact[n], D[n]) - E321        if abs(err) >= Fraction(8, fact[n+1]): bound_ok = False322        if (err > 0) != (n % 2 == 1): encl_ok = False   # sign(eps_n - e) = (-1)^{n+1}323    log_e(f"(a) |n!/D_n - e| < 8/(n+1)! for n=2..60 : {'PASS' if bound_ok else 'FAIL'}")324    log_e(f"    alternating enclosure eps_even < e < eps_odd : {'PASS' if encl_ok else 'FAIL'}")325    # 13:D3 (p13019)326    family("e", "R1")327    chk("|n!/D_n - e| < 8/(n+1)! for n = 2..60", bound_ok)328    chk("alternating enclosure eps_even < e < eps_odd for n = 2..60", encl_ok)329330    # (b) minimal n whose framed rational rounds to the IEEE-754 double nearest e331    target = math.e332    n_double = None333    for n in range(2, 40):334        if float(Fraction(fact[n], D[n])) == target:335            n_double = n; break336    log_e(f"(b) minimal n with float(n!/D_n) == IEEE-double(e): n = {n_double}")337    log_e(f"    n!/D_n at n={n_double}: {fact[n_double]}/{D[n_double]}")338    log_e(f"    |eps_17 - e| = {float(abs(Fraction(fact[17],D[17])-E)):.3e},  half-ulp(e) = {math.ulp(math.e)/2:.3e}")339    # 13:D4 (p13020), 13:I1 (p13048)340    family("e", "R2")341    chk("minimal n with float(n!/D_n) == binary64(e) is n = 18", n_double == 18)342    chk("eps_18 = 6402373705728000/2355301661033953", Fraction(fact[18], D[18]) == Fraction(6402373705728000, 2355301661033953))343    chk("eps_17 fails readout equality", float(Fraction(fact[17], D[17])) != target)344    chk("|eps_17 - e| exceeds the half-ulp 2^-52 of e and |eps_18 - e| lies below it", abs(Fraction(fact[17], D[17]) - E) > Fraction(1, 2**52) > abs(Fraction(fact[18], D[18]) - E))345346    # (b') the plateau thresholds below 200 bits, exact: n0(k) from the enclosure bound e(e+1)/(n+1)! < 2^-k; the readout347    # ro_k on [2, 4) (grid 2^(2-k), ties to even); d_k(e) the distance from e to the nearest rounding boundary; the enclosure348    # certifies the plateau from the least N with e(e+1)/(N+1)! < d_k(e), the readouts below N compared exactly; n1(k) the least349    # n from which they agree. E is e to 300 terms, |E - e| < 1/301!, below every d_k(e) met (certified per k).350    def ro_bits(x, k):351        step = Fraction(1, 2**(k-2)); q = x / step; m = q.numerator // q.denominator; r = q - m352        if r > Fraction(1, 2) or (r == Fraction(1, 2) and m % 2 == 1): m += 1353        return m * step354    ee = E * (E + 1); tail = Fraction(1, math.factorial(301))355    ok_ulp = ok_cert = ok_bound = True; n1_exc = {}; first_n0 = None356    for k in range(2, 200):357        n0 = 2358        while ee / fact[n0+1] >= Fraction(1, 2**k): n0 += 1359        step = Fraction(1, 2**(k-2)); tq = E / step; frac = tq - (tq.numerator // tq.denominator)360        d = abs(frac - Fraction(1, 2)) * step            # distance from E to the nearest boundary (m + 1/2)*step361        if d <= tail: ok_cert = False362        N = 2363        while ee / fact[N+1] >= d - tail: N += 1364        rE = ro_bits(E, k)365        if abs(ro_bits(Fraction(fact[n0], D[n0]), k) - rE) > step: ok_ulp = False366        agree = [ro_bits(Fraction(fact[n], D[n]), k) == rE for n in range(n0, N + 1)]367        n1 = N368        for i in range(len(agree) - 1, -1, -1):369            if agree[i]: n1 = n0 + i370            else: break371        if n1 > n0 + 1: ok_bound = False372        if n1 == n0 + 1: n1_exc[k] = n0373        if k == 6: first_n0 = (n0, ro_bits(Fraction(fact[n0], D[n0]), k), rE)374    log_e(f"(b') plateau thresholds k < 200: within one ulp at n0 {'PASS' if ok_ulp else 'FAIL'}; n1 <= n0 + 1 {'PASS' if ok_bound else 'FAIL'}; n1 = n0 + 1 at k = {sorted(n1_exc)}")375    chk("ro_k(eps_n0) lies within one ulp of ro_k(e) for every 2 <= k < 200", ok_ulp)376    chk("the plateau threshold n1(k) <= n0(k) + 1 for every 2 <= k < 200, the comparison exact below the certified index", ok_bound)377    chk("n1 = n0 + 1 exactly at k = 6, 21, 28, 113, 161, 177, 194, 198", sorted(n1_exc) == [6, 21, 28, 113, 161, 177, 194, 198])378    chk("k = 6: n0 = 5, eps_5 = 30/11 reads out as 44/16, e as 43/16", first_n0 == (5, Fraction(44, 16), Fraction(43, 16)) and Fraction(fact[5], D[5]) == Fraction(30, 11))379    chk("the 300-term e is within every rounding cell met: d_k(e) > 1/301! for 2 <= k < 200", ok_cert)380381    # (c) 100-digit determination382    n100 = None383    for n in range(2, NMAX):384        if Fraction(8, fact[n+1]) < Fraction(1, 10**100):385            n100 = n; break386    d_chain = frac_to_dec(Fraction(fact[n100], D[n100]), 130)387    d_e     = frac_to_dec(E, 130)388    match = 0389    s1, s2 = str(d_chain), str(d_e)390    for a, b in zip(s1, s2):391        if a == b: match += 1392        else: break393    log_e(f"(c) 100-digit precision reached at n = {n100}; leading agreement with e: {match-1} chars")394    log_e(f"    e            = {str(d_e)[:105]}")395    log_e(f"    {n100}!/D_{n100}     = {str(d_chain)[:105]}")396    family("e", "R3")397    chk("least n with 8/(n+1)! < 10^-100 is n = 70", n100 == 70)398    chk("70!/D_70 agrees with e in at least 101 leading characters (100 significant digits)", match - 1 >= 101)399400    # (d) feasibility comparison with (1+1/n)^n401    log_e("(d) precision-cost comparison (steps n needed for 2^-k):")402    # 13:D5 (p13021)403    family("e", "R4")404    for k in (24, 53, 333):405        # derangement chain: 8/(n+1)! < 2^-k406        n1 = 2407        while Fraction(8, fact[n1+1]) >= Fraction(1, 2**k): n1 += 1408        # compound chain: |(1+1/n)^n - e| ~ e/(2n) -> n ~ e*2^{k-1}409        n2 = math.e * 2**(k-1)410        log_e(f"    k={k:3d}:  derangement chain n = {n1:3d}   |   (1+1/n)^n chain n ~ {n2:.2e}")411        chk(f"derangement chain reaches 2^-{k} at n = {(24, 53, 333).index(k) and [11, 18, 70][(24, 53, 333).index(k)] or 11}", n1 == [11, 18, 70][(24, 53, 333).index(k)])412413    # (e) the derangement selector (13:C2, 12 October 2026): the rencontres polynomial P_n(X) = sum_{sigma in S_n} X^fix(sigma)414    # equals n! sum_{k<=n} (X-1)^k/k! — the census of S_n by enumeration for n <= 8 (40 320 permutations at n = 8), the415    # polynomial identity sum_k C(n,k) D_{n-k} X^k = n! sum_j (X-1)^j/j! in exact integer coefficients for n <= 60, and the416    # normalised kernel Z_n = P_n/n! with Z_n(1) = 1, Z_n(0) = D_n/n!, so the unit registration ratio R_n = n!/D_n is the chain of (a)417    from itertools import permutations418    def rencontres_coeffs(n):                        # coefficients of n! sum_{j<=n} (X-1)^j/j! in powers of X, exact integers419        c = [Fraction(0)] * (n + 1)420        for j in range(n + 1):421            w = Fraction(fact[n], fact[j])           # n!/j!, an integer422            for i in range(j + 1):                   # (X-1)^j = sum_i C(j,i) X^i (-1)^(j-i)423                c[i] += w * comb(j, i) * (-1) ** (j - i)424        assert all(x.denominator == 1 for x in c)425        return [int(x) for x in c]426    census_ok = True427    for n in range(0, 9):428        census = [0] * (n + 1)429        for sg in permutations(range(n)):430            census[sum(1 for i in range(n) if sg[i] == i)] += 1431        if census != [comb(n, k) * D[n - k] for k in range(n + 1)] or census != rencontres_coeffs(n): census_ok = False432    identity_ok = all([comb(n, k) * D[n - k] for k in range(n + 1)] == rencontres_coeffs(n) for n in range(0, 61))433    kernel_ok = True434    for n in range(2, 61):435        c = rencontres_coeffs(n); Z1, Z0 = Fraction(sum(c), fact[n]), Fraction(c[0], fact[n])      # Z_n(1), Z_n(0)436        if not (Z1 == 1 and Z0 == Fraction(D[n], fact[n]) and Z1 / Z0 == Fraction(fact[n], D[n])): kernel_ok = False437    log_e(f"(e) rencontres census of S_n (n <= 8) against C(n,k) D_(n-k) and n! sum (X-1)^j/j!: {'PASS' if census_ok else 'FAIL'}; P_3(X) = {rencontres_coeffs(3)[::-1]} (descending)")438    # 13:C2 (p13013)439    family("e", "P1")440    chk("the fixed-point census of S_n for n <= 8 (enumerated) is C(n,k) D_(n-k) and equals the coefficients of n! sum_(j<=n) (X-1)^j/j!", census_ok)441    chk("sum_k C(n,k) D_(n-k) X^k = n! sum_(j<=n) (X-1)^j/j! as an exact polynomial identity for n <= 60", identity_ok)442    chk("P_3(X) = X^3 + 3X + 2 and P_4(X) = X^4 + 6X^2 + 8X + 9", rencontres_coeffs(3) == [2, 3, 0, 1] and rencontres_coeffs(4) == [9, 8, 6, 0, 1])443    chk("Z_n(1) = 1 and Z_n(0) = D_n/n!, so R_n = Z_n(1)/Z_n(0) = n!/D_n is the chain eps_n of (a), n = 2..60", kernel_ok)444445    # (f) the tail form (13:D2, 12 October 2026): D_n = n!/e + (-1)^n delta_n with delta_n in (1/(n+2), 1/(n+1)), so D_n is446    # the nearest integer to n!/e — decided for n = 2..60 with e bracketed by the chain itself (E_lo < e < E_hi from the447    # enclosure of (a) at n = 109, 110, width < 10^-170), each bound of delta_n holding for the whole bracket448    E_lo, E_hi = min(Fraction(fact[109], D[109]), Fraction(fact[110], D[110])), max(Fraction(fact[109], D[109]), Fraction(fact[110], D[110]))449    tail_ok, nearest_ok = True, True450    for n in range(2, 61):451        lo, hi = Fraction(fact[n]) / E_hi, Fraction(fact[n]) / E_lo          # the bracket of n!/e452        d_lo, d_hi = (-1) ** n * (D[n] - hi), (-1) ** n * (D[n] - lo)       # the bracket of delta_n = (-1)^n (D_n - n!/e)453        d_lo, d_hi = min(d_lo, d_hi), max(d_lo, d_hi)454        if not (Fraction(1, n + 2) < d_lo and d_hi < Fraction(1, n + 1)): tail_ok = False455        if not (abs(D[n] - lo) < Fraction(1, 2) and abs(D[n] - hi) < Fraction(1, 2)): nearest_ok = False456    log_e(f"(f) tail form D_n = n!/e + (-1)^n delta_n, delta_n in (1/(n+2), 1/(n+1)), n = 2..60: {'PASS' if tail_ok else 'FAIL'}; D_n the nearest integer to n!/e: {'PASS' if nearest_ok else 'FAIL'}")457    # 13:D2 (p13018)458    family("e", "T1")459    chk("the chain brackets e at n = 109, 110 with width below 10^-170", 0 < E_hi - E_lo < Fraction(1, 10 ** 170) and E_lo < E < E_hi)460    chk("delta_n = (-1)^n (D_n - n!/e) lies strictly in (1/(n+2), 1/(n+1)) for n = 2..60, for the whole bracket of e", tail_ok)461    chk("D_n is the nearest integer to n!/e for n = 2..60 (|D_n - n!/e| < 1/2 on the whole bracket)", nearest_ok)462    chk("the signed tail alternates: D_n - n!/e has the sign (-1)^n for n = 2..60", all(((D[n] - Fraction(fact[n]) / E) > 0) == (n % 2 == 0) for n in range(2, 61)))463464    # =========================================================465    log_e(); log_e("="*72); log_e("[2] FRAME-INTERNAL ARITHMETIC (per-prime structure)")466    log_e("="*72)467468    def d_mod(p, upto):469        """derangement residues D_0..D_upto mod p via recurrence"""470        d = [1 % p, 0]471        s = -1472        for n in range(2, upto + 1):473            s = -s474            d.append((n * d[n-1] + s) % p)475        return d476477    def kurepa_mod(p):478        f, s = 1, 1479        for k in range(1, p):480            f = f * k % p481            s = (s + f) % p482        return s483484    def alt_fact_mod(p):485        f, s, sg = 1, 1, 1486        for k in range(1, p):487            f = f * k % p488            sg = -sg489            s = (s + sg * f) % p490        return s % p491492    sample_primes = [5, 13, 29, 101, 257, 1009, 10007]493    log_e(f"{'p':>6} | D_(p-1)==!p | antiperiod D_(n+p)==-D_n | wall residue k(e)= -(!p)^-1 | blind set Z0(p)")494    for p in sample_primes:495        d = d_mod(p, 2*p)496        kp = kurepa_mod(p)497        wall = (d[p-1] == kp)498        anti = all(d[n+p] == (-d[n]) % p for n in range(0, p))499        z0 = [n for n in range(p) if d[n] == 0]500        ke = (-pow(kp, -1, p)) % p if kp else None501        # 13:F1 (p13027), 13:F2 (p13028)502        family("e", "W1")503        f_ = [1]504        for k in range(1, p): f_.append(f_[-1] * k % p)505        chk(f"p={p}: Wilson reflection 1/k! == -(-1)^k (p-1-k)! for all k < p", all(pow(f_[k], -1, p) == (-(-1)**k * f_[p-1-k]) % p for k in range(p)))506        chk(f"p={p}: wall identity D_(p-1) == K(p)", wall)507        chk(f"p={p}: K(p) != 0, the terminal residue -(K(p))^-1 exists", kp != 0)508        # 13:F3 (p13029)509        family("e", "W2")510        chk(f"p={p}: antiperiodicity D_(n+p) == -D_n over a full period", anti)511        chk(f"p={p}: 1 in Z0(p) and p-1 not in Z0(p)", 1 in z0 and (p-1) not in z0)512        # the tail identity and the collision reading of the blind set (13:F8; the Y1 push of 2026-09-18)513        # 13:F8 (p13066)514        family("e", "W5")515        Kj = [0]516        for k in range(p): Kj.append((Kj[-1] + f_[k]) % p)          # Kj[j] = K(j) = sum_{k<j} k! mod p517        chk(f"p={p}: tail identity j! D_(p-1-j) == (-1)^j (K(p) - K(j)) for all j < p", all(f_[j] * d[p-1-j] % p == ((-1)**j * (Kj[p] - Kj[j])) % p for j in range(p)))518        chk(f"p={p}: blind set = collisions of the partial sums: Z0(p) = {{p-1-j : K(j) == K(p), j <= p-2}}", set(z0) == {p-1-j for j in range(p-1) if Kj[j] == Kj[p]})519        chk(f"p={p}: the forced collision K(p-2) == K(p) (n = 1) and K(0) = 0 != K(p) (Kurepa at p)", Kj[p-2] == Kj[p] and Kj[p] != 0)520        log_e(f"{p:6d} |   {'PASS' if wall else 'FAIL'}      |          {'PASS' if anti else 'FAIL'}           |    {ke!s:>8}                | {z0 if len(z0)<=8 else str(z0[:8])+'...'} (|Z0|={len(z0)})")521522    # series duals523    log_e()524    log_e("Series duals at the wall:  S- := sum (-1)^k/k! == -!p ;  S+ := sum 1/k! == -A(p)")525    dual_ok = True526    prod_examples = []527    for p in [7, 13, 29, 101, 257, 1009]:528        invf, f = [1], 1529        for k in range(1, p):530            f = f * k % p531            invf.append(pow(f, -1, p))532        Sm = sum(((-1)**k) * invf[k] for k in range(p)) % p533        Sp = sum(invf[k] for k in range(p)) % p534        kp, ap = kurepa_mod(p), alt_fact_mod(p)535        if Sm != (-kp) % p or Sp != (-ap) % p: dual_ok = False536        prod_examples.append((p, (kp * ap) % p))537    log_e(f"  duals verified for p in [7,13,29,101,257,1009]: {'PASS' if dual_ok else 'FAIL'}")538    log_e(f"  group-law breaking at the wall: (!p * A(p)) mod p = {prod_examples}  (would be == 1 if exp(1)exp(-1)=1 survived)")539    # 13:F4 (p13030)540    family("e", "W3")541    chk("series duals S- == -K(p), S+ == -A(p) for p in [7, 13, 29, 101, 257, 1009]", dual_ok)542    chk("K(p) A(p) mod p = 3, 0, 21, 28, 132, 258 (never 1)", [v for _, v in prod_examples] == [3, 0, 21, 28, 132, 258])543    chk("asymmetry at p = 13: A(13) == 0 and K(13) == 10", alt_fact_mod(13) == 0 and kurepa_mod(13) == 10)544    # the p = 13 picture: the residue line of eps_n, n = 2..12545    # 13:F6 (p13032)546    family("e", "W4")547    d13 = d_mod(13, 13)548    line13 = [None if d13[n] == 0 else (math.factorial(n) % 13) * pow(d13[n], -1, 13) % 13 for n in range(2, 13)]549    log_e(f"  p=13 residue line [eps_n]_13, n=2..12: {['blind' if v is None else v for v in line13]}")550    chk("p=13: [eps_n]_13 for n = 2..12 is 2, 3, 7, 11, 1, 6, blind, 2, 8, 10, 9", line13 == [2, 3, 7, 11, 1, 6, None, 2, 8, 10, 9])551    chk("p=13: D_8 = 14833 = 13 * 1141 (the blind scale)", D[8] == 14833 == 13 * 1141)552    chk("p=13: the wall invariant 9 = -(K(13))^-1 = [6/5] = [11/7]", (-pow(kurepa_mod(13), -1, 13)) % 13 == 9 == 6 * pow(5, -1, 13) % 13 == 11 * pow(7, -1, 13) % 13)553554    # blind-scale statistics555    log_e()556    primes_stat = [p for p in sieve_e(5000) if 1000 < p < 5000]557    sizes = []558    for p in primes_stat:559        d = d_mod(p, p - 1)560        z = sum(1 for n in range(p) if (d[n] if n < len(d) else None) == 0)561        sizes.append(z)562    from collections import Counter563    cnt = Counter(sizes)564    mean = sum(sizes)/len(sizes)565    log_e(f"Blind-scale statistics over {len(primes_stat)} primes in (1000,5000):")566    log_e(f"  |Z0(p)| distribution: {dict(sorted(cnt.items()))},  mean = {mean:.3f}")567    log_e(f"  (n=1 is always blind since D_1=0; excess over 1 has mean {mean-1:.3f}, Poisson(1) predicts 1.000)")568    # 13:F5 (p13031)569    family("e", "S1")570    chk("501 primes in (1000, 5000)", len(primes_stat) == 501)571    chk("|Z0(p)| distribution 1:169, 2:209, 3:92, 4:22, 5:6, 6:2, 7:1", dict(cnt) == {1: 169, 2: 209, 3: 92, 4: 22, 5: 6, 6: 2, 7: 1})572    chk("mean |Z0(p)| = 1000/501 = 1.996", Fraction(sum(sizes), len(sizes)) == Fraction(1000, 501))573574    # =========================================================575    log_e(); log_e("="*72); log_e("[3] ANGULAR UNIT e_t = g^{i_t}: radian calibration over p < 10^6")576    log_e("="*72)577    LIM = 1000000578    primes = sieve_e(LIM)579    PI_LO, PI_HI = pi_bracket(80)                 # certified rational bracket of pi, width < 1e-100580    best = []581    count_001 = 0582    n_p1mod4 = 0583    near = []                                     # candidates for the best list, screened in floating point, decided exactly584    exact_hit = False                             # 13:H3: an integer label with 2 pi r = p - 1 would make pi rational; none below 10^6 (decided on the bracket)585    for p in primes:586        if p % 4 != 1: continue587        n_p1mod4 += 1588        lo_r, hi_r = Fraction(p - 1) / (2 * PI_HI), Fraction(p - 1) / (2 * PI_LO)589        if math.ceil(lo_r) <= math.floor(hi_r): exact_hit = True590        # sqrt(-1) mod p via a quadratic nonresidue591        a = 2592        while pow(a, (p-1)//2, p) != p - 1: a += 1593        r = pow(a, (p-1)//4, p)594        for rr in (r, p - r):595            # |theta - 1| < 1/100 with theta = 2 pi rr/(p-1): decided as the rational inequality596            # 99 (p-1) < 200 pi rr < 101 (p-1), pi replaced by its bracket (margin >> width)597            if 99*(p-1) < 200*PI_LO*rr and 200*PI_HI*rr < 101*(p-1): count_001 += 1598            delta_f = abs(2*math.pi*rr/(p-1) - 1.0)          # screening only599            if delta_f < 0.001: near.append((p, rr))600    # exact ordering of the candidates: delta = |2 pi rr/(p-1) - 1| with pi the bracket midpoint (the ordering is601    # decided at margins >> 1e-100)602    PI_MID = (PI_LO + PI_HI) / 2603    best = sorted((abs(2*PI_MID*rr/(p-1) - 1), p, rr) for p, rr in near)604    log_e(f"primes p==1 (mod 4) below 10^6: {n_p1mod4}")605    log_e(f"frames with |theta_t - 1 rad| < 0.01: {count_001}  (equidistribution predicts ~ {2*0.01/(2*math.pi)*n_p1mod4*2:.0f} over both orientations)")606    log_e("best radian-calibrated frames (p, i_t-label r, theta_t, |theta_t - 1|):")607    for delta, p, rr in best[:6]:608        log_e(f"  p={p:7d}  r={rr:7d}  theta={float(2*PI_MID*rr/(p-1)):.9f}  delta={float(delta):.3e}")609    # 13:H4 (p13069), 13:H3 (p13045)610    family("e", "A1")611    chk("39175 primes p == 1 (mod 4) below 10^6", n_p1mod4 == 39175)612    chk("exact calibration is impossible on every shell below 10^6: no integer label r with 2 pi r = p - 1, decided on the Machin bracket", not exact_hit)613    chk("exactly 260 shells with |theta - 1| < 0.01 (exact rational decision, pi as a Machin bracket)", count_001 == 260)614    chk("the best six (p, lambda(i)) in order: (778933,123966), (497773,79227), (962237,153133), (261017,41537), (345181,54945), (891997,141986)",615        [(p, rr) for _, p, rr in best[:6]] == [(778933, 123966), (497773, 79227), (962237, 153133), (261017, 41537), (345181, 54945), (891997, 141986)])616    chk("the best shell is calibrated below 4e-5 and above 3.9e-5", Fraction(39, 10**6) < best[0][0] < Fraction(4, 10**5))617618    # =========================================================619    log_e(); log_e("="*72); log_e("[4] NULL EXPERIMENT: pointwise g^{i_t} == k(e)-wall-residue?")620    log_e("="*72)621    def primitive_root(p):622        fac = []623        m = p - 1; d = 2624        while d*d <= m:625            if m % d == 0:626                fac.append(d)627                while m % d == 0: m //= d628            d += 1629        if m > 1: fac.append(m)630        for g in range(2, p):631            if all(pow(g, (p-1)//q, p) != 1 for q in fac):632                return g633        return None634635    solvable, tied, total = 0, 0, 0636    exp_free, var_free, exp_tied, var_tied = 0.0, 0.0, 0.0, 0.0637    corr_x, corr_y = [], []638    for p in [q for q in sieve_e(3000) if q % 4 == 1 and q > 5]:639        kp = kurepa_mod(p)640        if kp == 0: continue641        Ep = (-pow(kp, -1, p)) % p            # wall residue of e642        g0 = primitive_root(p)643        # index table644        ind = {}645        x = 1646        for m in range(p - 1):647            ind[x] = m648            x = x * g0 % p649        s = ind[Ep]650        a = 2651        while pow(a, (p-1)//2, p) != p - 1: a += 1652        r = pow(a, (p-1)//4, p)653        ok = (math.gcd(s, p-1) == math.gcd(r, p-1)) or (math.gcd(s, p-1) == math.gcd(p - r, p-1))654        # the two conventions (12 October 2026): the orientation of the quarter-turn free (either root of -1, the index655        # multiplier u unrestricted — the gcd criterion above), or tied to the generator as in Section 2, i = g^(-kappa), where656        # g = g0^u keeps the orientation for u = 1 (mod 4) and reverses it for u = 3 (mod 4); the reachable indices of each,657        # and the success probability of a uniform target index (the expectation under independence)658        i0 = pow(g0, (p-1) - (p-1)//4, p)659        reach_free, reach_tied = set(), set()660        for u in range(1, p-1):661            if math.gcd(u, p-1) != 1: continue662            reach_free.add(u * i0 % (p-1)); reach_free.add(u * (p - i0) % (p-1))663            reach_tied.add(u * i0 % (p-1) if u % 4 == 1 else u * (p - i0) % (p-1))664        assert ok == (s in reach_free)665        okt = s in reach_tied666        qf, qt = len(reach_free) / (p-1), len(reach_tied) / (p-1)667        exp_free += qf; var_free += qf * (1-qf); exp_tied += qt; var_tied += qt * (1-qt)668        solvable += ok; tied += okt; total += 1669        corr_x.append(min(r, p-r)/(p-1)); corr_y.append(s/(p-1))670    # Pearson correlation671    mx = sum(corr_x)/len(corr_x); my = sum(corr_y)/len(corr_y)672    num = sum((x-mx)*(y-my) for x, y in zip(corr_x, corr_y))673    den = math.sqrt(sum((x-mx)**2 for x in corr_x) * sum((y-my)**2 for y in corr_y))674    log_e(f"p==1(4), 5<p<3000: some frame g with g^(i_t) == -(!p)^(-1) exists for {solvable}/{total} primes ({100*solvable/total:.1f}%), "675          f"the orientation free: expected {exp_free:.1f} +- {math.sqrt(var_free):.1f} under independence (mean rate {100*exp_free/total:.1f}%)")676    log_e(f"with the orientation tied to the generator (i = g^-kappa): {tied}/{total} primes ({100*tied/total:.1f}%), "677          f"expected {exp_tied:.1f} +- {math.sqrt(var_tied):.1f} (mean rate {100*exp_tied/total:.1f}%)")678    log_e(f"Pearson correlation between angular address of i_t and of the wall residue: {num/den:+.4f}")679    log_e("=> solvability is gcd-coincidence noise, no invariant identification; the two objects are distinct projections.")680    # 13:H5 (p13047)681    family("e", "N1")682    chk("210 primes p == 1 (mod 4), 5 < p < 3000, with K(p) != 0", total == 210)683    chk("the index condition is solvable for exactly 94 of them (the orientation of the quarter-turn free)", solvable == 94)684    chk("expected under independence: 96.0 +- 6.6, the mean gcd-coincidence rate 45.7 %",685        round(exp_free, 1) == 96.0 and round(math.sqrt(var_free), 1) == 6.6 and round(100 * exp_free / total, 1) == 45.7)686    chk("with the orientation tied to the generator (i = g^-kappa, Section 2): solvable for exactly 54 of them", tied == 54)687    chk("expected under the same independence: 57.9 +- 6.2", round(exp_tied, 1) == 57.9 and round(math.sqrt(var_tied), 1) == 6.2)688    _DIAG[0] = f"Pearson correlation of the angular addresses {num/den:+.3f} (the paper's [approx] reading)"689690691    flush("e", order=["R1", "R2", "R3", "R4", "P1", "T1", "W1", "W2", "W3", "W4", "W5", "S1", "A1", "N1"], details={"N1": _DIAG[0]})692693# ------------------------------------------------------------------------------------------------------------694# block pi — validate_pi.py: FRC emergence of pi: verification suite.695# Sections:696#   [1] Radial chains: Wallis central-binomial pair (v_n, w_n) enclosure; Machin readout.697#   [2] Wall arithmetic: legibility window [1,2t], empty blind set, hard wall, residue -2.698#   [3] Second order: Morley congruence, W mod p^2 = Fermat-quotient formula, pi-Wieferich search.699#   [4] Lucas revivals and digitwise self-similarity.700#   [5] arcsin(1/2) chain: wall vanishing T(p) == 0 test and second-order pattern hunt.701#   [6] Angular exactness and Gauss-sum saturation.702#703#704# Package form (2026-09): the memorandum script as written, its PASS/FAIL lines and the paper's stated figures turned705# into registry predicates (families pi.R1-R4, pi.W1-W2, pi.S1-S3, pi.L1, pi.V1); the arcsin tail bound (pi.R4) and706# the universality of the first revival (pi.L1) are added. Every predicate is exact (integers, exact rationals; pi as707# the exact Machin bracket PI). The Fermat-quotient moments of [3], the pattern hunt of [5] and the Gauss sums of [6]708# are the paper's [approx] diagnostics: printed, deciding nothing.709OUT_pi = []710def log_pi(s=""):711    print(s); OUT_pi.append(str(s))712713def sieve_pi(N):714    c = bytearray(N + 1)715    for i in range(2, int(N**0.5) + 1):716        if not c[i]: c[i*i::i] = b'\x01' * len(c[i*i::i])717    return [i for i in range(3, N + 1) if not c[i]]718719# reference pi to ~420 digits via exact Machin720def pi_fraction(N=320):721    def arctan_inv(x, N):722        s = Fraction(0)723        for k in range(N + 1):724            s += Fraction((-1)**k, (2*k + 1) * x**(2*k + 1))725        return s726    return 16*arctan_inv(5, N) - 4*arctan_inv(239, N)727PI = pi_fraction(320)728729def Sg(N): return sum(Fraction((-1)**k, 2*k+1) for k in range(N+1))730731def block_pi():732    """Block pi — the Wallis, Machin and arcsin chains of π, the legibility window and the −2 terminus, Morley and the π-Wieferich search, Lucas revivals, the first-order arcsin vanishing (pi.R1 … pi.V1)."""733    # =========================================================734    log_pi("="*72); log_pi("[1] RADIAL CHAINS")735    log_pi("="*72)736    # Wallis pair737    def vw(n):738        C = comb(2*n, n)739        w = Fraction(16**n, n * C * C)740        v = Fraction(2 * 16**n, (2*n + 1) * C * C)741        return v, w742    encl = True; ident = True; mono = True743    prev_v, prev_w = vw(1)744    for n in range(1, 301):745        v, w = vw(n)746        if not (v < PI < w): encl = False747        if w - v != w / (2*n + 1): ident = False748        if n > 1 and not (v > prev_v and w < prev_w): mono = False749        if n > 1: prev_v, prev_w = v, w750    log_pi(f"(a) Wallis pair: v_n < pi < w_n for n=1..300: {'PASS' if encl else 'FAIL'}; "751        f"exact width identity w-v=w/(2n+1): {'PASS' if ident else 'FAIL'}; strict monotonicity: {'PASS' if mono else 'FAIL'}")752    v50, w50 = vw(50)753    log_pi(f"    example n=50: width = {float(w50-v50):.4e} ~ pi/101 = {math.pi/101:.4e}  (O(1/n) rate)")754    # 13:C5 (p13016), 13:E2 (p13023)755    family("pi", "R1")756    chk("v_n < pi < w_n for n = 1..300 (Machin bracket)", encl)757    chk("exact width identity w_n - v_n = w_n/(2n+1)", ident)758    chk("v_n strictly increasing, w_n strictly decreasing", mono)759    # p=13 showcase residue line (round-01: anchors the manuscript's displayed line)760    p13_line = [pow(16, n, 13) * pow(n * pow(math.comb(2*n, n), 2, 13) % 13, -1, 13) % 13 for n in range(1, 7)]761    log_pi(f"    p=13 residue line [w_n]_13, n=1..6: {p13_line} " + ("PASS" if p13_line == [4, 5, 10, 9, 6, 11] else "FAIL"))762    # 13:G8 (p13041)763    family("pi", "R2")764    chk("p=13: [w_n]_13 for n = 1..6 is 4, 5, 10, 9, 6, 11", p13_line == [4, 5, 10, 9, 6, 11])765766    # Machin chain: M_N with both arctans truncated at N767    def machin(N):768        def S(x, N): return sum(Fraction((-1)**k, (2*k+1)*x**(2*k+1)) for k in range(N+1))769        return 16*S(5, N) - 4*S(239, N)770    n_double = None771    for N in range(2, 25):772        if float(machin(N)) == math.pi:773            n_double = N; break774    log_pi(f"(b) Machin chain: minimal N with float(M_N) == IEEE-double(pi): N = {n_double}")775    err_bound_ok = all(abs(machin(N) - PI) < Fraction(17, (2*N+3)*5**(2*N+3)) for N in range(2, 40))776    log_pi(f"    error bound |M_N - pi| < 17/((2N+3) 5^(2N+3)) for N=2..39: {'PASS' if err_bound_ok else 'FAIL'}")777    N100 = None778    for N in range(60, 90):779        if abs(machin(N) - PI) < Fraction(1, 10**101): N100 = N; break780    log_pi(f"    100-decimal-digit determination at N = {N100}")781    # enclosure via alternating pairing782    def S_pair(x, N):  # (lower, upper) alternating bounds for arctan(1/x)783        s = Fraction(0); terms = [Fraction((-1)**k, (2*k+1)*x**(2*k+1)) for k in range(N+2)]784        partial = []785        run = Fraction(0)786        for t in terms:787            run += t; partial.append(run)788        lo = min(partial[N], partial[N+1]); hi = max(partial[N], partial[N+1])789        return lo, hi790    l5, u5 = S_pair(5, 12); l239, u239 = S_pair(239, 4)791    lo, hi = 16*l5 - 4*u239, 16*u5 - 4*l239792    log_pi(f"    fQ-native certified enclosure at (12,4) terms: pi in ({float(lo):.16f}, {float(hi):.16f}): {'PASS' if lo < PI < hi else 'FAIL'}")793    # 13:E4 (p13025), 13:E5 (p13026), 13:I2 (p13049)794    family("pi", "R3")795    chk("minimal N with float(M_N) == binary64(pi) is N = 10", n_double == 10)796    chk("binary64(pi) = 884279719003555/281474976710656 = float(M_10)", Fraction(math.pi) == Fraction(884279719003555, 281474976710656) and float(machin(10)) == math.pi)797    chk("|M_N - pi| < 17/((2N+3) 5^(2N+3)) for N = 2..39", err_bound_ok)798    chk("100 decimal digits fixed at N = 71 (|M_71 - pi| < 10^-101, first N with that property from 60)", N100 == 71)799    chk("certified enclosure at (12, 4) terms contains pi and has width below 2^-52", lo < PI < hi and hi - lo < Fraction(1, 2**52))800    chk("Gregory brackets 4 S_(2m+1) < pi < 4 S_(2m) of width 4/(4m+3), m <= 40", all(4*Sg(2*m+1) < PI < 4*Sg(2*m) and 4*Sg(2*m) - 4*Sg(2*m+1) == Fraction(4, 4*m+3) for m in range(1, 41)))801    # arcsin chain tail bound (Definition arcsin): |pi - s_n| < 4^-n/(2n+3) for 1 <= n <= 120802    # 13:E3 (p13024)803    family("pi", "R4")804    s_run = Fraction(0); tail_ok = True805    for n in range(0, 121):806        s_run += Fraction(comb(2*n, n), (2*n+1) * 16**n)807        if n >= 1 and not abs(PI - 3*s_run) < Fraction(1, 4**n * (2*n+3)): tail_ok = False808    chk("arcsin chain s_n = 3 sum C(2k,k)/((2k+1)16^k): |pi - s_n| < 4^-n/(2n+3) for n = 1..120", tail_ok)809810    # feasibility table811    log_pi("(c) precision-cost (terms n for 2^-k):")812    for k in (24, 53, 333):813        nL = 2**k            # Leibniz/Wallis ~ error 1/n814        nM = math.ceil((k*math.log(2)/math.log(5) - 3)/2)   # Machin ~ 5^{-2N}815        nA = math.ceil(k/2)  # arcsin(1/2) ~ 4^{-n}816        log_pi(f"    k={k:3d}: Wallis/Leibniz ~ {nL:.1e}   Machin ~ {nM}   arcsin(1/2) ~ {nA}   (e-chain reference: {[11,18,70][(24,53,333).index(k)]})")817818    # =========================================================819    log_pi(); log_pi("="*72); log_pi("[2] WALL ARITHMETIC (mod p)")820    log_pi("="*72)821    def Cmod_track(p, nmax):822        """return list over n=1..nmax of (vp, unit) for C(2n,n): p-adic valuation and unit part mod p"""823        res = []824        e, u = 0, 1  # C(2,1)=2 handled in loop start from n=0: C(0,0)=1825        # iterate C(2(n+1),n+1) = C(2n,n)*2*(2n+1)/(n+1)826        cur_e, cur_u = 0, 1827        for n in range(0, nmax):828            num = 2*(2*n + 1); den = n + 1829            while num % p == 0: num //= p; cur_e += 1830            while den % p == 0: den //= p; cur_e -= 1831            cur_u = cur_u * (num % p) % p832            cur_u = cur_u * pow(den % p, -1, p) % p833            res.append((cur_e, cur_u))  # this is C(2(n+1), n+1)834        return res835836    samples = [5, 13, 29, 101, 257, 1009, 10007]837    log_pi(f"{'p':>6} | window [1,2t] all legible | (2t,p) all blind | k(w_wall) | v blind at wall (den=p)")838    for p in samples:839        t2 = (p - 1)//2840        trk = Cmod_track(p, p - 1)841        legible = all(trk[n-1][0] == 0 for n in range(1, t2 + 1))842        blind   = all(trk[n-1][0] >= 1 for n in range(t2 + 1, p))843        e_, u_ = trk[t2 - 1]  # C(2*2t choose 2t) = C(p-1, (p-1)/2)844        wall = pow(16, t2, p) * pow(t2 % p, -1, p) * pow(u_*u_ % p, -1, p) % p if e_ == 0 else None845        log_pi(f"{p:6d} |          {'PASS' if legible else 'FAIL'}          |      {'PASS' if blind else 'FAIL'}      |   {(wall - p) if wall and wall > p//2 else wall}      |  2*2t+1 = {2*t2+1} = p: {'YES' if 2*t2+1==p else 'NO'}")846        # 13:G1 (p13034), 13:G2 (p13035)847        family("pi", "W1")848        chk(f"p={p}: p does not divide C(2n,n) for 1 <= n <= m", legible)849        chk(f"p={p}: p divides C(2n,n) for m < n < p", blind)850        chk(f"p={p}: [w_m] == -2 == m^-1 == 4 pi_A", wall == p - 2 == pow(t2, -1, p) == 4 * t2 % p)851        chk(f"p={p}: the lower chain is blind exactly at the wall, 2m+1 = p", 2*t2 + 1 == p)852853    # also p == 3 mod 4854    p3ok = True855    for p in [7, 23, 1031]:856        m = (p-1)//2857        C = 1858        for n in range(m): C = C * 2*(2*n+1) % p * pow(n+1, -1, p) % p859        wall = pow(16, m, p) * pow(m, -1, p) * pow(C*C % p, -1, p) % p860        if wall != p - 2: p3ok = False861    log_pi(f"universal wall residue -2 also for p == 3 (mod 4) [7, 23, 1031]: {'PASS' if p3ok else 'FAIL'}")862    family("pi", "W2")863    chk("[w_m] == -2 for p == 3 (mod 4), p in {7, 23, 1031}", p3ok)864865    # =========================================================866    log_pi(); log_pi("="*72); log_pi("[3] SECOND ORDER: Morley, Fermat quotient, pi-Wieferich")867    log_pi("="*72)868    morley_ok = True; formula_ok = True869    rows = []870    for p in [13, 29, 37, 41, 53, 101]:871        m = (p-1)//2872        C = comb(p-1, m)873        # Morley mod p^3874        if ((-1)**m * C - pow(4, p-1, p**3)) % p**3 != 0: morley_ok = False875        W = Fraction(16**m, m * C * C)876        k2 = W.numerator * pow(W.denominator, -1, p*p) % (p*p)877        q4 = (pow(4, p-1, p*p) - 1)//p % p878        pred = (-2 + 2*p*(q4 - 1)) % (p*p)879        if k2 != pred: formula_ok = False880        rows.append((p, (k2 - p*p) , q4))881    log_pi(f"Morley: (-1)^m C(p-1,m) == 4^(p-1) (mod p^3) for p in [13..101]: {'PASS' if morley_ok else 'FAIL'}")882    log_pi(f"Second-order wall formula W == -2 + 2p(q_p(4)-1) (mod p^2): {'PASS' if formula_ok else 'FAIL'}")883    log_pi(f"  (p, W mod p^2 as negative lift, q_p(4) mod p): {rows[:4]}")884    # 13:G4 (p13037)885    family("pi", "S1")886    chk("Morley (-1)^m C(p-1,m) == 4^(p-1) (mod p^3) for p in {13, 29, 37, 41, 53, 101}", morley_ok)887    chk("w_m == -2 + 2p(q_p(4) - 1) (mod p^2) for the same primes", formula_ok)888    chk("q_p(4) == 2 q_p(2) (mod p) for the same primes", all((pow(4, p-1, p*p) - 1)//p % p == 2 * ((pow(2, p-1, p*p) - 1)//p) % p for p in [13, 29, 37, 41, 53, 101]))889890    # pi-Wieferich search: 4^(p-1) == 1 + p (mod p^2)891    hits = []892    LIM = 1000000893    for p in sieve_pi(LIM):894        if pow(4, p-1, p*p) == (1 + p) % (p*p):895            hits.append(p)896    log_pi(f"pi-Wieferich primes (W == -2 mod p^2), p < 10^6: {hits if hits else 'none'}  "897        f"(heuristic expectation ~ sum 1/p ~ {sum(1/q for q in sieve_pi(LIM)):.2f} ... over the range)")898    family("pi", "S2")899    chk("pi-Wieferich primes below 10^6 are exactly {5, 45827}", hits == [5, 45827])900    # the recorded search of the standalone programs pi_wieferich.c and pi_wieferich128.c (results_wieferich.txt, 18 September901    # 2026): the recorded ranges tile [3, 2^35) and the recorded hits are 5 and 45827 alone902    rec = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "results_wieferich.txt"), encoding="utf-8").read()903    import re as _re904    spans = sorted((int(a), int(b)) for a, b in _re.findall(r"^range \[(\d+), (\d+)\):", rec, _re.M))905    rhits = [int(x) for x in _re.findall(r"^HIT p=(\d+)", rec, _re.M)]906    tiled = bool(spans) and spans[0][0] == 3 and spans[-1][1] == 2 ** 35 and all(spans[i][1] == spans[i + 1][0] for i in range(len(spans) - 1))907    counted = sum(int(n) for n in _re.findall(r"pi-Wieferich hits (\d+)", rec))908    log_pi(f"recorded search (results_wieferich.txt): ranges {spans}, hits {rhits}, hits counted {counted}")909    chk("the recorded search tiles [3, 2^35) with the hits 5 and 45827 alone (results_wieferich.txt: three runs, 2 hits counted)",910        tiled and rhits == [5, 45827] and counted == 2)911912    # the pi-Wieferich condition in three forms, Eisenstein's harmonic form and the OEIS A355959 form (13:G10)913    # 13:G10 (p13067)914    family("pi", "S3")915    S3 = [5, 13, 29, 37, 41, 53, 101, 45827]916    def fq(a, p): return (pow(a, p-1, p*p) - 1)//p             # the Fermat quotient q_p(a) reduced to [0, p)917    def fq_exact(a, p): return (a**(p-1) - 1)//p               # the Fermat quotient as the integer (a^(p-1) - 1)/p918    def H(m, p): return sum(pow(j, -1, p) for j in range(1, m+1)) % p919    chk("q_p(4) = 2 q_p(2) + p q_p(2)^2 exactly (as integers) for p in {5, 13, 29, 37, 41, 53, 101, 45827}",920        all(fq_exact(4, p) == 2*fq_exact(2, p) + p*fq_exact(2, p)**2 for p in S3))921    chk("the three forms agree on the same primes: 4^(p-1) == 1 + p (mod p^2) iff q_p(4) == 1 iff 2 q_p(2) == 1 (mod p)",922        all((pow(4, p-1, p*p) == (1 + p) % (p*p)) == (fq(4, p) % p == 1) == (2*fq(2, p) % p == 1) for p in S3))923    chk("Eisenstein: 2 q_p(2) == -H_((p-1)/2) (mod p), H_m = sum_(j<=m) 1/j, on the same primes",924        all(2*fq(2, p) % p == (-H((p-1)//2, p)) % p for p in S3))925    chk("OEIS A355959: (p+2)^(p-1) == 1 (mod p^2) iff 2 q_p(2) == 1 (mod p) on the same primes, holding at 5 and 45827 only",926        all((pow(p+2, p-1, p*p) == 1) == (2*fq(2, p) % p == 1) for p in S3)927        and [p for p in S3 if pow(p+2, p-1, p*p) == 1] == [5, 45827])928929    # distribution of q_p(4) mod p (uniformity check via normalized mean/var on p<20000)930    vals = []931    for p in sieve_pi(20000):932        q4 = (pow(4, p-1, p*p) - 1)//p % p933        vals.append(q4/p)934    mean = sum(vals)/len(vals); var = sum((x-mean)**2 for x in vals)/len(vals)935    log_pi(f"q_p(4)/p over p<20000: mean {mean:.4f} (uniform: 0.5), variance {var:.4f} (uniform: 0.0833)")936937    # =========================================================938    log_pi(); log_pi("="*72); log_pi("[4] LUCAS REVIVALS AND SELF-SIMILARITY")939    log_pi("="*72)940    p = 13941    Cv = comb(26, 13) % p942    vres = 2 * pow(16, 13, p) * pow((27 % p) * Cv * Cv % p, -1, p) % p943    log_pi(f"p=13: C(26,13) mod 13 = {Cv} = C(2,1)*C(0,0); k(v_13) = {vres} (predicted universal 8: {'PASS' if vres==8 else 'FAIL'})")944    C28 = comb(28, 14) % 13945    log_pi(f"p=13, n=14=(1,1)_13: C(28,14) mod 13 = {C28} = C(2,1)^2 = 4 digitwise (Lucas): {'PASS' if C28==4 else 'FAIL'}")946    # density of legible n in [p, p^2): digits all <= (p-1)/2 and n not = 0 mod p etc.947    p = 13; cnt = 0948    for n in range(p, p*p):949        a, b = divmod(n, p)950        if a <= 6 and b <= 6 and n % p != 0 and comb(2*n, n) % p != 0: cnt += 1951    log_pi(f"p=13: legible w-scales in [13,169): {cnt}/156 = {cnt/156:.3f} (prediction ((2t+1)/p)^2-ish = {(7/13)**2:.3f}, minus n==0 mod p column)")952    # 13:G5 (p13038)953    family("pi", "L1")954    chk("p=13: C(26,13) == 2 == C(2,1) C(0,0) (mod 13) and [v_13] == 8", Cv == 2 and vres == 8)955    chk("p=13: C(28,14) == 4 == C(2,1)^2 (mod 13), Lucas digitwise", C28 == 4)956    chk("p=13: two-digit revival scales in [13, 169) number exactly (2 kappa)^2 = 36", cnt == 36)957    chk("[v_p] == 8 universally, p in {5, 7, 11, 13, 29, 37}", all(2 * pow(16, q, q) * pow((2*q+1) % q * pow(comb(2*q, q), 2, q) % q, -1, q) % q == 8 % q for q in [5, 7, 11, 13, 29, 37]))958959    # =========================================================960    log_pi(); log_pi("="*72); log_pi("[5] ARCSIN(1/2) CHAIN WALL: T(p) and second order")961    log_pi("="*72)962    def T_mod(p):963        m = (p-1)//2964        s = 0; C = 1  # C(0,0)965        inv16 = pow(16, -1, p)966        pw = 1967        for k in range(m):968            if k > 0:969                C = C * 2*(2*k-1) % p * pow(k, -1, p) % p970                pw = pw * inv16 % p971            s = (s + C * pw % p * pow(2*k+1, -1, p)) % p972        return s973    zero_all = True; tested = 0974    for p in sieve_pi(3000):975        if p == 3:976            if T_mod(p) == 0: zero_all = False; log_pi("  T(3) == 0 (expected nonzero: sole exception)!")977            continue978        if T_mod(p) != 0: zero_all = False; log_pi(f"  T({p}) != 0 !")979        tested += 1980    log_pi(f"T(p) := sum_(k=0)^(m-1) C(2k,k)/((2k+1)16^k) == 0 (mod p) for all odd primes 5<=p<3000 (p=3 is the sole exception, T(3)={T_mod(3)}): {'PASS' if zero_all else 'FAIL'} ({tested} primes)")981    # 13:G6 (p13039)982    family("pi", "V1")983    chk("sigma_p == 0 (mod p) for all 428 odd primes 5 <= p < 3000", zero_all and tested == 428)984    chk("sigma_3 = 1, the sole exception", T_mod(3) == 1)985986    # second order: sigma/p mod p, pattern hunt987    def euler_numbers_mod(p, upto):988        E = [0]*(upto+1); E[0] = 1989        for n in range(1, upto//2 + 1):990            s = 0991            for j in range(n):992                s += comb(2*n, 2*j) * E[2*j]993            E[2*n] = (-s) % p994        return E995    def two_squares(p):996        # p = a^2 + b^2, p == 1 mod 4; a odd normalized a == 1 mod 4? return (a,b) a odd997        for a in range(1, int(p**0.5)+1):998            b2 = p - a*a999            b = int(b2**0.5)1000            if b*b == b2: 1001                if a % 2 == 1: return a, b1002                return b, a1003        return None1004    log_pi("second-order invariant T2(p) := (sigma/p) mod p, sigma the exact rational partial sum:")1005    log_pi(f"{'p':>4} | T2 | E_(p-3) | q_p(2) | a (p=a^2+b^2) | T2/E | T2/q2")1006    pat = []1007    for p in [q for q in sieve_pi(200) if q > 3]:1008        m = (p-1)//21009        sigma = sum(Fraction(comb(2*k,k), (2*k+1)*16**k) for k in range(m))1010        num, den = sigma.numerator, sigma.denominator1011        chk(f"p={p}: the rational sigma_p has numerator divisible by p and denominator prime to p", num % p == 0 and den % p != 0)1012        T2 = (num//p) * pow(den, -1, p) % p1013        E = euler_numbers_mod(p, p-3)[p-3]1014        q2 = (pow(2, p-1, p*p) - 1)//p % p1015        ts = two_squares(p) if p % 4 == 1 else None1016        rE = T2 * pow(E, -1, p) % p if E else None1017        rq = T2 * pow(q2, -1, p) % p if q2 else None1018        pat.append((p, T2, E, q2, ts, rE, rq))1019        log_pi(f"{p:4d} | {T2:3d} | {E:5d} | {q2:4d} | {str(ts):>10} | {str(rE):>4} | {str(rq):>4}")10201021    # also Leibniz wall U(p) and ratio to E_(p-3)1022    log_pi()1023    log_pi("Leibniz wall U(p) := sum_(k=0)^(m-1) (-1)^k/(2k+1) mod p, ratio to E_(p-3):")1024    for p in [q for q in sieve_pi(120) if q > 3]:1025        m = (p-1)//21026        U = sum((-1)**k * pow(2*k+1, -1, p) for k in range(m)) % p1027        E = euler_numbers_mod(p, p-3)[p-3]1028        r = U * pow(E, -1, p) % p if E else None1029        log_pi(f"  p={p:3d}: U={U:3d}  E_(p-3)={E:4d}  U/E={r}")10301031    # =========================================================1032    log_pi(); log_pi("="*72); log_pi("[6] ANGULAR EXACTNESS AND GAUSS SUM")1033    log_pi("="*72)1034    for p in [13, 101, 1009, 10007]:1035        re = sum(math.cos(2*math.pi*(x*x % p)/p) for x in range(p))1036        im = sum(math.sin(2*math.pi*(x*x % p)/p) for x in range(p))1037        log_pi(f"  p={p:6d}: Gauss sum = {re:.6f} + {im:.2e} i;  sqrt(p) = {math.sqrt(p):.6f}  (class: +sqrt(p) if p==1 mod 4, +i*sqrt(p) if p==3 mod 4)")1038    log_pi("angle of -1 on the multiplicative circle: 2pi * (2t)/(4t) = pi exactly, every frame, both chiralities.")103910401041    flush("pi", order=["R1", "R2", "R3", "R4", "W1", "W2", "S1", "S2", "S3", "L1", "V1"])10421043# ------------------------------------------------------------------------------------------------------------1044# block pi2 — validate_pi2.py: Follow-up: quarter-wall two-squares invariant, arcsin supercongruence, proof sanity checks.1045#1046# Package form (2026-09): the script as written, its PASS/FAIL lines turned into registry predicates1047# (families pi2.A1-A2, pi2.B1-B2, pi2.C1-C4, pi2.D1-D2). Exact arithmetic only: integers and fractions.Fraction.1048OUT_pi2 = []1049def log_pi2(s=""):1050    print(s); OUT_pi2.append(str(s))10511052def sieve_pi2(N):1053    c = bytearray(N + 1)1054    for i in range(2, int(N**0.5) + 1):1055        if not c[i]: c[i*i::i] = b'\x01' * len(c[i*i::i])1056    return [i for i in range(3, N + 1) if not c[i]]10571058def two_squares(p):1059    for a in range(1, int(p**0.5) + 1):1060        b = int((p - a*a)**0.5)1061        if a*a + b*b == p:1062            return (a, b) if a % 2 == 1 else (b, a)1063    return None10641065def block_pi2():1066    """Block pi2 — the two-squares quarter invariant, Sun's supercongruence and the Bernoulli law, the proof ingredients, the blind-range Euler congruence, the revival to p² (pi2.A1 … pi2.D2)."""1067    # ---- [A] Quarter-wall: Gauss congruence C(2t,t) == 2a (mod p), a odd, a == 1 (mod 4);1068    #          and k(w_t) == -(a^2)^{-1} == (b^2)^{-1} (mod p)1069    log_pi2("[A] QUARTER-WALL TWO-SQUARES INVARIANT, all p == 1 (mod 4), p < 3000")1070    ok_g, ok_w, cnt = True, True, 01071    for p in [q for q in sieve_pi2(3000) if q % 4 == 1]:1072        t = (p - 1)//41073        a, b = two_squares(p)1074        astar = a if a % 4 == 1 else -a           # Gauss normalization1075        C = 11076        for n in range(t): C = C * 2*(2*n+1) % p * pow(n+1, -1, p) % p   # C(2t,t) mod p1077        if C != (2*astar) % p: ok_g = False; log_pi2(f"  Gauss fails at p={p}")1078        # 13:G3 (p13036)1079        chk(f"p={p}: C(2t,t) == 2a* (mod p)", C == (2*astar) % p, "pi2", "A1")1080        kw = pow(16, t, p) * pow(t % p, -1, p) % p * pow(C*C % p, -1, p) % p1081        if kw != (-pow(a*a % p, -1, p)) % p or kw != pow(b*b % p, -1, p) % p:1082            ok_w = False; log_pi2(f"  wall formula fails at p={p}")1083        chk(f"p={p}: [w_t] == -(a^2)^-1 == (b^2)^-1 (mod p)", kw == (-pow(a*a % p, -1, p)) % p == pow(b*b % p, -1, p) % p, "pi2", "A2")1084        cnt += 11085    log_pi2(f"  Gauss congruence C(2t,t) == 2a* (mod p): {'PASS' if ok_g else 'FAIL'} ({cnt} primes)")1086    log_pi2(f"  k(w_t) == -(a^2)^(-1) == (b^2)^(-1) (mod p): {'PASS' if ok_w else 'FAIL'}")1087    p = 13; a, b = two_squares(13)1088    log_pi2(f"  example p=13 = {a}^2+{b}^2: k(w_3) = {pow(16,3,13)*pow(3,-1,13)*pow(pow(comb(6,3),2,13),-1,13)%13} = -(9)^(-1) = 10")1089    chk("211 primes p == 1 (mod 4) below 3000", cnt == 211, "pi2", "A1")1090    chk("p=13 = 3^2 + 2^2: [w_3] = 10 = -(9)^-1", (a, b) == (3, 2) and pow(16,3,13)*pow(3,-1,13)*pow(pow(comb(6,3),2,13),-1,13) % 13 == 10 == (-pow(9, -1, 13)) % 13, "pi2", "A2")10911092    # ---- [B] arcsin supercongruence: sigma == 0 (mod p^2) and third-order invariant1093    log_pi2(); log_pi2("[B] ARCSIN CHAIN sigma_p := sum_(k=0)^((p-3)/2) C(2k,k)/((2k+1)16^k)")1094    ok2 = True1095    rows = []1096    def euler_mod(p, upto):1097        E = [0]*(upto+1); E[0] = 1 % p1098        for n in range(1, upto//2 + 1):1099            s = sum(comb(2*n, 2*j) * E[2*j] for j in range(n))1100            E[2*n] = (-s) % p1101        return E1102    def bernoulli_upto(N):1103        """Exact Bernoulli numbers B_0..B_N via sum C(n+1,k) B_k = 0."""1104        B = [Fraction(0)] * (N + 1); B[0] = Fraction(1)1105        for n in range(1, N + 1):1106            B[n] = -sum(Fraction(comb(n+1, k)) * B[k] for k in range(n)) / (n + 1)1107        return B1108    BERN = bernoulli_upto(297)1109    for p in [q for q in sieve_pi2(300) if q >= 5]:1110        m = (p-1)//21111        sig = sum(Fraction(comb(2*k, k), (2*k+1)*16**k) for k in range(m))1112        num, den = sig.numerator, sig.denominator1113        if num % (p*p) != 0: ok2 = False; log_pi2(f"  p^2 fails at p={p}")1114        chk(f"p={p}: sigma_p == 0 (mod p^2)", num % (p*p) == 0 and den % p != 0, "pi2", "B1")1115        T3 = (num//(p*p)) % p * pow(den % p, -1, p) % p1116        Bpm3 = BERN[p-3]1117        bmod = Bpm3.numerator % p * pow(Bpm3.denominator % p, -1, p) % p1118        law = (-1)**((p+1)//2) * bmod % p * pow(36, -1, p) % p1119        rows.append((p, T3, law))1120        # 13:G7 (p13040)1121        chk(f"p={p}: sigma_p/p^2 == (-1)^((p+1)/2) B_(p-3)/36 (mod p)", T3 == law, "pi2", "B2")1122    log_pi2(f"  sigma == 0 (mod p^2) for all 5 <= p < 300: {'PASS' if ok2 else 'FAIL'}")1123    okB3 = all(T3 == law for _, T3, law in rows)1124    log_pi2(f"  third-order Bernoulli law sigma/p^2 == (-1)^((p+1)/2) B_(p-3)/36 (mod p), 5<=p<300: "1125        f"{'PASS' if okB3 else 'FAIL'} ({len(rows)} primes)  [mod-p^3 reading of Sun Conj. 5.1 + Wolstenholme refinement]")1126    log_pi2("  p, sigma/p^2 mod p, law: " + str([(p, t, l) for p, t, l in rows[:8]]))1127    chk("sixty primes 5 <= p < 300", len(rows) == 60, "pi2", "B2")11281129    # ---- [B3] the third-order law to 10^4, modular throughout (the method of reports/daily-push/t8-push-20260917)1130    def sigma_mod(p, M):1131        """sigma_p mod M (M a power of p): C(2k,k) by the recursion C(2k,k) = C(2k-2,k-1)*2(2k-1)/k, every divisor a p-unit for 2k+1 <= p-2."""1132        c = 1; s = 0; sixteen = 11133        for k in range((p - 1) // 2):1134            if k:1135                c = c * 2 * (2*k - 1) * pow(k, -1, M) % M1136                sixteen = sixteen * 16 % M1137            s = (s + c * pow((2*k + 1) * sixteen, -1, M)) % M1138        return s1139    def bern_mod_p(p, n):1140        """B_n mod p for even n <= p-3: sum_{j=1}^{p-1} j^n == p B_n (mod p^2) (Faulhaber, every denominator a p-unit)."""1141        M = p * p1142        S = sum(pow(j, n, M) for j in range(1, p)) % M1143        return (S // p) % p if S % p == 0 else None1144    okL = True; nL = 0; first = []; agree = True1145    for p in [q for q in sieve_pi2(10**4) if q >= 5]:1146        s3 = sigma_mod(p, p**3); b3 = bern_mod_p(p, p - 3)1147        lhs = (s3 // (p*p)) % p if s3 % (p*p) == 0 else None1148        rhs = None if b3 is None else (-1)**((p+1)//2) * b3 % p * pow(36, -1, p) % p1149        if lhs is None or rhs is None or lhs != rhs: okL = False1150        if p < 300:1151            Bpm3 = BERN[p-3]1152            if b3 != Bpm3.numerator % p * pow(Bpm3.denominator % p, -1, p) % p: agree = False1153        if len(first) < 10: first.append((p, lhs))1154        nL += 11155    log_pi2(f"  third-order law to 10^4, sigma_p mod p^3 and B_(p-3) mod p modular: {'PASS' if okL else 'FAIL'} ({nL} primes)")1156    chk("1227 primes 5 <= p < 10^4", nL == 1227, "pi2", "B3")1157    chk("sigma_p == 0 (mod p^2) and sigma_p/p^2 == (-1)^((p+1)/2) B_(p-3)/36 (mod p) for all 1227 primes 5 <= p < 10^4", okL, "pi2", "B3")1158    chk("B_(p-3) mod p by the power sum agrees with the exact Bernoulli recurrence for the sixty primes 5 <= p < 300", agree, "pi2", "B3")1159    chk("first rows (p, sigma_p/p^2 mod p) = (5,4) (7,3) (11,5) (13,6) (17,15) (19,2) (23,10) (29,21) (31,9) (37,2)",1160        first == [(5,4),(7,3),(11,5),(13,6),(17,15),(19,2),(23,10),(29,21),(31,9),(37,2)], "pi2", "B3")11611162    # ---- [C] proof sanity checks1163    log_pi2(); log_pi2("[C] PROOF INGREDIENTS")1164    okc = True1165    for p in [11, 13, 29, 37]:1166        m = (p-1)//21167        for k in range(m):1168            if comb(2*k, k) % p != pow(-4, k, p) * comb(m, k) % p: okc = False1169    log_pi2(f"  C(2k,k) == (-4)^k binom(m,k) (mod p), k < m: {'PASS' if okc else 'FAIL'}")1170    chk("C(2k,k) == (-4)^k C(m,k) (mod p) for k < m, p in {11, 13, 29, 37}", okc, "pi2", "C1")1171    okl = True1172    for p in [q for q in sieve_pi2(500) if q >= 5]:1173        m = (p-1)//21174        H = sum(pow(j, -1, p) for j in range(1, m+1)) % p1175        q2 = (pow(2, p-1, p*p) - 1)//p % p1176        if H != (-2*q2) % p: okl = False1177    log_pi2(f"  Lerch: H_((p-1)/2) == -2 q_p(2) (mod p): {'PASS' if okl else 'FAIL'}")1178    chk("Lerch H_((p-1)/2) == -2 q_p(2) (mod p) for 5 <= p < 500", okl, "pi2", "C2")1179    okI = True1180    for m in range(1, 61):1181        L = sum(Fraction(comb(m, j)*(-1)**j, 2*j+1) for j in range(m+1))1182        if L != Fraction(4**m * math.factorial(m)**2, math.factorial(2*m+1)): okI = False1183    log_pi2(f"  L_m = sum binom(m,j)(-1)^j/(2j+1) = 4^m (m!)^2/(2m+1)!, m<=60: {'PASS' if okI else 'FAIL'}")1184    chk("L_m = 4^m (m!)^2/(2m+1)! over Q for m <= 60", okI, "pi2", "C3")1185    okB = True1186    for m in range(1, 61):1187        y = Fraction(-1, 4)1188        B = sum(comb(m, k)*y**k/(k+m+1) for k in range(m+1))1189        A = sum(comb(m, k)*y**k/(2*k+1) for k in range(m+1))1190        L = sum(Fraction(comb(m, j)*(-1)**j, 2*j+1) for j in range(m+1))1191        # finite proof route (round-01): F = formal antiderivative of expanded (1-s^2)^m;1192        # A = 2(F(1/2)-F(0)), B = 2(F(1)-F(1/2)), 2L = 2(F(1)-F(0)) -- coefficient bookkeeping only1193        def F(s):1194            return sum(Fraction(comb(m, j)*(-1)**j) * s**(2*j+1) / (2*j+1) for j in range(m+1))1195        if A + B != 2*L: okB = False1196        if A != 2*(F(Fraction(1, 2)) - F(Fraction(0))): okB = False1197        if B != 2*(F(Fraction(1)) - F(Fraction(1, 2))): okB = False1198        if L != F(Fraction(1)) - F(Fraction(0)): okB = False1199    log_pi2(f"  key identity A_m + B_m = 2 L_m over Q, via formal-antiderivative bookkeeping, m<=60: {'PASS' if okB else 'FAIL'}")1200    chk("A_m + B_m = 2 L_m with A = 2(F(1/2)-F(0)), B = 2(F(1)-F(1/2)), L = F(1)-F(0) for m <= 60", okB, "pi2", "C4")12011202    # ---- [D] blind-range Euler congruence and refined revival residue1203    log_pi2(); log_pi2("[D] BLIND-RANGE AND REVIVAL REFINEMENTS")1204    okU = True1205    for p in [q for q in sieve_pi2(80) if q >= 5]:1206        m = (p-1)//21207        Us = sum(Fraction(comb(2*k, k), (2*k+1)*16**k) for k in range(m+1, p))1208        E = euler_mod(p, p-3)[p-3]1209        lhs = Us.numerator * pow(Us.denominator, -1, p*p) % (p*p)1210        rhs = p*E % (p*p) * pow(3, -1, p*p) % (p*p)1211        if lhs != rhs: okU = False1212    log_pi2(f"  blind-range sum == p*E_(p-3)/3 (mod p^2), 5<=p<80: {'PASS' if okU else 'FAIL'}  [known: van Hamme-Sun family]")1213    chk("blind-range sum == p E_(p-3)/3 (mod p^2) for 5 <= p < 80", okU, "pi2", "D1")1214    okR = True1215    for p in [5, 7, 11, 13, 29, 37]:1216        v = Fraction(2*16**p, (2*p+1)*comb(2*p, p)**2)1217        lhs = v.numerator * pow(v.denominator, -1, p*p) % (p*p)1218        q2 = (pow(2, p-1, p*p) - 1)//p % p1219        if lhs != (8 + 16*p*(2*q2 - 1)) % (p*p): okR = False1220    log_pi2(f"  first revival v_p == 8 + 16p(2q_p(2)-1) (mod p^2): {'PASS' if okR else 'FAIL'}  [via Wolstenholme]")1221    chk("v_p == 8 + 16p(2 q_p(2) - 1) (mod p^2) for p in {5, 7, 11, 13, 29, 37}", okR, "pi2", "D2")122212231224    flush("pi2", order=["A1", "A2", "B1", "B2", "B3", "C1", "C2", "C3", "C4", "D1", "D2"])12251226# ------------------------------------------------------------------------------------------------------------1227# block tow — validate_towers.py: validate_towers.py -- 13-epi fixed-shell towers, Cayley identities, orientation rule.1228# Exact arithmetic only: integers and fractions.Fraction. No floats anywhere.1229# External comparisons use the paper's own certified rational brackets.1230# Package form (2026-09): the paper's script as written, its micro-checks reported to the registry under the1231# families tow.E (e-tower), tow.P (π-tower), tow.C (Cayley), tow.O (orientation), tow.H (height run), tow.W (window1232# and pins); a failing micro-check_tow prints and fails its family without stopping the run.1233def check_tow(name, cond):1234    chk(name, cond)12351236def subfact(n):1237    d = [1, 0]1238    for k in range(2, n+1):1239        d.append(k*d[-1] + (-1)**k)1240    return d[n]12411242def C(n, k):1243    return math.comb(n, k)12441245def e_brackets(n):1246    """Certified bracket for e via the subfactorial chain: eps_{2m} < e < eps_{2m+1}."""1247    lo = F(math.factorial(2*n), subfact(2*n))1248    hi = F(math.factorial(2*n+1), subfact(2*n+1))1249    return lo, hi12501251def pi_brackets(m):1252    """Certified Machin brackets: L_m < pi < U_m (alternating pairing)."""1253    def AN(x, N):1254        return sum(F((-1)**k, (2*k+1)*x**(2*k+1)) for k in range(N+1))1255    L = 16*AN(5, 2*m+1) - 4*AN(239, 2*m)1256    U = 16*AN(5, 2*m)   - 4*AN(239, 2*m+1)1257    return L, U12581259def centered(x, p):1260    x %= p1261    return x - p if x > (p-1)//2 else x12621263def primitive_root(p):1264    fac = set()1265    n = p-11266    d = 21267    while d*d <= n:1268        while n % d == 0:1269            fac.add(d); n //= d1270        d += 11271    if n > 1: fac.add(n)1272    for g in range(2, p):1273        if all(pow(g, (p-1)//q, p) != 1 for q in fac):1274            return g1275    raise ValueError12761277def block_tow():1278  """Block tow — the fixed-shell towers, the Cayley map, orientation transport, the height run, the wrap-free window and the pins (tow.E … tow.W)."""1279  for p in (13, 29):1280    kap = (p-1)//41281    g = primitive_root(p)1282    i = pow(g, p-1-kap, p)              # i = g^{-kap}1283    # 13:H2 (p13044), 13:I4 (p13051), 13:J1 (p13052)1284    family("tow", "W")1285    check_tow(f"p={p}: i^2 = -1", (i*i) % p == p-1)1286    lam = i                              # label lift: canonical representative read as exponent1287    eP = pow(g, lam, p)1288    piA = (2*kap) % p12891290    # --- e-tower: q_m = ((mp)! + delta_m)/!(mp), shell reading = eP, external -> e ---1291    # 13:B3 (p13010), 13:F7 (p13033)1292    family("tow", "E")1293    for m in (1, 2, 3):1294        delta = centered(((-1)**m * eP) % p, p)1295        num = math.factorial(m*p) + delta1296        den = subfact(m*p)1297        check_tow(f"p={p} e-tower grade {m}: !(mp) == (-1)^m (mod p)", den % p == ((-1)**m) % p)1298        check_tow(f"p={p} e-tower grade {m}: shell reading = e_p", (num * pow(den % p, p-2, p)) % p == eP)1299        lo, hi = e_brackets(max(m*p//2, 12))1300        q = F(num, den)1301        # the paper's certificate: |q - e| < (mp)!/( !(mp) !(mp+1) ) + pi_F/!(mp)1302        width = F(math.factorial(m*p), den*subfact(m*p+1)) + F((p-1)//2, den)1303        check_tow(f"p={p} e-tower grade {m}: |q - e| inside certified window",1304              q > lo - width and q < hi + width)13051306    # --- pi-tower: (2*16^n + delta)/((2n+1) C(2n,n)^2), n = p^r ---1307    # 13:G9 (p13042)1308    family("tow", "P")1309    for r in (1, 2):1310        n = p**r1311        Cb = C(2*n, n)1312        Bn = (2*n+1) * Cb * Cb1313        check_tow(f"p={p} pi-tower r={r}: C(2n,n) == 2 (mod p)", Cb % p == 2)1314        check_tow(f"p={p} pi-tower r={r}: 16^n == 16 (mod p)", pow(16, n, p) == 16 % p)1315        delta = centered((piA * Bn - 2*pow(16, n)) % p, p)1316        check_tow(f"p={p} pi-tower r={r}: delta = centered(-34)", delta == centered(-34 % p, p))1317        num = 2*pow(16, n) + delta1318        check_tow(f"p={p} pi-tower r={r}: shell reading = pi_A",1319              (num * pow(Bn % p, p-2, p)) % p == piA)1320        if r == 1:1321            L, U = pi_brackets(8)1322            v = F(num, Bn)1323            check_tow(f"p={p} pi-tower r=1: member below the certified upper bracket", v < U)1324            check_tow(f"p={p} pi-tower r=1: member above v_n floor", v > F(2*4**n, (2*n+1)*Cb*Cb) - F(1))13251326    # --- Cayley identities, exhaustive where denominators are units ---1327    # 13:C3 (p13014)1328    family("tow", "C")1329    cnt = 01330    for x in range(p):1331        for y in range(p):1332            dx, dy = (1 - i*x) % p, (1 - i*y) % p1333            dxy = (1 - x*y) % p1334            if dx and dy and dxy:1335                num_c = ((x + y) * pow(dxy, p-2, p)) % p1336                dc = (1 - i*num_c) % p1337                if dc:1338                    Cx = ((1 + i*x) * pow(dx, p-2, p)) % p1339                    Cy = ((1 + i*y) * pow(dy, p-2, p)) % p1340                    Cc = ((1 + i*num_c) * pow(dc, p-2, p)) % p1341                    if (Cx*Cy) % p != Cc:1342                        check_tow(f"p={p} Cayley composition at ({x},{y})", False)1343                    cnt += 11344    check_tow(f"p={p} Cayley composition exhaustive ({cnt} pairs)", True)1345    check_tow(f"p={p} Cayley C(0)=1", ((1) * pow(1, p-2, p)) % p == 1)1346    check_tow(f"p={p} Cayley C(1)=i", ((1+i) * pow((1-i) % p, p-2, p)) % p == i)13471348    # --- orientation rule: u = 1 mod 4 preserves i, u = 3 mod 4 conjugates ---1349    # 13:B2 (p13009)1350    family("tow", "O")1351    for u in range(1, p-1):1352        if math.gcd(u, p-1) != 1:1353            continue1354        gp = pow(g, u, p)1355        ip = pow(gp, p-1-kap, p)1356        if u % 4 == 1:1357            if ip != i: check_tow(f"p={p} orientation u={u} (u=1 mod 4) preserves i", False)1358        else:1359            if ip != (p - i) % p: check_tow(f"p={p} orientation u={u} (u=3 mod 4) conjugates i", False)1360    check_tow(f"p={p} orientation-transport rule exhaustive over units", True)136113621363  _run_tail()13641365# --- two-level horizon law: minimal framed-rational heights ---1366def minheight(x, p):1367    best = p1368    for b in range(1, best+1):1369        a = (x*b) % p1370        a = a - p if a > p//2 else a1371        h = max(abs(a), b)1372        if h < best: best = h1373        if b > best: break1374    return best13751376# Full-population height run (round-02): all shells p == 1 (mod 4), p <= 8009;1377# smallest-primitive-root frame, H minimized over the two chirality units.1378def _sieve(N):1379    c = bytearray(N+1)1380    for i in range(2, _m.isqrt(N)+1):1381        if not c[i]: c[i*i::i] = b'\x01'*len(c[i*i::i])1382    return [i for i in range(3, N+1) if not c[i]]1383def _run_tail():1384  # 13:J2 (p13053), 13:J3 (p13054), 13:J5 (p13056)1385  family("tow", "H")1386  hs = []1387  pin_ok = True; band_ok = True1388  for p in [q for q in _sieve(8009) if q % 4 == 1]:1389      kap = (p-1)//41390      g = primitive_root(p)                      # iterates from 2: the smallest primitive root1391      e_neg = pow(g, pow(g, p-1-kap, p), p)      # exponential unit, oriented chirality  i = g^{-kap}1392      e_pos = pow(g, pow(g, kap, p), p)          # exponential unit, conjugate chirality i' = g^{+kap}1393      piA = (2*kap) % p1394      if not ((piA * (-2)) % p == 1 and minheight(piA, p) == 2): pin_ok = False1395      h = min(minheight(e_neg, p), minheight(e_pos, p))1396      if h > 2*_m.isqrt(p): band_ok = False1397      hs.append((p, h, _m.isqrt(p)))1398  check_tow("pi_A = [-1/2] pinned at height 2 on all 500 shells p <= 8009", pin_ok and len(hs) == 500)1399  check_tow("H(e_p) within the horizon band 2 sqrt(p) on every shell", band_ok)1400  small = sorted(p for p, h, s in hs if h <= 10)1401  exc2 = sorted(p for p, h, s in hs if h == 2)1402  ratios = sorted(h/s for p, h, s in hs)1403  med = (ratios[249] + ratios[250]) / 21404  _MED[0] = med1405  check_tow("small-height shells (H <= 10): exactly 68", len(small) == 68)1406  check_tow("height-2 shells exactly {13, 1933, 4177, 5857}", exc2 == [13, 1933, 4177, 5857])1407  # growth reading (exact counts; the generic-growth interpretation is DFI-equidistribution [import])1408  print(f"height run: 500 shells, median H/sqrt(p) = {med:.3f}, H<=10: {len(small)}, height-2: {exc2}")140914101411  # --- accessible window: wrap-free below sqrt(p); pinning relations; the half-turn tautology ---1412  family("tow", "W")1413  for p in (13, 29, 101, 997, 8009):1414      s = _m.isqrt(p)1415      check_tow(f"p={p}: wrap-free window (isqrt(p)^2 < p, sums below p)", s*s < p and 2*s < p)1416      kap = (p-1)//41417      piA = (2*kap) % p1418      check_tow(f"p={p}: pinning relation 2*pi_A + 1 == 0 (the shell calibration)", (2*piA + 1) % p == 0)1419      g = primitive_root(p)1420      check_tow(f"p={p}: half-turn tautology g^(2 kap) == -1 (the exact angular carrier of pi)", pow(g, 2*kap, p) == p - 1)1421      i = pow(g, p-1-kap, p)1422      check_tow(f"p={p}: quarter-turn pin i^2 + 1 == 0", (i*i + 1) % p == 0)1423  # --- the duality on the shells p = 13 and 29 (13:I3): pi angular-exact and wall-trivial, e angular-approximate and wall-wild ---1424  PI_LO2, PI_HI2 = pi_bracket(40)1425  # 13:I3 (p13050)1426  family("tow", "D")1427  for p in (13, 29):1428      kap = (p-1)//4; m = (p-1)//2; g = primitive_root(p)1429      cbm = [_m.comb(2*n, n) for n in range(p)]1430      wall = pow(16, m, p) * pow(m * pow(cbm[m], 2, p), -1, p) % p                    # [w_m] = 16^m/(m C_m^2)1431      a = next(x for x in range(1, p) if (p - x*x) >= 0 and _m.isqrt(p - x*x)**2 == p - x*x and x % 2 == 1)1432      quarter = pow(16, kap, p) * pow(kap * pow(cbm[kap], 2, p), -1, p) % p          # [w_kap]1433      check_tow(f"p={p}: pi angular-exact, g^(2 kap) == -1", pow(g, 2*kap, p) == p - 1)1434      check_tow(f"p={p}: pi wall-trivial, terminus [w_m] == -2", wall == p - 2)1435      check_tow(f"p={p}: pi quarter invariant [w_kap] == -(a^2)^-1, p = a^2 + b^2", quarter == (-pow(a*a, -1, p)) % p)1436      check_tow(f"p={p}: pi blind set below the wall empty (p does not divide C(2n,n) for n <= m)", all(cbm[n] % p != 0 for n in range(1, m+1)))1437      f, K, A, D_mod, blind = 1, 1, 1, [1, 0], []1438      for k in range(1, p): f = f * k % p; K = (K + f) % p; A = (A + (-1)**k * f) % p1439      for n in range(2, p - 1): D_mod.append((n * (D_mod[-1] + D_mod[-2])) % p)1440      blind = [n for n in range(1, p - 1) if D_mod[n] == 0]1441      i = pow(g, p-1-kap, p)1442      exact = any(2 * r * PI_LO2 <= (p - 1) <= 2 * r * PI_HI2 for r in (i, p - i))1443      check_tow(f"p={p}: e angular-approximate, no exact calibration (2 pi lambda(i) != p - 1 on the bracket)", not exact)1444      check_tow(f"p={p}: e wall-wild, the terminus the Kurepa quantity K(p) != 0", K != 0)1445      check_tow(f"p={p}: e blind set Z_0(p) nonempty (n = 1 forced): {blind}", 1 in blind)1446      check_tow(f"p={p}: the group law broken at the wall, K(p) A(p) != 1", (K * A) % p != 1)1447  flush("tow", order=list("EPCOHWD"), details={"H": f"median H/sqrt(p) = {_MED[0]:.3f} (diagnostic)"})14481449_MED = [0.0]14501451# ------------------------------------------------------------------------------------------------------------1452# block kur — kurepa_wall.py: kurepa_wall.py — the Kurepa wall of e in the derangement form, for all odd primes p < 2.5e5 (block kur.K1).1453#1454# Compiles and runs kurepa_wall.c (one O(p) pass per prime in 128-bit modular arithmetic, checking1455# !(p-1) == K(p) (mod p) and K(p) != 0 (mod p)), and turns its summary line into the registry predicate: 22 043 primes1456# checked, zero failures of either condition. A C compiler (cc/gcc/clang) is required; Colab and Debian/macOS provide1457# one. Without a compiler the block runs the same pass in pure Python to the bound 2·10⁴ (2 261 primes) and records1458# the reduced bound in the check's detail.1459N_kur = 25000014601461def compile_c_kur():1462    for cc in ("cc", "gcc", "clang"):1463        if shutil.which(cc):1464            exe = os.path.join(HERE, "kurepa_wall")1465            r = subprocess.run([cc, "-O2", os.path.join(HERE, "kurepa_wall.c"), "-o", exe], capture_output=True, text=True)1466            if r.returncode == 0:1467                return exe1468            print(f"    {cc} failed: {r.stderr.strip()[:200]}")1469    return None14701471def python_pass_kur(N_kur):1472    comp = bytearray(N_kur + 1)1473    for i in range(2, int(N_kur**0.5) + 1):1474        if not comp[i]: comp[i*i::i] = b"\x01" * len(comp[i*i::i])1475    checked = kfail = cfail = 01476    for p in range(3, N_kur + 1, 2):1477        if comp[p]: continue1478        f = s = D = 1; sign = -11479        for k in range(1, p):1480            f = f * k % p1481            s = (s + f) % p1482            D = (D * k + sign) % p1483            sign = -sign1484        if s == 0: kfail += 11485        if D != s: cfail += 11486        checked += 11487    return checked, kfail, cfail14881489def block_kur():1490    """Block kur — !(p−1) ≡ K(p) and K(p) ≢ 0 for all 22 043 odd primes p < 2.5·10⁵ (kurepa_wall.c, compiled on the fly; a pure-python pass to a smaller bound without a compiler) (kur.K1)."""1491    t = time.time()1492    exe = compile_c_kur()1493    if exe:1494        out = subprocess.run([exe, str(N_kur)], capture_output=True, text=True).stdout1495        print("    " + out.strip().replace("\n", "\n    "))1496        m = re.search(r"primes_checked=(\d+) kurepa_vanishing=(\d+) wall_congruence_failures=(\d+)", out)1497        checked, kfail, cfail = (int(x) for x in m.groups()) if m else (0, 1, 1)1498        # 13:J6 (p13057)1499        family("kur", "K1")1500        chk(f"{checked} odd primes p < {N_kur}: exactly 22043", checked == 22043)1501        chk("K(p) != 0 (mod p) for every prime checked", kfail == 0)1502        chk("!(p-1) == K(p) (mod p) for every prime checked", cfail == 0)1503        flush("kur", details={"K1": f"kurepa_wall.c, N = {N_kur}, {time.time() - t:.0f} s"})1504    else:1505        n = 200001506        checked, kfail, cfail = python_pass_kur(n)1507        family("kur", "K1")1508        chk(f"{checked} odd primes p < {n}: exactly 2261", checked == 2261)1509        chk("K(p) != 0 (mod p) for every prime checked", kfail == 0)1510        chk("!(p-1) == K(p) (mod p) for every prime checked", cfail == 0)1511        flush("kur", details={"K1": f"NO C COMPILER: pure-Python pass to N_kur = {n} only ({checked} primes), {time.time() - t:.0f} s"})15121513# ------------------------------------------------------------------------------------------------------------1514# block frm — frame_invariants.py: frame_invariants.py — the three frame invariants of the two constants on the 4 783 shells p ≡ 1 (mod 4) below 10⁵1515# (block frm.I1; 13:I5).1516#1517# Compiles and runs frame_invariants.c (one O(p) pass per prime: K(p) mod p, q_p(4) mod p, and p = a² + b² with a ≡ 11518# (mod 4), b > 0 even) and decides, in exact arithmetic, the paper's stated figures: the shell count; the triples at1519# p = 13, 233, 30089 and at the Wieferich prime 1093 (q_1093(4) = 0); the 4×4×4 contingency table of (K/p, q_p(4)/p,1520# φ/π) against the quarter bins — the first two by rational comparison, the third by the sign of a and |a| against b —1521# and its chi-square as an exact rational; the marginal means and variances of K/p and q_p(4)/p by certified rational1522# bounds of width 10⁻¹², rounded as the paper states them. The moments of φ/π and of a/√p (the arcsine law) use float arctangents and are1523# printed as [approx] readings; the independence reading (chi-square 59.6 on 63 df) is the paper's [approx] statement1524# and decides nothing. A C compiler is required; without one the block runs the same pass in pure Python to the bound1525# 2·10⁴ and records the reduced bound.1526N_frm = 10000015271528def compile_c_frm():1529    for cc in ("cc", "gcc", "clang"):1530        if shutil.which(cc):1531            exe = os.path.join(HERE, "frame_invariants")1532            r = subprocess.run([cc, "-O2", os.path.join(HERE, "frame_invariants.c"), "-o", exe, "-lm"], capture_output=True, text=True)1533            if r.returncode == 0:1534                return exe1535            print(f"    {cc} failed: {r.stderr.strip()[:200]}")1536    return None15371538def python_pass_frm(N_frm):1539    comp = bytearray(N_frm + 1)1540    for i in range(2, int(N_frm**0.5) + 1):1541        if not comp[i]: comp[i*i::i] = b"\x01" * len(comp[i*i::i])1542    rows = []1543    for p in range(5, N_frm, 4):1544        if comp[p]: continue1545        f = s = 11546        for k in range(1, p):1547            f = f * k % p; s = (s + f) % p1548        q4 = (pow(4, p - 1, p * p) - 1) // p1549        a = b = 01550        for x in range(1, math.isqrt(p) + 1, 2):1551            y = math.isqrt(p - x * x)1552            if y * y == p - x * x: a, b = x, y; break1553        rows.append((p, s, q4, a if a % 4 == 1 else -a, b))1554    return rows15551556def parse(out):1557    return [tuple(int(t) for t in line.split()) for line in out.strip().splitlines() if line.strip()]15581559def quarter_bin(x):                       # x a Fraction in [0, 1): the index of its quarter1560    return min(3, int(4 * x))15611562def angle_bin(a, b):                      # φ = arg(a + b i) with b > 0: the quarter of φ/π, decided in integers1563    if a > 0: return 0 if b < a else 1    # φ < π/4  iff b < a ; else φ ∈ [π/4, π/2)  (a = b never: p odd)1564    return 2 if b > -a else 3             # φ ∈ [π/2, 3π/4) iff b > |a| ; else [3π/4, π)15651566S = 10 ** 1215671568def moments(xs):1569    """Certified bounds on the mean and the variance of the rationals num/den, as (lo, hi) pairs of Fractions:1570    each term is bracketed by floor(num*S/den)/S and its successor, so the mean lies in an interval of width 1/S and the1571    variance E[x^2] - E[x]^2 in an interval of width < 3/S; no float enters."""1572    n = len(xs)1573    lo1 = sum(num * S // den for num, den in xs); hi1 = lo1 + n1574    lo2 = sum(num * num * S // (den * den) for num, den in xs); hi2 = lo2 + n1575    m_lo, m_hi = Fr(lo1, n * S), Fr(hi1, n * S)1576    v_lo, v_hi = Fr(lo2, n * S) - m_hi ** 2, Fr(hi2, n * S) - m_lo ** 21577    return (m_lo, m_hi), (v_lo, v_hi)15781579def rounds_to(bounds, digits, target):1580    lo, hi = bounds1581    return round(lo, digits) == target and round(hi, digits) == target15821583def block_frm():1584    """Block frm — the triples (K(p)/p, q_p(4)/p, φ/π) on the 4 783 shells p ≡ 1 (mod 4) below 10⁵ (frame_invariants.c, compiled on the fly) (frm.I1)."""1585    t = time.time()1586    exe = compile_c_frm()1587    if exe:1588        r = subprocess.run([exe, str(N_frm)], capture_output=True, text=True)1589        rows, bound = parse(r.stdout), N_frm1590    else:1591        bound = 200001592        rows = python_pass_frm(bound)1593    n = len(rows)1594    byp = {p: (K, q, a, b) for p, K, q, a, b in rows}1595    # exact consistency of every triple1596    consistent = all(0 < K < p and 0 <= q < p and a % 4 == 1 and b % 2 == 0 and b > 0 and a * a + b * b == p for p, K, q, a, b in rows)1597    # the 4×4×4 table, exact1598    table = {}1599    for p, K, q, a, b in rows:1600        key = (quarter_bin(Fr(K, p)), quarter_bin(Fr(q, p)), angle_bin(a, b))1601        table[key] = table.get(key, 0) + 11602    E = Fr(n, 64)1603    chi = sum((Fr(table.get((i, j, k), 0)) - E) ** 2 / E for i in range(4) for j in range(4) for k in range(4))1604    (meanK, varK), (meanQ, varQ) = moments([(K, p) for p, K, q, a, b in rows]), moments([(q, p) for p, K, q, a, b in rows])1605    phi = [math.atan2(b, a) / math.pi for p, K, q, a, b in rows]1606    meanP = sum(phi) / n; varP = sum((x - meanP) ** 2 for x in phi) / n1607    cosv = [a / math.sqrt(p) for p, K, q, a, b in rows]1608    meanC = sum(cosv) / n; varC = sum((x - meanC) ** 2 for x in cosv) / n1609    print(f"    {n} shells p = 1 (mod 4) below {bound}; K/p: mean {float(meanK[0]):.4f} var {float(varK[0]):.4f}; q_p(4)/p: mean {float(meanQ[0]):.4f} var {float(varQ[0]):.4f} (certified rational bounds, width 1e-12)")1610    print(f"    phi/pi: mean {meanP:.4f} var {varP:.4f} [approx]; a/sqrt p: mean {meanC:.3f} var {varC:.3f} [approx] (arcsine law: 0, 1/2); uniform law: 1/2, 1/12 = 0.0833")1611    print(f"    4x4x4 table against the product of the quarter bins: chi-square {float(chi):.2f} on 63 df (exact rational {chi.numerator}/{chi.denominator}) [approx reading: consistent with independence]")1612    # 13:I5 (p13068)1613    family("frm", "I1")1614    if bound == N_frm:1615        chk("exactly 4783 shells p = 1 (mod 4) below 10^5", n == 4783)1616        chk("the triples at p = 13, 233, 30089 are (10, 6, -3+2i), (69, 14, 13+8i), (5260, 7666, -67+160i)",1617            byp.get(13) == (10, 6, -3, 2) and byp.get(233) == (69, 14, 13, 8) and byp.get(30089) == (5260, 7666, -67, 160))1618        chk("the Wieferich prime 1093 is the one shell with q_p(4) = 0 (q_p(2) = 0 gives q_p(4) = 2q_p(2) + p q_p(2)^2 = 0)",1619            byp.get(1093, (0, -1))[1] == 0 and sum(1 for p, K, q, a, b in rows if q == 0) == 1)1620        chk("means of K/p and q_p(4)/p round to 0.502 and 0.499, variances to 0.0828 and 0.0831 (certified rational bounds)",1621            rounds_to(meanK, 3, Fr(502, 1000)) and rounds_to(meanQ, 3, Fr(499, 1000)) and rounds_to(varK, 4, Fr(828, 10000)) and rounds_to(varQ, 4, Fr(831, 10000)))1622        chk("the 4x4x4 chi-square is an exact rational in [59.5, 59.7] (the paper's 59.6 on 63 df)", Fr(595, 10) <= chi <= Fr(597, 10))1623    else:1624        chk(f"NO C COMPILER: {n} shells p = 1 (mod 4) below {bound} (pure-Python pass)", n > 0)1625    chk("every triple consistent: 0 < K(p) < p, 0 <= q_p(4) < p, a = 1 (mod 4), b > 0 even, a^2 + b^2 = p", consistent)1626    chk("the 64 cells of the table sum to the shell count", sum(table.values()) == n)1627    flush("frm", details={"I1": f"frame_invariants.c, N = {bound}, {time.time() - t:.0f} s" if bound == N_frm else f"NO C COMPILER: pure-Python pass to N = {bound}, {time.time() - t:.0f} s"})16281629# ------------------------------------------------------------------------------------------------------------1630# block lab — the laboratory Carriers' own shells (12 October 2026, the Y1 and Y3 pushes). A prime of the totality is a1631# Subject's shell, p² < Ω (00:C4), so a certificate "over the shells of the totality" is read over the primes p < √Ω and is1632# decided, on a laboratory Carrier Ω ∈ {233, 30089, 2 408 561}, over that Carrier's own shells: the odd primes p ≤ ⌊√(Ω−1)⌋,1633# 5, 39, 243 of them. lab.Z1: Kurepa's certificate Z = #{shells : p | K(p)} is zero on each Carrier (every shell walked,1634# Σ(p−1) steps in all), and on the same shells the two range laws of Z2 (c) — the blind-set sizes1635# against Poisson(1) (F5) and the index condition against independence (H5). lab.Z2: the calibration minimum over the1636# Carrier's shells p ≡ 1 (mod 4), decided exactly with π as a Machin bracket, against the nearest-hit scale π/(2N) of H4; the1637# three frame invariants on the same shells, exact, their correlations an [approx] reading. Pure Python, under a second.1638LAB_OMEGA = (233, 30089, 2408561)16391640def lab_shells(Om):1641    """The odd primes p with p² < Ω: the Subject shells a Carrier of cardinality Ω holds."""1642    return [p for p in sieve_e(math.isqrt(Om - 1)) if p > 2]16431644def lab_kurepa(p):1645    """K(p) = Σ_{k<p} k! mod p, the walk of p − 1 steps."""1646    f = s = 11647    for k in range(1, p):1648        f = f * k % p; s = (s + f) % p1649    return s16501651def lab_blind_size(p):1652    """|Z₀(p)| = #{1 ≤ n ≤ p−2 : !n ≡ 0 (mod p)}, the derangement numbers !n = (n−1)(!(n−1) + !(n−2)) walked mod p."""1653    d0, d1, z = 1, 0, 01654    for n in range(1, p - 1):1655        if d1 == 0: z += 11656        d0, d1 = d1, n * (d1 + d0) % p1657    return z16581659def lab_phi(n):1660    r = m = n; d = 21661    while d * d <= m:1662        if m % d == 0:1663            r -= r // d1664            while m % d == 0: m //= d1665        d += 11666    return r - r // m if m > 1 else r16671668def lab_primitive_root(p):1669    fac, m, d = [], p - 1, 21670    while d * d <= m:1671        if m % d == 0:1672            fac.append(d)1673            while m % d == 0: m //= d1674        d += 11675    if m > 1: fac.append(m)1676    return next(g for g in range(2, p) if all(pow(g, (p - 1) // q, p) != 1 for q in fac))16771678def lab_quarter_turn(p):1679    """The quarter-turn i = a^{(p−1)/4} for the least quadratic non-residue a (p ≡ 1 mod 4), as e.A1 and e.N1 take it."""1680    a = 21681    while pow(a, (p - 1) // 2, p) != p - 1: a += 11682    return pow(a, (p - 1) // 4, p)16831684def lab_index_condition(p):1685    """H5 on one shell: gcd(s, p−1) ∈ {gcd(λ(i), p−1), gcd(λ(−i), p−1)} for s the index of −K(p)⁻¹, with the probability of the1686    coincidence for s uniform on Z/(p−1) (the independence model): (φ((p−1)/d₁) + φ((p−1)/d₂)) / (p−1)."""1687    Ep = (-pow(lab_kurepa(p), -1, p)) % p; g = lab_primitive_root(p)1688    ind, x = {}, 11689    for m in range(p - 1): ind[x] = m; x = x * g % p1690    r = lab_quarter_turn(p); d1, d2 = math.gcd(r, p - 1), math.gcd(p - r, p - 1)1691    prob = Fraction(lab_phi((p - 1) // d1) + (lab_phi((p - 1) // d2) if d2 != d1 else 0), p - 1)1692    return math.gcd(ind[Ep], p - 1) in (d1, d2), prob16931694def block_lab():1695    """Block lab — the certificates of block Z on the laboratory Carriers' own shells (lab.Z1, lab.Z2)."""1696    from collections import Counter1697    t = time.time()1698    shells = {Om: lab_shells(Om) for Om in LAB_OMEGA}1699    kur = {Om: {p: lab_kurepa(p) for p in shells[Om]} for Om in LAB_OMEGA}1700    Z = {Om: sum(1 for v in kur[Om].values() if v == 0) for Om in LAB_OMEGA}1701    cost = {Om: sum(p - 1 for p in shells[Om]) for Om in LAB_OMEGA}1702    sizes = {Om: Counter(lab_blind_size(p) for p in shells[Om]) for Om in LAB_OMEGA}1703    h5 = {Om: [(p, *lab_index_condition(p)) for p in shells[Om] if p % 4 == 1 and p > 5] for Om in LAB_OMEGA}1704    for Om in LAB_OMEGA:1705        sh = shells[Om]; mean = Fraction(sum(k * v for k, v in sizes[Om].items()), len(sh))1706        solv = sum(ok for _, ok, _ in h5[Om]); exp = sum(pr for _, _, pr in h5[Om]); sd = math.sqrt(sum(pr * (1 - pr) for _, _, pr in h5[Om]))1707        print(f"    Omega = {Om}: {len(sh)} shells (odd primes below {math.isqrt(Om - 1) + 1}, the largest {sh[-1]}); Kurepa's certificate Z = {Z[Om]} "1708              f"over Σ(p−1) = {cost[Om]} walk steps; |Z0(p)|: {dict(sorted(sizes[Om].items()))}, mean {float(mean):.3f} (Poisson(1): 2) [approx]; "1709              f"index condition solvable on {solv} of {len(h5[Om])} shells p = 1 (mod 4), p > 5, against {float(exp):.1f} ± {sd:.1f} under independence [approx]")1710    own = {Om: lab_kurepa(Om) for Om in LAB_OMEGA}        # the Carriers' own residues: shells of the physical totality1711    print(f"    the laboratory Omega walked as shells of the physical totality: K(233) = {own[233]}, K(30089) = {own[30089]}, K(2408561) = {own[2408561]} (mod Omega)")1712    # 13:Z1 (p13061)1713    family("lab", "Z1")1714    chk("the shells of Omega = 233, 30089, 2408561 are the 5, 39, 243 odd primes below 16, 174, 1552, the largest 13, 173, 1549",1715        [len(shells[Om]) for Om in LAB_OMEGA] == [5, 39, 243] and [shells[Om][-1] for Om in LAB_OMEGA] == [13, 173, 1549])1716    chk("K(p) != 0 (mod p) on every shell of each Carrier: Z = 0, 0, 0", all(Z[Om] == 0 for Om in LAB_OMEGA))1717    chk("the certificates cost 34, 3046, 172452 walk steps", [cost[Om] for Om in LAB_OMEGA] == [34, 3046, 172452])1718    chk("K(233) = 69, K(30089) = 5260, K(2408561) = 435994 (mod Omega)", (own[233], own[30089], own[2408561]) == (69, 5260, 435994))1719    chk("|Z0(p)| over the 243 shells of 2408561: 1:87, 2:93, 3:43, 4:15, 5:3, 6:1, 7:1 (mean 490/243); over the 39 of 30089: 1:19, 2:11, 3:8, 4:1 (mean 69/39)",1720        dict(sizes[2408561]) == {1: 87, 2: 93, 3: 43, 4: 15, 5: 3, 6: 1, 7: 1} and dict(sizes[30089]) == {1: 19, 2: 11, 3: 8, 4: 1})1721    chk("the index condition is solvable on 48 of the 116 shells p = 1 (mod 4), p > 5, of 2408561; on 7 of the 17 of 30089; on the one shell (13) of 233",1722        [(sum(ok for _, ok, _ in h5[Om]), len(h5[Om])) for Om in LAB_OMEGA] == [(1, 1), (7, 17), (48, 116)])1723    # the calibration minimum over the Carrier's shells p = 1 (mod 4), decided exactly: for each shell and each chirality1724    # delta = |2 pi r/(p-1) - 1| is bracketed with pi's Machin bracket; the least candidate is decided when its upper bound1725    # lies below every other candidate's lower bound1726    PI_LO, PI_HI = pi_bracket(80); PI_MID = (PI_LO + PI_HI) / 21727    argmin, delta, N4 = {}, {}, {}1728    for Om in LAB_OMEGA:1729        cands = []1730        for p in shells[Om]:1731            if p % 4 != 1: continue1732            r = lab_quarter_turn(p)1733            for rr in (r, p - r):1734                lo = min(abs(2 * PI_LO * rr / (p - 1) - 1), abs(2 * PI_HI * rr / (p - 1) - 1))1735                hi = max(abs(2 * PI_LO * rr / (p - 1) - 1), abs(2 * PI_HI * rr / (p - 1) - 1))1736                if (2 * PI_LO * rr / (p - 1) - 1) * (2 * PI_HI * rr / (p - 1) - 1) < 0: lo = Fraction(0)1737                cands.append((lo, hi, p, rr))1738        N4[Om] = len(cands) // 21739        best = min(cands, key=lambda c: c[1]); decided = all(best[1] < c[0] for c in cands if c is not best)1740        argmin[Om] = (best[2], best[3]) if decided else None; delta[Om] = abs(2 * PI_MID * best[3] / (best[2] - 1) - 1)1741        print(f"    Omega = {Om}: {N4[Om]} shells p = 1 (mod 4); the calibration minimum at (p, lambda) = {argmin[Om]}, delta = {float(delta[Om]):.4e}, "1742              f"the nearest-hit scale pi/(2N) = {math.pi / (2 * N4[Om]):.4e}, ratio {float(delta[Om]) / (math.pi / (2 * N4[Om])):.3f} [approx]")1743    # the three frame invariants on the 117 shells of 2 408 561 (as frm.I1 takes them: K(p) mod p, q_p(4) mod p, p = a^2 + b^2, a = 1 (mod 4), b > 0)1744    rows = []1745    for p in shells[2408561]:1746        if p % 4 != 1: continue1747        q4 = (pow(4, p - 1, p * p) - 1) // p1748        a, b = next((x, math.isqrt(p - x * x)) for x in range(1, math.isqrt(p) + 1, 2) if math.isqrt(p - x * x) ** 2 == p - x * x)1749        rows.append((p, kur[2408561][p], q4, a if a % 4 == 1 else -a, b))1750    byp = {p: (K, q, a, b) for p, K, q, a, b in rows}1751    n = len(rows); T = [(K / p, q / p, math.atan2(b, a) / math.pi) for p, K, q, a, b in rows]1752    m = [sum(t[k] for t in T) / n for k in range(3)]; v = [sum((t[k] - m[k]) ** 2 for t in T) / n for k in range(3)]1753    corr = lambda i, j: sum((t[i] - m[i]) * (t[j] - m[j]) for t in T) / n / math.sqrt(v[i] * v[j])1754    c = (corr(0, 1), corr(0, 2), corr(1, 2)); sig = 1 / math.sqrt(n)1755    print(f"    the frame invariants on the {n} shells p = 1 (mod 4) of 2408561: means {m[0]:.3f}, {m[1]:.3f}, {m[2]:.3f}; variances {v[0]:.4f}, {v[1]:.4f}, {v[2]:.4f} "1756          f"(uniform 1/12 = 0.0833); correlations (K,q) {c[0]:+.3f}, (K,phi) {c[1]:+.3f}, (q,phi) {c[2]:+.3f} at sigma = 1/sqrt(n) = {sig:.3f} [approx]")1757    # 13:Z2 (p13063)1758    family("lab", "Z2")1759    chk("the shells p = 1 (mod 4) of the three Carriers number 2, 18, 117", [N4[Om] for Om in LAB_OMEGA] == [2, 18, 117])1760    chk("the calibration minimum is decided (the least bracket below every other) at (p, lambda) = (13, 5), (37, 6), (397, 63)",1761        [argmin[Om] for Om in LAB_OMEGA] == [(13, 5), (37, 6), (397, 63)])1762    chk("delta = 1.618 (in (1.61, 1.62)), 4.720e-2 (in (0.0471, 0.0473)), 4.023e-4 (in (4.02e-4, 4.03e-4))",1763        Fraction(161, 100) < delta[233] < Fraction(162, 100) and Fraction(471, 10**4) < delta[30089] < Fraction(473, 10**4) and Fraction(402, 10**6) < delta[2408561] < Fraction(403, 10**6))1764    chk("every triple on the 117 shells consistent: 0 < K(p) < p, 0 <= q_p(4) < p, a = 1 (mod 4), b > 0 even, a^2 + b^2 = p",1765        all(0 < K < p and 0 <= q < p and a % 4 == 1 and b % 2 == 0 and b > 0 and a * a + b * b == p for p, K, q, a, b in rows))1766    chk("the triples at p = 13 and p = 397 are (10, 6, -3+2i) and (152, 383, -19+6i)", byp.get(13) == (10, 6, -3, 2) and byp.get(397) == (152, 383, -19, 6))1767    flush("lab", order=["Z1", "Z2"], details={"Z1": f"pure Python, {time.time() - t:.1f} s", "Z2": f"correlations {c[0]:+.3f}, {c[1]:+.3f}, {c[2]:+.3f} at sigma {sig:.3f} (the paper's [approx] reading)"})17681769if __name__ == "__main__":1770    import time1771    want = [a.lower() for a in sys.argv[1:]] or list(BLOCK)1772    bad = [b for b in want if b not in BLOCK]1773    if bad: sys.exit(f"no block {', '.join(bad)}: the blocks are {', '.join(BLOCK)}")1774    t0 = time.time()1775    for b in want:1776        t = time.time(); print(f"— block {b}"); _run_block(b); print(f"    [block {b}: {time.time() - t:.1f} s]")1777    ok = summary(write=(want == list(BLOCK)))1778    print(f"{len(RESULTS)} family checks, {len(MICRO)} exact micro-checks, {time.time() - t0:.1f} s" + ("; results.json written" if want == list(BLOCK) else ""))1779    sys.exit(0 if ok else 1)1780