src/1-algebra — algebra.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/src/1-algebra/#<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; 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": "",          # 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",                       # 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; 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 59286def block_B():287    """Block B — the framed numbers, the charts and the horizon (EXACT; B3 CHART in exact rationals): B1–B6."""288    print("block B — the framed numbers, the charts and the horizon")289290    # B1 — the window law, the three thresholds as iff over every H291    ok, det = True, []292    for p in [13, 17, 29]:293        for H in range(1, p):294            W = range(-H, H + 1)295            inj = len({z % p for z in W}) == len(W)296            sums = len({(x + y) % p for x in W for y in W}) == len({x + y for x in W for y in W})297            prods = len({(x * y) % p for x in W for y in W}) == len({x * y for x in W for y in W})298            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)299        det.append(f"p={p}: H = 1..{p-1} swept")300    # 1:D2 (p01013)301    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))302303    # B2 — scale-periodicity: in the field yes, in Q no304    ok, det = True, []305    for p, g in [(13, 11), (13, 2), (17, 3), (29, 2)]:306        assert g in primitive_roots(p)307        grid = lambda n: tuple((x * pow(g, -n, p)) % p for x in range(p))308        ok &= all(grid(n) == grid(n + p - 1) for n in range(2 * (p - 1)))309        ok &= all(grid(n) != grid(n + d) for n in range(p - 1) for d in range(1, p - 1))   # exact period p−1310        qgrid = lambda n: tuple(Fraction(x, g ** n) for x in range(p))311        ok &= all(qgrid(n) != qgrid(n + p - 1) for n in range(3))312        det.append(f"({p},{g}): period {p-1} in F_p, none in Q")313    # 1:D4 (p01015)314    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",315          ok, "; ".join(det))316317    # B3 — the chart of the grid: the trade-off holds; its range at (13, 2)318    ok, det = True, []319    p, g, H = 13, 2, 12320    for n in range(0, 8):321        step = Fraction(1, g ** n)322        for num in range(-40, 41):                                    # r = num/8 · H g^{−n}, a sample within the range323            r = Fraction(num, 40) * H * step324            x = round(r / step)                                        # x = round(r gⁿ), |x| ≤ H325            ok &= (abs(x) <= H) and (abs(r - x * step) <= step / 2)326    r, tol = Fraction(33, 10), Fraction(1, 16)327    near = [(x, n) for n in range(3, 13) for x in range(p) if abs(Fraction(x, g ** n) - r) < tol]328    ok &= (near == []) and all(Fraction(x, g ** n) <= Fraction(3, 2) for n in range(3, 13) for x in range(p))329    best = min((abs(Fraction(x, g ** n) - r), x, n) for n in range(13) for x in range(p))330    ok &= (best[1:] == (7, 1) and best[0] == Fraction(1, 5))331    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")332    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",333          ok, "; ".join(det), kind="CHART")334335    # B4 — the complex chart is not an extension: zero divisors iff −1 is a square336    ok, det = True, []337    for p in [5, 13, 17, 29, 37]:338        u = sqrt_neg_one(p)[0]339        # in F_p[X]/(X²+1): elements a + bX, product (a+bX)(c+dX) = (ac − bd) + (ad + bc)X340        mul = lambda A, B: ((A[0] * B[0] - A[1] * B[1]) % p, (A[0] * B[1] + A[1] * B[0]) % p)341        ok &= (mul((u, 1), ((-u) % p, 1)) == (0, 0)) and (u, 1) != (0, 0) and ((-u) % p, 1) != (0, 0)342        det.append(f"p={p}: (u+X)(u−X)=0 with u={u}")343    for p in CONTROLS:344        ok &= (sqrt_neg_one(p) == [])                                    # X²+1 irreducible: the quotient is the field F_{p²}345        mul = lambda A, B: ((A[0] * B[0] - A[1] * B[1]) % p, (A[0] * B[1] + A[1] * B[0]) % p)346        elems = [(a, b) for a in range(p) for b in range(p) if (a, b) != (0, 0)]347        ok &= all(mul(A, B) != (0, 0) for A in elems for B in elems)    # no zero divisors348    # 1:E2 (p01019)349    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)",350          ok, "; ".join(det) + f"; fields at p ∈ {CONTROLS} (no zero divisors, exhaustive)")351352    # B5 — no element of additive order two; the half-turn353    ok, det = True, []354    for p in [q for q in range(3, 200) if is_prime(q)]:355        ok &= ([s for s in range(p) if (2 * s) % p == 0] == [0])356        if p % 4 == 1:357            k = kappa(p)358            ok &= ((2 * (2 * k + 1)) % p == 1) and (2 * k + 1 == (p + 1) // 2)359    # 1:Z1 (p01020)360    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")361362    # B6 — the Euclidean step count against the bound ⌊log₂ p⌋+1 (recorded under V1)363    def steps(a, b):364        k = 0365        while b:366            a, b = b, a % b367            if b: k += 1368        return k369    fib = [1, 1]370    while fib[-1] < 2000: fib.append(fib[-1] + fib[-2])371    first = None372    for p in [q for q in range(5, 2000) if is_prime(q)]:373        worst = max(steps(fib[i + 1], fib[i]) for i in range(1, len(fib) - 1) if fib[i + 1] < p)374        if worst > p.bit_length():375            first = p; break376    ok = (steps(55, 34) == 7 and (59).bit_length() == 6 and steps(987, 610) == 13 and (1009).bit_length() == 10 and first == 59)377    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 = 4378    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",379          ok, "a ≥ b convention, consecutive Fibonacci pairs below p")380381# ------------------------------------------------------------------------------------------------------------382# block C: the conjecture of the conclusion, clause by clause (EXACT; C2–C4 CHART)383# The conclusion conjectures that the finite substrate supports polynomial equation solving, limit-like384# approximation and ε-approximation of continuous symmetries.  Predicates 1:E3, D7 and E4–E8 decide it:385#386#   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 of387#               distinct roots — every monic polynomial of degree ≤ 3 over F_13 (2197) and F_17 (4913), by388#               brute-force roots against the Euclidean gcd            [exact]389#   C2  (1:D7)  limit-like approximation across shells: for r ∈ {π, e, √2, 33/10, 1/3} and ε = 10^−k390#               (k = 2..8) the construction n = ⌈log₂(1/ε)⌉, x = round(r·2ⁿ), p = the first prime ≡ 1 (mod 4)391#               beyond 2|x|+1 gives |r − x/2ⁿ| < ε with |x| ≤ 2κ      [chart: r read as a rational to 60 digits]392#   C3  (1:E4)  the circle net: for N = p − 1 on every shell the rounding k(θ) = ⌊Nθ/2π + ½⌋ has angle error393#               ≤ π/N, chord error |e^{2πik/N} − e^{iθ}| ≤ π/N, and group-law defect |k(θ₁)+k(θ₂)−k(θ₁+θ₂)| ≤ 1,394#               on a grid of 4N angles plus 20000 random pairs           [chart]395#   C4  (1:E5)  the SO(3) obstruction: the covering radii of the finite rotation groups in the rotation-angle396#               metric — tetrahedral π/2, octahedral arccos((2√2−1)/4) ≈ 62.80°, icosahedral397#               ε₀ = arccos((3√5−1)/8) ≈ 44.48° (the deep hole of the 600-cell, cos = φ²/2√2), cyclic and398#               dihedral ≥ π/2 (they lie in an O(2)) — by exhaustive group closure and a sampled maximin;399#               no finite subgroup of SO(3) is an ε-net for ε < ε₀      [chart; the list of groups is A3, Klein]400#   C5  (1:E6)  the window resolves SO(3): the framed quaternions W_H⁴ = {q ∈ Z⁴ : |q_i| ≤ H}, normalised, are an401#               ε-net of SO(3) with ε ≤ 2·arcsin(1/H) (round H·s to the lattice: |Hs − q| ≤ 1, so the angle is402#               ≤ arcsin(1/H)); measured ε(1) ≈ 60.8°, ε(2) ≈ 41.0°, ε(3) ≈ 30.0°; 2·arcsin(1/3) = 38.9° < ε₀, so on403#               the shell p = 73 the window H = 3 out-resolves every finite subgroup; products of window quaternions404#               have entries in W_{4H²} and read back exactly from the shell when 8H² < p (20000 random pairs on405#               F_73, H = 3)                                              [chart; the read-back is exact]406#   C6  (1:E7)  the Lie-algebra layer, one axis: every window rotation ρ_q, q = a + v, equals (aI − S_v)⁻¹(aI + S_v) — the407#               Cayley step of ηS_v at α = η/(νa), computed in F_p[η] — and I + 2N(v)⁻¹S_v² at a = 0; on every axis the408#               steps, the identity and the half-turn are the powers of one element of order p − (N(v)|p), every409#               iterate stays among them, and at most 2H² + 2H are window rotations (p = 13, 17, 37, 73, 101) [exact]410#   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 about411#               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 and412#               q^m reads back from the shell iff 2‖q^m‖∞ < p (H ≤ 3, m ≤ 6; p = 13, 101, 1009)                   [exact]413414# ---------------------------------------------------------------- polynomial arithmetic over F_p (lists, low degree first)415def p_trim(a):416    while a and a[-1] == 0: a.pop()417    return a418def p_mod(a, b, p):419    a = a[:]; db = len(b) - 1; inv = pow(b[-1], -1, p)420    while len(a) - 1 >= db and a:421        c = (a[-1] * inv) % p; s = len(a) - 1 - db422        for i, bi in enumerate(b): a[s + i] = (a[s + i] - c * bi) % p423        p_trim(a)424    return a425def p_gcd(a, b, p):426    a, b = p_trim(a[:]), p_trim(b[:])427    while b: a, b = b, p_mod(a, b, p)428    if a: inv = pow(a[-1], -1, p); a = [(c * inv) % p for c in a]429    return a430def p_eval(a, x, p):431    r = 0432    for c in reversed(a): r = (r * x + c) % p433    return r434435def block_C():436    """Block C — the conjecture of the conclusion, clause by clause (EXACT; C2–C5 CHART): C1–C5."""437    print("block C — the conjecture of the conclusion, clause by clause")438439    # C1 — solving: roots exist iff gcd(f, X^p − X) ≠ 1; deg gcd counts the distinct roots440    ok, det = True, []441    for p, dmax in [(13, 3), (17, 3)]:442        xp_minus_x = [0] * (p + 1); xp_minus_x[p] = 1; xp_minus_x[1] = (xp_minus_x[1] - 1) % p443        n = 0444        for d in range(1, dmax + 1):445            for coeffs in itertools.product(range(p), repeat=d):446                f = list(coeffs) + [1]                                  # monic of degree d447                roots = {x for x in range(p) if p_eval(f, x, p) == 0}448                g = p_gcd(f, xp_minus_x, p)449                ok &= ((len(g) - 1) == len(roots)) and ((len(roots) > 0) == (len(g) > 1)); n += 1450        det.append(f"p={p}: {n} monic polynomials of degree ≤ {dmax}")451    # 1:E3 (p01026)452    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)",453          ok, "; ".join(det))454455    # C2 — limit-like approximation across the tower of shells456    getcontext().prec = 60457    reals = {"π": Fraction(Decimal("3.14159265358979323846264338327950288419716939937510582097494")),458             "e": Fraction(Decimal("2.71828182845904523536028747135266249775724709369995957496697")),459             "√2": Fraction(Decimal("1.41421356237309504880168872420969807856967187537694807317668")),460             "33/10": Fraction(33, 10), "1/3": Fraction(1, 3)}461    def next_shell(m):                                                  # the first prime ≡ 1 (mod 4) beyond m462        q = m + 1463        while not (q % 4 == 1 and is_prime(q)): q += 1464        return q465    ok, det = True, []466    for name, r in reals.items():467        for k in range(2, 9):468            eps = Fraction(1, 10 ** k); n = math.ceil(math.log2(10 ** k)); x = round(r * 2 ** n)469            p = next_shell(2 * abs(x) + 1); kap = (p - 1) // 4470            ok &= (abs(r - Fraction(x, 2 ** n)) < eps) and (abs(x) <= 2 * kap) and p % 4 == 1 and is_prime(p)471        det.append(f"{name}: ε=10⁻⁸ → n={n}, p={p}")472    # 1:D7 (p01027)473    check("C2", "for every r and ε some shell p = 4κ+1 carries x/2ⁿ, |x| ≤ 2κ, within ε of r (five reals, ε = 10⁻²..10⁻⁸)",474          ok, "; ".join(det), kind="CHART")475476    # C3 — the circle net and the group-law defect on every shell477    ok, det = True, []478    rnd = random.Random(1)479    for p in SHELLS:480        N = p - 1481        k = lambda th: math.floor(N * th / (2 * math.pi) + 0.5)482        worst_angle = worst_chord = 0.0; worst_defect = 0483        thetas = [2 * math.pi * j / (4 * N) for j in range(4 * N)] + [rnd.uniform(0, 2 * math.pi) for _ in range(2000)]484        for th in thetas:485            kk = k(th); worst_angle = max(worst_angle, abs(th - 2 * math.pi * kk / N))486            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))))487        for _ in range(20000):488            a, b = rnd.uniform(0, 2 * math.pi), rnd.uniform(0, 2 * math.pi)489            worst_defect = max(worst_defect, abs(k(a) + k(b) - k(a + b)))490        ok &= (worst_angle <= math.pi / N + 1e-12) and (worst_chord <= math.pi / N + 1e-12) and (worst_defect <= 1)491        det.append(f"p={p}: angle {worst_angle:.4f} ≤ π/N={math.pi / N:.4f}, defect ≤ {worst_defect}")492    # 1:E4 (p01028)493    check("C3", "k(θ) = ⌊Nθ/2π + ½⌋, N = p−1: angle and chord error ≤ π/N, group-law defect ≤ 1, on every shell",494          ok, "; ".join(det), kind="CHART")495496    # C4 — the SO(3) obstruction: covering radii of the finite rotation groups497    def qmul(a, b):498        w1, x1, y1, z1 = a; w2, x2, y2, z2 = b499        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)500    def rot(axis, ang):501        n = math.sqrt(sum(c * c for c in axis)); s = math.sin(ang / 2) / n502        return (math.cos(ang / 2), s * axis[0], s * axis[1], s * axis[2])503    def closure(gens):504        G, frontier = [(1.0, 0.0, 0.0, 0.0)], [(1.0, 0.0, 0.0, 0.0)]505        while frontier:506            new = []507            for a in frontier:508                for g in gens:509                    b = qmul(a, g)510                    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)511            frontier = new512        return G513    def dist(q, G):                                                     # rotation-angle distance from q to the group514        return 2 * math.acos(min(1.0, max(abs(sum(u * v for u, v in zip(q, g))) for g in G)))515    def covering_radius(G, samples=20000, refine=60):516        best = 0.0517        pts = []518        for _ in range(samples):519            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))520        pts.sort(key=lambda q: -dist(q, G))521        for q in pts[:12]:522            cur, step = dist(q, G), 0.05523            for _ in range(refine):524                cands = []525                for _ in range(48):526                    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))527                c = max(cands, key=lambda q: dist(q, G)); dc = dist(c, G)528                if dc > cur: q, cur = c, dc529                else: step *= 0.7530            best = max(best, cur)531        return best532    phi = (1 + 5 ** 0.5) / 2533    groups = {"tetrahedral": closure([rot((0, 0, 1), math.pi), rot((1, 1, 1), 2 * math.pi / 3)]),534              "octahedral": closure([rot((0, 0, 1), math.pi / 2), rot((1, 1, 1), 2 * math.pi / 3)]),535              "icosahedral": closure([rot((0, 0, 1), math.pi), rot((1, 1, 1), 2 * math.pi / 3), rot((0, 1, phi), 2 * math.pi / 5)])}536    exact = {"tetrahedral": math.pi / 2, "octahedral": math.acos((2 * 2 ** 0.5 - 1) / 4), "icosahedral": math.acos((3 * 5 ** 0.5 - 1) / 8)}537    ok, det = True, []538    ok &= [len(G) for G in groups.values()] == [24, 48, 120]           # the binary groups: 2 × 12, 24, 60539    for name, G in groups.items():540        r = covering_radius(G); ok &= abs(r - exact[name]) < 2e-3541        det.append(f"{name}: {math.degrees(r):.2f}° (exact {math.degrees(exact[name]):.2f}°)")542    # cyclic and dihedral groups lie in an O(2): rotations about an axis and half-turns about axes normal to it; the543    # 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)544    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))545    d_o2 = 2 * math.acos(max(math.hypot(q[0], q[3]), math.hypot(q[1], q[2])))546    ok &= abs(d_o2 - math.pi / 2) < 1e-12547    eps0 = min(exact.values()); ok &= abs(eps0 - exact["icosahedral"]) < 1e-15 and eps0 < math.pi / 2548    det.append(f"cyclic/dihedral (in O(2)): ≥ {math.degrees(d_o2):.2f}°; ε₀ = {math.degrees(eps0):.4f}° = {eps0:.6f} rad")549    # 1:E5 (p01029)550    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 ε < ε₀",551          ok, "; ".join(det), kind="CHART")552553    # C5 — the window resolves SO(3): the normalised framed quaternions are a 2·arcsin(1/H)-net, composing exactly on the shell554    ok, det = True, []555    for H in [1, 2, 3]:556        pts = [q for q in itertools.product(range(-H, H + 1), repeat=4) if any(q)]557        seen, W = set(), []558        for q in pts:                                                   # one representative per rotation (q and −q)559            n = math.sqrt(sum(c * c for c in q)); u = tuple(c / n for c in q)560            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)))561            if key not in seen: seen.add(key); W.append(u)562        bound = 2 * math.asin(1 / H)563        samples = []564        for _ in range(3000):565            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))566        ds = [dist(q, W) for q in samples]567        ok &= max(ds) <= bound + 1e-12                                   # the elementary bound, on every sample568        # the rounding witness itself, on every sample: q = round(H s) is within arcsin(1/H)569        for s_ in samples:570            q = tuple(round(H * c) for c in s_)571            if any(q):572                n = math.sqrt(sum(c * c for c in q)); cosang = abs(sum(a * b for a, b in zip(s_, q))) / n573                ok &= 2 * math.acos(min(1.0, cosang)) <= bound + 1e-12574        measured = covering_radius(W, samples=4000, refine=40)575        ok &= measured <= bound + 1e-9576        det.append(f"H={H}: {len(W)} rotations, measured {math.degrees(measured):.1f}° ≤ bound {math.degrees(bound):.1f}°")577    ok &= 2 * math.asin(1 / 3) < math.acos((3 * 5 ** 0.5 - 1) / 8)     # H = 3 beats every finite subgroup578    # exact composition on the shell p = 73 > 8·3² = 72: the product of two window quaternions reads back from F_73579    p, H = 73, 3580    def qmul_int(a, b):581        w1, x1, y1, z1 = a; w2, x2, y2, z2 = b582        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)583    read = lambda r: r - p if r > (p - 1) // 2 else r                   # the window reading of a residue584    n_pairs = 0585    for _ in range(20000):586        a = tuple(rnd.randint(-H, H) for _ in range(4)); b = tuple(rnd.randint(-H, H) for _ in range(4))587        exact = qmul_int(a, b); modp = qmul_int(tuple(c % p for c in a), tuple(c % p for c in b))588        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 += 1589    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")590    # 1:E6 (p01030)591    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",592          ok, "; ".join(det), kind="CHART")593594    # ---------------------------------------------------------------------------------------------------- the Lie-algebra layer595    # 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.596    def qm(x, y):597        a1, b1, c1, d1 = x; a2, b2, c2, d2 = y598        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)599    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))600    def minv(A, p):601        (a, b, c), (d, e, f), (g, h, i) = A602        det = (a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g)) % p603        if det == 0: return None604        di = pow(det, -1, p)605        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))606        return tuple(tuple((x * di) % p for x in r) for r in adj)607    I3 = ((1, 0, 0), (0, 1, 0), (0, 0, 1))608    def rho(q, p):609        n = sum(c * c for c in q) % p; ni = pow(n, -1, p); cols = []610        for e in ((0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)):611            r = qm(qm(q, e), (q[0], -q[1], -q[2], -q[3])); cols.append([(c * ni) % p for c in r[1:]])612        return tuple(tuple(cols[j][i] for j in range(3)) for i in range(3))613    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)))614    def cayley(a, v, p):615        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)616        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)617    def halfturn(v, p):618        S = Smat(v, p); S2 = mm(S, S, p); ni = pow(sum(c * c for c in v) % p, -1, p)619        return tuple(tuple((I3[i][j] + 2 * ni * S2[i][j]) % p for j in range(3)) for i in range(3))620    def legendre(n, p):621        n %= p622        return 0 if n == 0 else (1 if pow(n, (p - 1) // 2, p) == 1 else -1)623    def nonsquare(p): return next(x for x in range(2, p) if legendre(x, p) == -1)624    def eta_cayley_ok(a, v, p):625        # 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/a626        nu = nonsquare(p); inv = lambda x: pow(x % p, -1, p)627        alpha = (0, inv(nu * a))                                  # 0 + (1/(νa)) η628        S = Smat(v, p)629        for i in range(3):630            for j in range(3):631                h = (0, S[i][j])                                  # η S_ij632                prod = ((alpha[0] * h[0] + nu * alpha[1] * h[1]) % p, (alpha[0] * h[1] + alpha[1] * h[0]) % p)633                if prod != ((S[i][j] * inv(a)) % p, 0): return False634        return True635    def primitive(v):636        g = 0637        for c in v: g = math.gcd(g, abs(c))638        v = tuple(c // g for c in v)639        for c in v:                                               # sign: first nonzero coordinate positive640            if c: return v if c > 0 else tuple(-x for x in v)641    def group_order(M, p):642        k, X = 1, M643        while X != I3: X = mm(X, M, p); k += 1644        return k645646    ok, det = True, []647    n_rot = n_axes = 0648    for p, H in ((13, 1), (17, 1), (37, 2), (73, 3), (101, 3)):649        assert p > 8 * H * H650        window = [q for q in itertools.product(range(-H, H + 1), repeat=4) if any(q[1:])]651        axes = {}652        for q in window:653            a, v = q[0], q[1:]654            R = rho(q, p); n_rot += 1655            ok &= R == (cayley(a % p, v, p) if a % p else halfturn(v, p))656            if a % p: ok &= eta_cayley_ok(a % p, v, p)657            axes.setdefault(primitive(v), set()).add(R)658        for v0, rots in axes.items():659            n_axes += 1660            N = sum(c * c for c in v0) % p661            G = {I3, halfturn(v0, p)} | {C for C in (cayley(a, v0, p) for a in range(1, p)) if C is not None}662            order = p - legendre(N, p)663            ok &= len(G) == order                                                         # the order p − (N(v)|p)664            gen = next((M for M in G if group_order(M, p) == order), None)                # cyclic: an element of full order665            powers, X = {I3}, gen666            while gen is not None and X != I3: powers.add(X); X = mm(X, gen, p)667            ok &= gen is not None and powers == G                                         # whose powers are exactly the set668            for R in rots:                                                                # every iterate of a window rotation stays in it669                X = R670                for _ in range(order): ok &= X in G; X = mm(X, R, p)671            ok &= len(rots) <= 2 * H * H + 2 * H                                         # at most 2H² + 2H window rotations672    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")673    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")674    # 1:E7 (p01032)675    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",676          ok, "; ".join(det))677678    ok, det = True, []679    rnd = random.Random(20261004)680    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])681    def dot(u, w): return sum(x * y for x, y in zip(u, w))682    n_comp = n_half = n_br = 0683    for p in (13, 29, 73, 233):684        inv = lambda x: pow(x % p, -1, p)685        tries = 0686        while n_comp < 4000 * (1 + (13, 29, 73, 233).index(p)) and tries < 100000:687            tries += 1688            u = tuple(rnd.randrange(p) for _ in range(3)); w = tuple(rnd.randrange(p) for _ in range(3))689            uw = dot(u, w) % p690            if (1 + dot(u, u)) % p == 0 or (1 + dot(w, w)) % p == 0: continue691            s = (1 - uw) % p692            m = tuple((u[i] + w[i] + cross(u, w)[i]) % p for i in range(3))693            lhs = mm(rho((1,) + u, p), rho((1,) + w, p), p)694            if s:695                c = tuple((x * inv(s)) % p for x in m)696                if (1 + dot(c, c)) % p == 0: continue697                ok &= lhs == rho((1,) + c, p); n_comp += 1698            elif any(m) and dot(m, m) % p:699                ok &= lhs == halfturn(m, p); n_half += 1700            br = tuple((x - y) % p for x, y in zip(qm((0,) + u, (0,) + w), qm((0,) + w, (0,) + u)))701            ok &= br == (0,) + tuple((2 * x) % p for x in cross(u, w)); n_br += 1702        # the half-turn case on purpose: w chosen with u·w = 1703        for _ in range(300):704            u = tuple(rnd.randrange(p) for _ in range(3))705            if not u[0] % p: continue706            w1, w2 = rnd.randrange(p), rnd.randrange(p)707            w = (((1 - u[1] * w1 - u[2] * w2) * inv(u[0])) % p, w1, w2)708            m = tuple((u[i] + w[i] + cross(u, w)[i]) % p for i in range(3))709            if (1 + dot(u, u)) % p == 0 or (1 + dot(w, w)) % p == 0 or not any(m) or dot(m, m) % p == 0: continue710            ok &= mm(rho((1,) + u, p), rho((1,) + w, p), p) == halfturn(m, p); n_half += 1711    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)")712    n_pow = n_rb = 0713    for H in (1, 2, 3):714        for q in itertools.product(range(-H, H + 1), repeat=4):715            if not any(q): continue716            r = (1, 0, 0, 0)717            for m_ in range(1, 7):718                r = qm(r, q); n_pow += 1719                ok &= max(map(abs, r)) <= (2 * H) ** m_720                for p in (13, 101, 1009):721                    reads = all(((c % p) - p if (c % p) > (p - 1) // 2 else (c % p)) == c for c in r)722                    ok &= reads == (2 * max(map(abs, r)) < p); n_rb += 1723    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")724    # 1:E8 (p01033)725    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",726          ok, "; ".join(det))727728if __name__ == "__main__":729    import time730    want = [a.upper() for a in sys.argv[1:]] or sorted(BLOCK)731    bad = [b for b in want if b not in BLOCK]732    if bad: sys.exit(f"no block {', '.join(bad)}: the blocks are {', '.join(sorted(BLOCK))}")733    t0 = time.time()734    for b in want:735        t = time.time(); _run_block(b); print(f"    [block {b}: {time.time() - t:.1f} s]")736    ok = summary(write=(want == sorted(BLOCK)))737    kinds = {}738    for r in RESULTS: kinds[r["kind"]] = kinds.get(r["kind"], 0) + 1739    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 ""))740    sys.exit(0 if ok else 1)741