algebra.py — py

1"""2algebra.py — the validation package of "Relativistic Algebra over Finite Ring Continuum" (Akhtman, Axioms 2025, 14, 636;3doi 10.3390/axioms14080636), the paper 1-algebra of the FRC corpus (finite-ring-space/src/1-algebra), added with the paper's4predicate ledger (Appendix A, 16 September 2026; one script since 24 September 2026, the three block scripts merged).5========================================================================================================================67One script, three blocks, eighteen checks, standard library only. Each check names the predicate(s) of the paper's ledger8it witnesses (LEDGER below; predicates cited as 1:XN) under a `# 1:XN (<key>)` marker, and the ledger's source column9links the marker in return (finitering.space/py/algebra.html#<key>). The paper \\label(s) a check decides are10in the block banners. Eight master-ledger predicates of the corpus cite this paper (00:A8, B1, B8, B11, C9, C10, C25, D13); the11predicates of the paper ledger that the master carries are listed in the site generator.1213    python3 algebra.py            every block, results.json written; exit 1 if a check fails (≈ 25 s)14    python3 algebra.py C          one block (A, B or C); no results.json15    from frc_1_algebra import predicate; predicate("1:B2")     one predicate: its block runs once per session1617Blocks:  A  the shell, its frame and the orbital complex          EXACT          (1:B2–B4, C2, C4)18         B  the framed numbers, the charts and the horizon         EXACT / CHART  (1:D2, D4, D6, E2, Z1, Z3; B6 decides no predicate)19         C  the conjecture of the conclusion, clause by clause     EXACT / CHART  (1:E3, D7, E4–E8)2021Everything the blocks share:22  * the shell datum: p = 4κ+1 prime, the primitive roots of F_p, the oriented quarter-turn i = −g^κ;23  * exact arithmetic in F_p (python int, reduced mod p) and in Q (fractions.Fraction);24  * the PASS/FAIL registry every block reports into (results.json).2526Kinds: EXACT checks are integer-pinned computations in F_p or exact rationals (a pass is a proof on the27tested instances); CHART checks decide a [chart] predicate — a statement about the reading of the shell against28Q or R — in exact rationals.29"""30import os, json, sys, math, random, itertools31from fractions import Fraction32from decimal import Decimal, getcontext3334SCRIPT = os.path.splitext(os.path.basename(__file__))[0]        # "algebra": the one script, the name results.json and the site pages carry3536# The shells: p = 4κ+1 prime; the controls: p ≡ 3 (mod 4).37SHELLS = [5, 13, 17, 29, 37, 41, 173]38CONTROLS = [7, 11, 19, 23]3940# ----------------------------------------------------------------------------- registry41RESULTS = []4243# The paper's predicate ledger (Appendix A, predicates cited as 1:XN): the predicate(s) each check witnesses.44LEDGER = {45    "A1": "1:B2", "A2": "1:B3", "A3": "1:B4", "A4": "1:C2", "A5": "1:C4",46    "B1": "1:D2", "B2": "1:D4", "B3": "1:D6", "B4": "1:E2", "B5": "1:Z1", "B6": "", "B7": "1:Z3",   # B6 (the Euclidean step count) decides no predicate of the ledger47    "C1": "1:E3", "C2": "1:D7", "C3": "1:E4", "C4": "1:E5", "C5": "1:E6", "C6": "1:E7", "C7": "1:E8",48}4950BLOCK = {"A": "the shell, its frame and the orbital complex",             # check-id prefix -> the block (the function block_<letter> below)51         "B": "the framed numbers, the charts and the horizon",52         "C": "the conjecture of the conclusion, clause by clause"}5354# the deciding check of each witnessed predicate: the one whose verdict decides the predicate's statement (the other checks that55# touch it are corroboration, listed by predicate() from the records)56PREDICATES = {57    "1:B2": "A1", "1:B3": "A2", "1:B4": "A3", "1:C2": "A4", "1:C4": "A5",58    "1:D2": "B1", "1:D4": "B2", "1:E2": "B4", "1:Z1": "B5", "1:Z3": "B7",         # D6: B3 witnesses it in LEDGER but carries no marker (as before the merge)59    "1:E3": "C1", "1:D7": "C2", "1:E4": "C3", "1:E5": "C4", "1:E6": "C5", "1:E7": "C6", "1:E8": "C7",60}61_RAN = set()                                            # blocks already run in this session (predicate() runs each once)6263def check(pid, label, ok, detail="", kind="EXACT"):64    """Record one predicate check. pid = package check id; LEDGER[pid] = the paper statement(s) decided."""65    ok = bool(ok)66    rows = LEDGER.get(pid, "")67    RESULTS.append({"id": pid, "rows": rows, "block": pid[0], "script": SCRIPT, "label": label, "ok": ok, "detail": detail, "kind": kind})68    print(f"  [{'PASS' if ok else 'FAIL'}] {pid:4s} {kind:6s} [{rows}] {label}" + (f"  --  {detail}" if detail else ""))69    return ok7071def summary(write=True):72    n_ok = sum(r["ok"] for r in RESULTS)73    print(f"\nSUMMARY: {n_ok}/{len(RESULTS)} checks passed" + ("" if n_ok == len(RESULTS) else "  <-- FAILURES"))74    if write:75        with open("results.json", "w") as f:76            json.dump(RESULTS, f, indent=1)77    return n_ok == len(RESULTS)7879def markers():80    """predicate label -> (script file, line) of its marker `# <paper>:<label> (<key>)`: the line of the check that decides it."""81    import re82    out = {}83    for i, line in enumerate(open(os.path.abspath(__file__), encoding="utf-8"), 1):84        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)85        if m:86            for lab in re.findall(r"\d+:[A-Z]+\d+[a-z]?", m.group(1)): out.setdefault(lab, (SCRIPT + ".py", i))87    return out8889def _run_block(letter):90    if letter not in _RAN:91        globals()[f"block_{letter}"](); _RAN.add(letter)9293def predicate(label, lines=14):94    """Verify one ledger predicate: run the block of the check that decides it (once per session), print that check's source95    (from its marker) and every record that cites the predicate, and return True iff all pass."""96    pid = PREDICATES.get(label)97    if pid is None:98        print(f"{label}: no python witness (see the predicate's Lean witness or its source)"); return None99    citing = {i[0] for i, rows in LEDGER.items() if label in [t.strip() for t in rows.split(",")]}100    for b in sorted({pid[0]} | citing): _run_block(b)          # the deciding block and every block whose checks cite the predicate101    mk = markers().get(label)102    if mk:103        src = open(os.path.abspath(__file__), encoding="utf-8").read().split("\n")104        print(f"— {mk[0]}:{mk[1]} (the check that decides {label}: {pid})")105        for j in range(mk[1] - 1, min(mk[1] - 1 + lines, len(src))): print(f"{j + 1:5d}  {src[j]}")106    recs = [r for r in RESULTS if label in [t.strip() for t in r["rows"].split(",")]]107    ok = all(r["ok"] for r in recs)108    for r in recs:109        role = "(deciding)" if r["id"] == pid else "(corroborating)"110        print(f"  [{'PASS' if r['ok'] else 'FAIL'}] {r['id']:4s} {role:16s} {r['detail'][:150]}")111    print(f"{label}: {'VERIFIED' if ok and recs else 'FAILED'} — {len(recs)} record(s)")112    return ok113114def verify_all():115    """Run every block (those already run in this session are not re-run) and print the summary; True iff every check passed."""116    for b in sorted(BLOCK): _run_block(b)117    return summary(write=False)118119# ----------------------------------------------------------------------------- the shell datum120def is_prime(n):121    if n < 2: return False122    d = 2123    while d * d <= n:124        if n % d == 0: return False125        d += 1126    return True127128def prime_factors(n):129    f, d = set(), 2130    while d * d <= n:131        while n % d == 0:132            f.add(d); n //= d133        d += 1134    if n > 1: f.add(n)135    return f136137def primitive_roots(p):138    n = p - 1139    fac = prime_factors(n)140    return [g for g in range(2, p) if all(pow(g, n // q, p) != 1 for q in fac)]141142def kappa(p):143    assert p % 4 == 1 and is_prime(p), p144    return (p - 1) // 4145146def quarter_turn(p, g):147    """i = −g^κ, the oriented quarter-turn of the frame (τ; 0, 1, g)."""148    return (-pow(g, kappa(p), p)) % p149150def sqrt_neg_one(p):151    return sorted(x for x in range(1, p) if (x * x) % p == p - 1)152153# ------------------------------------------------------------------------------------------------------------154# block A: the shell, its frame and the orbital complex (EXACT, integer-pinned)155# Paper statements decided (Sections 2–3; ledger predicates 1:B2–B4, 1:C2, 1:C4):156#157#   A1  thm:symmetric-completeness   the fourth roots of unity are Q4 = {1, i, −1, −i}; under the Klein four-group158#       (1:B2)                       ⟨x ↦ −x, x ↦ x⁻¹⟩, Q4 is the union of the two size-2 orbits {±1}, {±i}, and159#                                    F_p^× \ Q4 splits into exactly κ−1 orbits of size 4; for p ≡ 3 (mod 4) no i160#   A2  (1:B3)                       i = −g^κ satisfies i² = −1 for every primitive root g; {−g^κ, g^κ} are the two161#                                    square roots of −1; on F_13 with g = 2, i = 5162#   A3  def:framed-field,            φ_{a,b}(x) = a + bx is a ring isomorphism (F_p, +, ·) → (F_p, ⊕, ⊗) on every163#       lem:affine-invariance (1:B4) frame (a, b), b ≠ 0; the unit of the relabelled field is a + b, and b is the164#                                    unit only when a = 0165#   A4  prop:frame-group (1:C2)      ⟨T_a, S_m⟩ = Aff(F_p) of order p(p−1), simply transitive on the frames166#   A5  def:orbit-sphere (1:C4)      M_n(a) = M_{n+2κ}(−a), L_a(m) = L_{−a}(m+2κ); p−1 distinct meridian lists in167#                                    2κ great circles; 2κ latitude pairs; vertex count (p−1)²/2 + 1168#169# Shells: p ∈ {5, 13, 17, 29, 37, 41, 173} for A1–A2; p ∈ {13, 17} for A3–A5 (exhaustive over frames and pairs).170def klein_orbits(p):171    seen, orbits = set(), []172    for x in range(1, p):173        if x in seen: continue174        inv = pow(x, -1, p)175        orb = {x, (-x) % p, inv, (-inv) % p}176        seen |= orb; orbits.append(orb)177    return orbits178179def block_A():180    """Block A — the shell, its frame and the orbital complex (EXACT, integer-pinned): A1–A5."""181    print("block A — the shell, its frame and the orbital complex")182183    # A1 — symmetry completeness184    ok, det = True, []185    for p in SHELLS:186        k = kappa(p); orbits = klein_orbits(p)187        Q4 = {1, p - 1} | set(sqrt_neg_one(p))188        small = [o for o in orbits if len(o) < 4]189        big = [o for o in orbits if len(o) == 4]190        ok &= (len(Q4) == 4 and set().union(*small) == Q4 and sorted(len(o) for o in small) == [2, 2])191        ok &= (len(big) == k - 1 and len(big) + len(small) == k + 1)192        ok &= all(pow(x, 4, p) == 1 for x in Q4) and sum(1 for x in range(1, p) if pow(x, 4, p) == 1) == 4193        det.append(f"p={p}: Q4={sorted(Q4)}, {len(big)} orbits of size 4 = κ−1")194    for p in CONTROLS:195        orbits = klein_orbits(p)196        ok &= (sqrt_neg_one(p) == [] and sorted(len(o) for o in orbits) == [2] + [4] * ((p - 3) // 4))197    # 1:B2 (p01004)198    check("A1", "Q4 = {1,i,−1,−i} the unique order-4 subgroup, two Klein orbits of size 2, κ−1 orbits of size 4; no i for p ≡ 3 (mod 4)",199          ok, "; ".join(det) + f"; controls p ∈ {CONTROLS}")200201    # A2 — the oriented quarter-turn202    ok, det = True, []203    for p in SHELLS:204        k = kappa(p); roots = set(sqrt_neg_one(p)); n = 0205        for g in primitive_roots(p):206            i = quarter_turn(p, g)207            ok &= ((i * i) % p == p - 1) and ({i, pow(g, k, p)} == roots) and (pow(g, 2 * k, p) == p - 1); n += 1208        det.append(f"p={p}: {n} primitive roots")209    ok &= (quarter_turn(13, 2) == 5)210    # 1:B3 (p01005)211    check("A2", "i = −g^κ, i² = −1 for every primitive root; {−g^κ, g^κ} the two square roots of −1; g^{2κ} = −1; F_13, g = 2: i = 5",212          ok, "; ".join(det))213214    # A3 — the affine frame: ring isomorphism, the unit a + b215    ok, n = True, 0216    for p in [13, 17]:217        for a in range(p):218            for b in range(1, p):219                inv = lambda X: ((X - a) * pow(b, -1, p)) % p220                phi = lambda x: (a + b * x) % p221                oplus = lambda X, Y: (a + b * (inv(X) + inv(Y))) % p222                otimes = lambda X, Y: (a + b * (inv(X) * inv(Y))) % p223                ok &= all(phi((x + y) % p) == oplus(phi(x), phi(y)) and phi((x * y) % p) == otimes(phi(x), phi(y))224                          for x in range(p) for y in range(p))225                ok &= (len({phi(x) for x in range(p)}) == p)226                ok &= all(otimes((a + b) % p, y) == y for y in range(p))227                ok &= ((a == 0) == all(otimes(b, y) == y for y in range(p)))228                n += 1229    # 1:B4 (p01006)230    check("A3", "φ_{a,b} = a + bx a ring isomorphism onto (F_p, ⊕, ⊗) on every frame; unit a + b; b the unit iff a = 0",231          ok, f"{n} frames on p ∈ {{13, 17}}, all pairs (x, y)")232233    # A4 — the frame group234    ok, det = True, []235    for p in [13, 17]:236        frames = [(a, b) for a in range(p) for b in range(1, p)]237        comp = lambda f, g: ((f[0] + f[1] * g[0]) % p, (f[1] * g[1]) % p)     # f∘g as affine maps x ↦ a + b x238        ok &= (len(frames) == p * (p - 1))239        ok &= all(comp(f, g) in set(frames) for f in frames for g in frames)  # closure240        # simple transitivity: for each (φ, φ') exactly one ψ with ψ∘φ = φ'241        cnt = 0242        for f in frames:243            for f2 in frames:244                sols = [psi for psi in frames if comp(psi, f) == f2]245                ok &= (len(sols) == 1); cnt += 1246        det.append(f"p={p}: |Aff| = {len(frames)} = p(p−1), {cnt} frame pairs, one carrier each")247    # 1:C2 (p01009)248    check("A4", "⟨T_a, S_m⟩ = Aff(F_p), order p(p−1), simply transitive on the frames (a, b)", ok, "; ".join(det))249250    # A5 — the orbital complex: involutions and counts251    ok, det = True, []252    for p in [13, 17]:253        k = kappa(p); g = primitive_roots(p)[0]; n2 = 2 * k254        M = lambda n, a: (a * pow(g, n, p)) % p255        L = lambda a, m: (a * pow(g, m, p)) % p256        ok &= all(M(n, a) == M(n + n2, (-a) % p) for n in range(p - 1) for a in range(p))257        ok &= all(L(a, m) == L((-a) % p, m + n2) for a in range(1, p) for m in range(p - 1))258        mer = {tuple(M(n, a) for a in range(p)) for n in range(p - 1)}259        circles = {frozenset({tuple(M(n, a) for a in range(p)), tuple(M(n + n2, a) for a in range(p))}) for n in range(p - 1)}260        latpairs = {frozenset({tuple(L(a, m) for m in range(p - 1)), tuple(L((-a) % p, m) for m in range(p - 1))}) for a in range(1, p)}261        ok &= (len(mer) == p - 1 and len(circles) == 2 * k and len(latpairs) == 2 * k)262        V = (p - 1) * (2 * k) + 1263        ok &= (V == (p - 1) ** 2 // 2 + 1)264        det.append(f"p={p}: {len(mer)} meridian lists, {len(circles)} great circles, {len(latpairs)} latitude pairs, |V| = {V}")265    # 1:C4 (p01011)266    check("A5", "M_n(a) = M_{n+2κ}(−a), L_a(m) = L_{−a}(m+2κ); p−1 meridian lists in 2κ circles; 2κ latitude pairs; |V| = (p−1)²/2 + 1",267          ok, "; ".join(det))268269# ------------------------------------------------------------------------------------------------------------270# block B: the framed numbers, the charts and the horizon (EXACT; B3 CHART in exact rationals)271# Paper statements decided (Sections 4–5; ledger predicates 1:D2, 1:D4, 1:D6, 1:E2, 1:Z1, 1:Z3; B6 decides no predicate):272#273#   B1  def:integers (1:D2)          the window W_H = {|z| ≤ H}: z ↦ z mod p injective iff 2H < p; sums read back274#                                    iff 4H < p; products read back whenever 2H² < p (sufficient; the product set is275#                                    sparse, so the bound is not sharp) — swept over every H276#   B2  thm:scale-periodicity (1:D4) the residue grids G_n = (x g^{−n})_x are (p−1)-periodic in n; the rational277#                                    grids (x/g^n) with g lifted to Z are not278#   B3  prop:r-rationals, thm:approx the chart of the grid [chart]: every r with |r| ≤ H g^{−n} lies within 1/(2gⁿ)279#       (1:D6)                       of some x/gⁿ, x ∈ W_H (exact rationals); at (13, 2) no grid point of step ≤ 1/8280#                                    lies within 1/16 of 33/10 — range and resolution trade off at fixed window281#   B4  prop:Cp-field (1:E2)         F_p[X]/(X²+1) has zero divisors on the shell ((u+X)(u−X) = 0, factors nonzero);282#                                    X²+1 has no root, hence the quotient is a field, exactly for p ≡ 3 (mod 4)283#   B5  thm:no-south-pole (1:Z1)     2s = 0 ⇒ s = 0 on every odd prime; 2·(2κ+1) = 1: the half-turn 2⁻¹ = 2κ+1284#   B6  (no predicate)               the Euclidean step count against the bound k ≤ ⌊log₂ p⌋+1 (lem:euclid-bound):285#                                    p = 59 (55, 34) needs 7 > 6; p = 1009 (987, 610) needs 13 > 10; first excess at 59286#   B7  thm:no-equator (1:Z3)        4s = 0 ⇒ s = 0; a fraction a/b (b ≢ 0) at a quarter position of the additive cycle287#                                    (p/4, p/2, 3p/4 mod p) has value 0, and none with a, b in a window 2H < p is there;288#                                    the equator p/4 = κ + 1/4 lies between L_κ and L_{κ+1}; κ ≡ −4⁻¹, κ + 1 ≡ 3·4⁻¹289def block_B():290    """Block B — the framed numbers, the charts and the horizon (EXACT; B3 CHART in exact rationals): B1–B6."""291    print("block B — the framed numbers, the charts and the horizon")292293    # B1 — the window law, the three thresholds as iff over every H294    ok, det = True, []295    for p in [13, 17, 29]:296        for H in range(1, p):297            W = range(-H, H + 1)298            inj = len({z % p for z in W}) == len(W)299            sums = len({(x + y) % p for x in W for y in W}) == len({x + y for x in W for y in W})300            prods = len({(x * y) % p for x in W for y in W}) == len({x * y for x in W for y in W})301            ok &= (inj == (2 * H < p)) and (sums == (4 * H < p)) and ((2 * H * H < p) <= prods)   # products: sufficient, not necessary (the product set is sparse)302        det.append(f"p={p}: H = 1..{p-1} swept")303    # 1:D2 (p01013)304    check("B1", "z ↦ z mod p injective on W_H iff 2H < p; sums read back iff 4H < p; products read back whenever 2H² < p", ok, "; ".join(det))305306    # B2 — scale-periodicity: in the field yes, in Q no307    ok, det = True, []308    for p, g in [(13, 11), (13, 2), (17, 3), (29, 2)]:309        assert g in primitive_roots(p)310        grid = lambda n: tuple((x * pow(g, -n, p)) % p for x in range(p))311        ok &= all(grid(n) == grid(n + p - 1) for n in range(2 * (p - 1)))312        ok &= all(grid(n) != grid(n + d) for n in range(p - 1) for d in range(1, p - 1))   # exact period p−1313        qgrid = lambda n: tuple(Fraction(x, g ** n) for x in range(p))314        ok &= all(qgrid(n) != qgrid(n + p - 1) for n in range(3))315        det.append(f"({p},{g}): period {p-1} in F_p, none in Q")316    # 1:D4 (p01015)317    check("B2", "G_n = (x g^{−n})_x is (p−1)-periodic in the field, with exact period p−1; the rational grids x/gⁿ are not periodic",318          ok, "; ".join(det))319320    # B3 — the chart of the grid: the trade-off holds; its range at (13, 2)321    ok, det = True, []322    p, g, H = 13, 2, 12323    for n in range(0, 8):324        step = Fraction(1, g ** n)325        for num in range(-40, 41):                                    # r = num/8 · H g^{−n}, a sample within the range326            r = Fraction(num, 40) * H * step327            x = round(r / step)                                        # x = round(r gⁿ), |x| ≤ H328            ok &= (abs(x) <= H) and (abs(r - x * step) <= step / 2)329    r, tol = Fraction(33, 10), Fraction(1, 16)330    near = [(x, n) for n in range(3, 13) for x in range(p) if abs(Fraction(x, g ** n) - r) < tol]331    ok &= (near == []) and all(Fraction(x, g ** n) <= Fraction(3, 2) for n in range(3, 13) for x in range(p))332    best = min((abs(Fraction(x, g ** n) - r), x, n) for n in range(13) for x in range(p))333    ok &= (best[1:] == (7, 1) and best[0] == Fraction(1, 5))334    det.append(f"trade-off on {8*81} sampled reals at (13,2), H=12; range: no x/2^n (n ≥ 3, x < 13) within 1/16 of 33/10, best 7/2 at error 1/5")335    check("B3", "|r| ≤ H g^{−n} ⇒ some x/gⁿ, x ∈ W_H, within 1/(2gⁿ); at (13, 2) no point of step ≤ 1/8 within 1/16 of r = 33/10: range and resolution trade off at fixed window",336          ok, "; ".join(det), kind="CHART")337338    # B4 — the complex chart is not an extension: zero divisors iff −1 is a square339    ok, det = True, []340    for p in [5, 13, 17, 29, 37]:341        u = sqrt_neg_one(p)[0]342        # in F_p[X]/(X²+1): elements a + bX, product (a+bX)(c+dX) = (ac − bd) + (ad + bc)X343        mul = lambda A, B: ((A[0] * B[0] - A[1] * B[1]) % p, (A[0] * B[1] + A[1] * B[0]) % p)344        ok &= (mul((u, 1), ((-u) % p, 1)) == (0, 0)) and (u, 1) != (0, 0) and ((-u) % p, 1) != (0, 0)345        det.append(f"p={p}: (u+X)(u−X)=0 with u={u}")346    for p in CONTROLS:347        ok &= (sqrt_neg_one(p) == [])                                    # X²+1 irreducible: the quotient is the field F_{p²}348        mul = lambda A, B: ((A[0] * B[0] - A[1] * B[1]) % p, (A[0] * B[1] + A[1] * B[0]) % p)349        elems = [(a, b) for a in range(p) for b in range(p) if (a, b) != (0, 0)]350        ok &= all(mul(A, B) != (0, 0) for A in elems for B in elems)    # no zero divisors351    # 1:E2 (p01019)352    check("B4", "F_p[X]/(X²+1) has zero divisors on the shell (p ≡ 1 mod 4); it is a field exactly when p ≡ 3 (mod 4)",353          ok, "; ".join(det) + f"; fields at p ∈ {CONTROLS} (no zero divisors, exhaustive)")354355    # B5 — no element of additive order two; the half-turn356    ok, det = True, []357    for p in [q for q in range(3, 200) if is_prime(q)]:358        ok &= ([s for s in range(p) if (2 * s) % p == 0] == [0])359        if p % 4 == 1:360            k = kappa(p)361            ok &= ((2 * (2 * k + 1)) % p == 1) and (2 * k + 1 == (p + 1) // 2)362    # 1:Z1 (p01020)363    check("B5", "2s = 0 ⇒ s = 0 for every odd prime < 200; on the shell 2⁻¹ = 2κ+1 = (p+1)/2, the residue past the antipode", ok, "primes 3..199")364365    # B6 — the Euclidean step count against the bound ⌊log₂ p⌋+1 (recorded under V1)366    def steps(a, b):367        k = 0368        while b:369            a, b = b, a % b370            if b: k += 1371        return k372    fib = [1, 1]373    while fib[-1] < 2000: fib.append(fib[-1] + fib[-2])374    first = None375    for p in [q for q in range(5, 2000) if is_prime(q)]:376        worst = max(steps(fib[i + 1], fib[i]) for i in range(1, len(fib) - 1) if fib[i + 1] < p)377        if worst > p.bit_length():378            first = p; break379    ok = (steps(55, 34) == 7 and (59).bit_length() == 6 and steps(987, 610) == 13 and (1009).bit_length() == 10 and first == 59)380    ok &= (fib[3] == 3 and fib[3] < 2 ** 2)                               # F_4 = 3 < 2^{4−2}: the Fibonacci bound 2^{k−2} does not hold at k = 4381    check("B6", "the Euclidean step count against ⌊log₂p⌋+1: (55,34) needs 7 > ⌊log₂59⌋+1 = 6; (987,610) needs 13 > 10; first excess at p = 59; F_4 = 3 < 4",382          ok, "a ≥ b convention, consecutive Fibonacci pairs below p")383384    # B7 — no latitude on the equator: the quarter positions of the additive cycle carry only fractions of value 0385    ok, det = True, []386    for p in SHELLS:387        k = kappa(p); inv4 = pow(4, p - 2, p)388        ok &= ([s for s in range(p) if (4 * s) % p == 0] == [0])                                  # no element of additive order four389        ok &= (Fraction(p, 4) == k + Fraction(1, 4)) and (k % p == (-inv4) % p) and ((k + 1) % p == (3 * inv4) % p)390        H = (p - 1) // 2                                                                           # the widest window, 2H < p391        hits = 0392        for b in range(1, 3 * p):393            if b % p == 0: continue394            for a in range(-3 * p, 3 * p + 1):395                q4 = 4 * Fraction(a, b)396                if q4.denominator == 1 and q4 % p == 0 and (q4 // p) % 4 != 0:                     # a/b ≡ p/4, p/2 or 3p/4 (mod p)397                    hits += 1398                    ok &= (a % p == 0) and ((a * pow(b, p - 2, p)) % p == 0)                     # the fraction's value is 0399                    ok &= not (abs(a) <= H and abs(b) <= H)                                       # and no window fraction sits there400        ok &= hits > 0401        det.append(f"p={p}: {hits} quarter-position fractions a/b, |a| ≤ 3p, b < 3p, all of value 0")402    # 1:Z3 (p01034)403    check("B7", "4s = 0 ⇒ s = 0; every fraction at a quarter position of the additive cycle has value 0 and none of the window is there; the equator κ+1/4 between L_κ and L_{κ+1}, κ ≡ −4⁻¹, κ+1 ≡ 3·4⁻¹",404          ok, "; ".join(det))405406# ------------------------------------------------------------------------------------------------------------407# block C: the conjecture of the conclusion, clause by clause (EXACT; C2–C4 CHART)408# The conclusion conjectures that the finite substrate supports polynomial equation solving, limit-like409# approximation and ε-approximation of continuous symmetries.  Predicates 1:E3, D7 and E4–E8 decide it:410#411#   C1  (1:E3)  solving: f ∈ F_p[X] has a root in F_p iff gcd(f, X^p − X) ≠ 1, and deg gcd = the number of412#               distinct roots — every monic polynomial of degree ≤ 3 over F_13 (2197) and F_17 (4913), by413#               brute-force roots against the Euclidean gcd            [exact]414#   C2  (1:D7)  limit-like approximation across shells: for r ∈ {π, e, √2, 33/10, 1/3} and ε = 10^−k415#               (k = 2..8) the construction n = ⌈log₂(1/ε)⌉, x = round(r·2ⁿ), p = the first prime ≡ 1 (mod 4)416#               beyond 2|x|+1 gives |r − x/2ⁿ| < ε with |x| ≤ 2κ      [chart: r read as a rational to 60 digits]417#   C3  (1:E4)  the circle net: for N = p − 1 on every shell the rounding k(θ) = ⌊Nθ/2π + ½⌋ has angle error418#               ≤ π/N, chord error |e^{2πik/N} − e^{iθ}| ≤ π/N, and group-law defect |k(θ₁)+k(θ₂)−k(θ₁+θ₂)| ≤ 1,419#               on a grid of 4N angles plus 20000 random pairs           [chart]420#   C4  (1:E5)  the SO(3) obstruction: the covering radii of the finite rotation groups in the rotation-angle421#               metric — tetrahedral π/2, octahedral arccos((2√2−1)/4) ≈ 62.80°, icosahedral422#               ε₀ = arccos((3√5−1)/8) ≈ 44.48° (the deep hole of the 600-cell, cos = φ²/2√2), cyclic and423#               dihedral ≥ π/2 (they lie in an O(2)) — by exhaustive group closure and a sampled maximin;424#               no finite subgroup of SO(3) is an ε-net for ε < ε₀      [chart; the list of groups is A3, Klein]425#   C5  (1:E6)  the window resolves SO(3): the framed quaternions W_H⁴ = {q ∈ Z⁴ : |q_i| ≤ H}, normalised, are an426#               ε-net of SO(3) with ε ≤ 2·arcsin(1/H) (round H·s to the lattice: |Hs − q| ≤ 1, so the angle is427#               ≤ arcsin(1/H)); measured ε(1) ≈ 60.8°, ε(2) ≈ 41.0°, ε(3) ≈ 30.0°; 2·arcsin(1/3) = 38.9° < ε₀, so on428#               the shell p = 73 the window H = 3 out-resolves every finite subgroup; products of window quaternions429#               have entries in W_{4H²} and read back exactly from the shell when 8H² < p (20000 random pairs on430#               F_73, H = 3)                                              [chart; the read-back is exact]431#   C6  (1:E7)  the Lie-algebra layer, one axis: every window rotation ρ_q, q = a + v, equals (aI − S_v)⁻¹(aI + S_v) — the432#               Cayley step of ηS_v at α = η/(νa), computed in F_p[η] — and I + 2N(v)⁻¹S_v² at a = 0; on every axis the433#               steps, the identity and the half-turn are the powers of one element of order p − (N(v)|p), every434#               iterate stays among them, and at most 2H² + 2H are window rotations (p = 13, 17, 37, 73, 101) [exact]435#   C7  (1:E8)  two axes and range: ρ_(1+u)ρ_(1+w) = ρ_(1+u∘w), u∘w = (u + w + u×w)/(1 − u·w); the half-turn about436#               u + w + u×w at u·w = 1; uw − wu = 2u×w (F_13, F_29, F_73, F_233); the entries of q^m are at most (2H)^m and437#               q^m reads back from the shell iff 2‖q^m‖∞ < p (H ≤ 3, m ≤ 6; p = 13, 101, 1009)                   [exact]438439# ---------------------------------------------------------------- polynomial arithmetic over F_p (lists, low degree first)440def p_trim(a):441    while a and a[-1] == 0: a.pop()442    return a443def p_mod(a, b, p):444    a = a[:]; db = len(b) - 1; inv = pow(b[-1], -1, p)445    while len(a) - 1 >= db and a:446        c = (a[-1] * inv) % p; s = len(a) - 1 - db447        for i, bi in enumerate(b): a[s + i] = (a[s + i] - c * bi) % p448        p_trim(a)449    return a450def p_gcd(a, b, p):451    a, b = p_trim(a[:]), p_trim(b[:])452    while b: a, b = b, p_mod(a, b, p)453    if a: inv = pow(a[-1], -1, p); a = [(c * inv) % p for c in a]454    return a455def p_eval(a, x, p):456    r = 0457    for c in reversed(a): r = (r * x + c) % p458    return r459460def block_C():461    """Block C — the conjecture of the conclusion, clause by clause (EXACT; C2–C5 CHART): C1–C5."""462    print("block C — the conjecture of the conclusion, clause by clause")463464    # C1 — solving: roots exist iff gcd(f, X^p − X) ≠ 1; deg gcd counts the distinct roots465    ok, det = True, []466    for p, dmax in [(13, 3), (17, 3)]:467        xp_minus_x = [0] * (p + 1); xp_minus_x[p] = 1; xp_minus_x[1] = (xp_minus_x[1] - 1) % p468        n = 0469        for d in range(1, dmax + 1):470            for coeffs in itertools.product(range(p), repeat=d):471                f = list(coeffs) + [1]                                  # monic of degree d472                roots = {x for x in range(p) if p_eval(f, x, p) == 0}473                g = p_gcd(f, xp_minus_x, p)474                ok &= ((len(g) - 1) == len(roots)) and ((len(roots) > 0) == (len(g) > 1)); n += 1475        det.append(f"p={p}: {n} monic polynomials of degree ≤ {dmax}")476    # 1:E3 (p01026)477    check("C1", "f has a root in F_p iff gcd(f, X^p − X) ≠ 1; deg gcd(f, X^p − X) = number of distinct roots (exhaustive, degree ≤ 3)",478          ok, "; ".join(det))479480    # C2 — limit-like approximation across the tower of shells481    getcontext().prec = 60482    reals = {"π": Fraction(Decimal("3.14159265358979323846264338327950288419716939937510582097494")),483             "e": Fraction(Decimal("2.71828182845904523536028747135266249775724709369995957496697")),484             "√2": Fraction(Decimal("1.41421356237309504880168872420969807856967187537694807317668")),485             "33/10": Fraction(33, 10), "1/3": Fraction(1, 3)}486    def next_shell(m):                                                  # the first prime ≡ 1 (mod 4) beyond m487        q = m + 1488        while not (q % 4 == 1 and is_prime(q)): q += 1489        return q490    ok, det = True, []491    for name, r in reals.items():492        for k in range(2, 9):493            eps = Fraction(1, 10 ** k); n = math.ceil(math.log2(10 ** k)); x = round(r * 2 ** n)494            p = next_shell(2 * abs(x) + 1); kap = (p - 1) // 4495            ok &= (abs(r - Fraction(x, 2 ** n)) < eps) and (abs(x) <= 2 * kap) and p % 4 == 1 and is_prime(p)496        det.append(f"{name}: ε=10⁻⁸ → n={n}, p={p}")497    # 1:D7 (p01027)498    check("C2", "for every r and ε some shell p = 4κ+1 carries x/2ⁿ, |x| ≤ 2κ, within ε of r (five reals, ε = 10⁻²..10⁻⁸)",499          ok, "; ".join(det), kind="CHART")500501    # C3 — the circle net and the group-law defect on every shell502    ok, det = True, []503    rnd = random.Random(1)504    for p in SHELLS:505        N = p - 1506        k = lambda th: math.floor(N * th / (2 * math.pi) + 0.5)507        worst_angle = worst_chord = 0.0; worst_defect = 0508        thetas = [2 * math.pi * j / (4 * N) for j in range(4 * N)] + [rnd.uniform(0, 2 * math.pi) for _ in range(2000)]509        for th in thetas:510            kk = k(th); worst_angle = max(worst_angle, abs(th - 2 * math.pi * kk / N))511            worst_chord = max(worst_chord, abs(complex(math.cos(2 * math.pi * kk / N), math.sin(2 * math.pi * kk / N)) - complex(math.cos(th), math.sin(th))))512        for _ in range(20000):513            a, b = rnd.uniform(0, 2 * math.pi), rnd.uniform(0, 2 * math.pi)514            worst_defect = max(worst_defect, abs(k(a) + k(b) - k(a + b)))515        ok &= (worst_angle <= math.pi / N + 1e-12) and (worst_chord <= math.pi / N + 1e-12) and (worst_defect <= 1)516        det.append(f"p={p}: angle {worst_angle:.4f} ≤ π/N={math.pi / N:.4f}, defect ≤ {worst_defect}")517    # 1:E4 (p01028)518    check("C3", "k(θ) = ⌊Nθ/2π + ½⌋, N = p−1: angle and chord error ≤ π/N, group-law defect ≤ 1, on every shell",519          ok, "; ".join(det), kind="CHART")520521    # C4 — the SO(3) obstruction: covering radii of the finite rotation groups522    def qmul(a, b):523        w1, x1, y1, z1 = a; w2, x2, y2, z2 = b524        return (w1*w2 - x1*x2 - y1*y2 - z1*z2, w1*x2 + x1*w2 + y1*z2 - z1*y2, w1*y2 - x1*z2 + y1*w2 + z1*x2, w1*z2 + x1*y2 - y1*x2 + z1*w2)525    def rot(axis, ang):526        n = math.sqrt(sum(c * c for c in axis)); s = math.sin(ang / 2) / n527        return (math.cos(ang / 2), s * axis[0], s * axis[1], s * axis[2])528    def closure(gens):529        G, frontier = [(1.0, 0.0, 0.0, 0.0)], [(1.0, 0.0, 0.0, 0.0)]530        while frontier:531            new = []532            for a in frontier:533                for g in gens:534                    b = qmul(a, g)535                    if not any(sum((u - v) ** 2 for u, v in zip(b, c)) < 1e-18 for c in G): G.append(b); new.append(b)536            frontier = new537        return G538    def dist(q, G):                                                     # rotation-angle distance from q to the group539        return 2 * math.acos(min(1.0, max(abs(sum(u * v for u, v in zip(q, g))) for g in G)))540    def covering_radius(G, samples=20000, refine=60):541        best = 0.0542        pts = []543        for _ in range(samples):544            v = [rnd.gauss(0, 1) for _ in range(4)]; n = math.sqrt(sum(c * c for c in v)); pts.append(tuple(c / n for c in v))545        pts.sort(key=lambda q: -dist(q, G))546        for q in pts[:12]:547            cur, step = dist(q, G), 0.05548            for _ in range(refine):549                cands = []550                for _ in range(48):551                    v = [c + step * rnd.gauss(0, 1) for c in q]; n = math.sqrt(sum(c * c for c in v)); cands.append(tuple(c / n for c in v))552                c = max(cands, key=lambda q: dist(q, G)); dc = dist(c, G)553                if dc > cur: q, cur = c, dc554                else: step *= 0.7555            best = max(best, cur)556        return best557    phi = (1 + 5 ** 0.5) / 2558    groups = {"tetrahedral": closure([rot((0, 0, 1), math.pi), rot((1, 1, 1), 2 * math.pi / 3)]),559              "octahedral": closure([rot((0, 0, 1), math.pi / 2), rot((1, 1, 1), 2 * math.pi / 3)]),560              "icosahedral": closure([rot((0, 0, 1), math.pi), rot((1, 1, 1), 2 * math.pi / 3), rot((0, 1, phi), 2 * math.pi / 5)])}561    exact = {"tetrahedral": math.pi / 2, "octahedral": math.acos((2 * 2 ** 0.5 - 1) / 4), "icosahedral": math.acos((3 * 5 ** 0.5 - 1) / 8)}562    ok, det = True, []563    ok &= [len(G) for G in groups.values()] == [24, 48, 120]           # the binary groups: 2 × 12, 24, 60564    for name, G in groups.items():565        r = covering_radius(G); ok &= abs(r - exact[name]) < 2e-3566        det.append(f"{name}: {math.degrees(r):.2f}° (exact {math.degrees(exact[name]):.2f}°)")567    # cyclic and dihedral groups lie in an O(2): rotations about an axis and half-turns about axes normal to it; the568    # rotation q = (w,x,y,z) with w²+z² = x²+y² = ½ is at distance 2·arccos(1/√2) = π/2 from every element of O(2)569    q = (0.5 ** 0.5 * math.cos(0.3), 0.5 ** 0.5 * math.cos(1.1), 0.5 ** 0.5 * math.sin(1.1), 0.5 ** 0.5 * math.sin(0.3))570    d_o2 = 2 * math.acos(max(math.hypot(q[0], q[3]), math.hypot(q[1], q[2])))571    ok &= abs(d_o2 - math.pi / 2) < 1e-12572    eps0 = min(exact.values()); ok &= abs(eps0 - exact["icosahedral"]) < 1e-15 and eps0 < math.pi / 2573    det.append(f"cyclic/dihedral (in O(2)): ≥ {math.degrees(d_o2):.2f}°; ε₀ = {math.degrees(eps0):.4f}° = {eps0:.6f} rad")574    # 1:E5 (p01029)575    check("C4", "covering radii in SO(3): tetrahedral π/2, octahedral arccos((2√2−1)/4), icosahedral arccos((3√5−1)/8) = ε₀ ≈ 44.48°, cyclic/dihedral ≥ π/2; no finite subgroup is an ε-net for ε < ε₀",576          ok, "; ".join(det), kind="CHART")577578    # C5 — the window resolves SO(3): the normalised framed quaternions are a 2·arcsin(1/H)-net, composing exactly on the shell579    ok, det = True, []580    for H in [1, 2, 3]:581        pts = [q for q in itertools.product(range(-H, H + 1), repeat=4) if any(q)]582        seen, W = set(), []583        for q in pts:                                                   # one representative per rotation (q and −q)584            n = math.sqrt(sum(c * c for c in q)); u = tuple(c / n for c in q)585            key = tuple(round(c, 9) for c in (u if next(c for c in u if abs(c) > 1e-12) > 0 else tuple(-c for c in u)))586            if key not in seen: seen.add(key); W.append(u)587        bound = 2 * math.asin(1 / H)588        samples = []589        for _ in range(3000):590            v = [rnd.gauss(0, 1) for _ in range(4)]; n = math.sqrt(sum(c * c for c in v)); samples.append(tuple(c / n for c in v))591        ds = [dist(q, W) for q in samples]592        ok &= max(ds) <= bound + 1e-12                                   # the elementary bound, on every sample593        # the rounding witness itself, on every sample: q = round(H s) is within arcsin(1/H)594        for s_ in samples:595            q = tuple(round(H * c) for c in s_)596            if any(q):597                n = math.sqrt(sum(c * c for c in q)); cosang = abs(sum(a * b for a, b in zip(s_, q))) / n598                ok &= 2 * math.acos(min(1.0, cosang)) <= bound + 1e-12599        measured = covering_radius(W, samples=4000, refine=40)600        ok &= measured <= bound + 1e-9601        det.append(f"H={H}: {len(W)} rotations, measured {math.degrees(measured):.1f}° ≤ bound {math.degrees(bound):.1f}°")602    ok &= 2 * math.asin(1 / 3) < math.acos((3 * 5 ** 0.5 - 1) / 8)     # H = 3 beats every finite subgroup603    # exact composition on the shell p = 73 > 8·3² = 72: the product of two window quaternions reads back from F_73604    p, H = 73, 3605    def qmul_int(a, b):606        w1, x1, y1, z1 = a; w2, x2, y2, z2 = b607        return (w1*w2 - x1*x2 - y1*y2 - z1*z2, w1*x2 + x1*w2 + y1*z2 - z1*y2, w1*y2 - x1*z2 + y1*w2 + z1*x2, w1*z2 + x1*y2 - y1*x2 + z1*w2)608    read = lambda r: r - p if r > (p - 1) // 2 else r                   # the window reading of a residue609    n_pairs = 0610    for _ in range(20000):611        a = tuple(rnd.randint(-H, H) for _ in range(4)); b = tuple(rnd.randint(-H, H) for _ in range(4))612        exact = qmul_int(a, b); modp = qmul_int(tuple(c % p for c in a), tuple(c % p for c in b))613        ok &= all(read(c % p) == e for c, e in zip(modp, exact)) and all(abs(e) <= 4 * H * H for e in exact); n_pairs += 1614    det.append(f"2·arcsin(1/3) = {math.degrees(2 * math.asin(1/3)):.1f}° < ε₀; composition on F_{p}, H={H}: {n_pairs} pairs read back exactly")615    # 1:E6 (p01030)616    check("C5", "the normalised window quaternions W_H⁴ are a 2·arcsin(1/H)-net of SO(3) (H = 1, 2, 3), 38.9° < ε₀ at H = 3; their products read back exactly from the shell when 8H² < p",617          ok, "; ".join(det), kind="CHART")618619    # ---------------------------------------------------------------------------------------------------- the Lie-algebra layer620    # Quaternion rotations over F_p as 3×3 matrices: ρ_q(x) = q x q̄ / N(q) on pure x, S_v x = v × x, N(v) = v·v.621    def qm(x, y):622        a1, b1, c1, d1 = x; a2, b2, c2, d2 = y623        return (a1*a2 - b1*b2 - c1*c2 - d1*d2, a1*b2 + b1*a2 + c1*d2 - d1*c2, a1*c2 - b1*d2 + c1*a2 + d1*b2, a1*d2 + b1*c2 - c1*b2 + d1*a2)624    def mm(A, B, p): return tuple(tuple(sum(A[i][k] * B[k][j] for k in range(3)) % p for j in range(3)) for i in range(3))625    def minv(A, p):626        (a, b, c), (d, e, f), (g, h, i) = A627        det = (a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g)) % p628        if det == 0: return None629        di = pow(det, -1, p)630        adj = ((e*i - f*h, c*h - b*i, b*f - c*e), (f*g - d*i, a*i - c*g, c*d - a*f), (d*h - e*g, b*g - a*h, a*e - b*d))631        return tuple(tuple((x * di) % p for x in r) for r in adj)632    I3 = ((1, 0, 0), (0, 1, 0), (0, 0, 1))633    def rho(q, p):634        n = sum(c * c for c in q) % p; ni = pow(n, -1, p); cols = []635        for e in ((0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)):636            r = qm(qm(q, e), (q[0], -q[1], -q[2], -q[3])); cols.append([(c * ni) % p for c in r[1:]])637        return tuple(tuple(cols[j][i] for j in range(3)) for i in range(3))638    def Smat(v, p): x, y, z = v; return tuple(tuple(c % p for c in r) for r in ((0, -z, y), (z, 0, -x), (-y, x, 0)))639    def cayley(a, v, p):640        S = Smat(v, p); M = minv(tuple(tuple((a * I3[i][j] - S[i][j]) % p for j in range(3)) for i in range(3)), p)641        return None if M is None else mm(M, tuple(tuple((a * I3[i][j] + S[i][j]) % p for j in range(3)) for i in range(3)), p)642    def halfturn(v, p):643        S = Smat(v, p); S2 = mm(S, S, p); ni = pow(sum(c * c for c in v) % p, -1, p)644        return tuple(tuple((I3[i][j] + 2 * ni * S2[i][j]) % p for j in range(3)) for i in range(3))645    def legendre(n, p):646        n %= p647        return 0 if n == 0 else (1 if pow(n, (p - 1) // 2, p) == 1 else -1)648    def nonsquare(p): return next(x for x in range(2, p) if legendre(x, p) == -1)649    def eta_cayley_ok(a, v, p):650        # 8:C3 in F_p[η], η² = ν: α = η/(νa), H = ηS_v, so αH = (η²/(νa)) S_v; its η-part vanishes and its F_p-part is S_v/a651        nu = nonsquare(p); inv = lambda x: pow(x % p, -1, p)652        alpha = (0, inv(nu * a))                                  # 0 + (1/(νa)) η653        S = Smat(v, p)654        for i in range(3):655            for j in range(3):656                h = (0, S[i][j])                                  # η S_ij657                prod = ((alpha[0] * h[0] + nu * alpha[1] * h[1]) % p, (alpha[0] * h[1] + alpha[1] * h[0]) % p)658                if prod != ((S[i][j] * inv(a)) % p, 0): return False659        return True660    def primitive(v):661        g = 0662        for c in v: g = math.gcd(g, abs(c))663        v = tuple(c // g for c in v)664        for c in v:                                               # sign: first nonzero coordinate positive665            if c: return v if c > 0 else tuple(-x for x in v)666    def group_order(M, p):667        k, X = 1, M668        while X != I3: X = mm(X, M, p); k += 1669        return k670671    ok, det = True, []672    n_rot = n_axes = 0673    for p, H in ((13, 1), (17, 1), (37, 2), (73, 3), (101, 3)):674        assert p > 8 * H * H675        window = [q for q in itertools.product(range(-H, H + 1), repeat=4) if any(q[1:])]676        axes = {}677        for q in window:678            a, v = q[0], q[1:]679            R = rho(q, p); n_rot += 1680            ok &= R == (cayley(a % p, v, p) if a % p else halfturn(v, p))681            if a % p: ok &= eta_cayley_ok(a % p, v, p)682            axes.setdefault(primitive(v), set()).add(R)683        for v0, rots in axes.items():684            n_axes += 1685            N = sum(c * c for c in v0) % p686            G = {I3, halfturn(v0, p)} | {C for C in (cayley(a, v0, p) for a in range(1, p)) if C is not None}687            order = p - legendre(N, p)688            ok &= len(G) == order                                                         # the order p − (N(v)|p)689            gen = next((M for M in G if group_order(M, p) == order), None)                # cyclic: an element of full order690            powers, X = {I3}, gen691            while gen is not None and X != I3: powers.add(X); X = mm(X, gen, p)692            ok &= gen is not None and powers == G                                         # whose powers are exactly the set693            for R in rots:                                                                # every iterate of a window rotation stays in it694                X = R695                for _ in range(order): ok &= X in G; X = mm(X, R, p)696            ok &= len(rots) <= 2 * H * H + 2 * H                                         # at most 2H² + 2H window rotations697    det.append(f"{n_rot} window rotations on the shells 13, 17, 37, 73, 101 (8H² < p): ρ_q = (aI − S_v)⁻¹(aI + S_v) = the Cayley step of ηS_v at α = η/(νa), and I + 2N(v)⁻¹S_v² at a = 0")698    det.append(f"{n_axes} axes: the steps, the identity and the half-turn form a cyclic group of order p − (N(v)|p) holding every iterate; at most 2H² + 2H window rotations per axis")699    # 1:E7 (p01032)700    check("C6", "ρ_q = (aI − S_v)⁻¹(aI + S_v), the Cayley step of ηS_v at α = η/(νa), and the half-turn I + 2N(v)⁻¹S_v² at a = 0; with it the steps form a cyclic group of order p − (N(v)|p) holding every iterate; ≤ 2H² + 2H window rotations per axis",701          ok, "; ".join(det))702703    ok, det = True, []704    rnd = random.Random(20261004)705    def cross(u, w): return (u[1]*w[2] - u[2]*w[1], u[2]*w[0] - u[0]*w[2], u[0]*w[1] - u[1]*w[0])706    def dot(u, w): return sum(x * y for x, y in zip(u, w))707    n_comp = n_half = n_br = 0708    for p in (13, 29, 73, 233):709        inv = lambda x: pow(x % p, -1, p)710        tries = 0711        while n_comp < 4000 * (1 + (13, 29, 73, 233).index(p)) and tries < 100000:712            tries += 1713            u = tuple(rnd.randrange(p) for _ in range(3)); w = tuple(rnd.randrange(p) for _ in range(3))714            uw = dot(u, w) % p715            if (1 + dot(u, u)) % p == 0 or (1 + dot(w, w)) % p == 0: continue716            s = (1 - uw) % p717            m = tuple((u[i] + w[i] + cross(u, w)[i]) % p for i in range(3))718            lhs = mm(rho((1,) + u, p), rho((1,) + w, p), p)719            if s:720                c = tuple((x * inv(s)) % p for x in m)721                if (1 + dot(c, c)) % p == 0: continue722                ok &= lhs == rho((1,) + c, p); n_comp += 1723            elif any(m) and dot(m, m) % p:724                ok &= lhs == halfturn(m, p); n_half += 1725            br = tuple((x - y) % p for x, y in zip(qm((0,) + u, (0,) + w), qm((0,) + w, (0,) + u)))726            ok &= br == (0,) + tuple((2 * x) % p for x in cross(u, w)); n_br += 1727        # the half-turn case on purpose: w chosen with u·w = 1728        for _ in range(300):729            u = tuple(rnd.randrange(p) for _ in range(3))730            if not u[0] % p: continue731            w1, w2 = rnd.randrange(p), rnd.randrange(p)732            w = (((1 - u[1] * w1 - u[2] * w2) * inv(u[0])) % p, w1, w2)733            m = tuple((u[i] + w[i] + cross(u, w)[i]) % p for i in range(3))734            if (1 + dot(u, u)) % p == 0 or (1 + dot(w, w)) % p == 0 or not any(m) or dot(m, m) % p == 0: continue735            ok &= mm(rho((1,) + u, p), rho((1,) + w, p), p) == halfturn(m, p); n_half += 1736    det.append(f"composition ρ_(1+u)ρ_(1+w) = ρ_(1+u∘w): {n_comp} pairs; the half-turn at u·w = 1: {n_half}; [u,w] = 2u×w: {n_br} (F_13, F_29, F_73, F_233)")737    n_pow = n_rb = 0738    for H in (1, 2, 3):739        for q in itertools.product(range(-H, H + 1), repeat=4):740            if not any(q): continue741            r = (1, 0, 0, 0)742            for m_ in range(1, 7):743                r = qm(r, q); n_pow += 1744                ok &= max(map(abs, r)) <= (2 * H) ** m_745                for p in (13, 101, 1009):746                    reads = all(((c % p) - p if (c % p) > (p - 1) // 2 else (c % p)) == c for c in r)747                    ok &= reads == (2 * max(map(abs, r)) < p); n_rb += 1748    det.append(f"heights: {n_pow} powers q^m (H ≤ 3, m ≤ 6) with entries ≤ (2H)^m; read-back iff 2‖q^m‖∞ < p: {n_rb} cases on F_13, F_101, F_1009")749    # 1:E8 (p01033)750    check("C7", "ρ_(1+u)ρ_(1+w) = ρ_(1+u∘w), u∘w = (u+w+u×w)/(1−u·w); at u·w = 1 the half-turn about u+w+u×w; [u,w] = 2u×w; entries of q^m ≤ (2H)^m, read back iff 2‖q^m‖∞ < p",751          ok, "; ".join(det))752753if __name__ == "__main__":754    import time755    want = [a.upper() for a in sys.argv[1:]] or sorted(BLOCK)756    bad = [b for b in want if b not in BLOCK]757    if bad: sys.exit(f"no block {', '.join(bad)}: the blocks are {', '.join(sorted(BLOCK))}")758    t0 = time.time()759    for b in want:760        t = time.time(); _run_block(b); print(f"    [block {b}: {time.time() - t:.1f} s]")761    ok = summary(write=(want == sorted(BLOCK)))762    kinds = {}763    for r in RESULTS: kinds[r["kind"]] = kinds.get(r["kind"], 0) + 1764    print("by kind: " + ", ".join(f"{k} {v}" for k, v in sorted(kinds.items())) + f"; {len(RESULTS)} checks in {time.time() - t0:.1f} s" + ("; results.json written" if want == sorted(BLOCK) else ""))765    sys.exit(0 if ok else 1)766