From 44c094454b174bba943a0dfbb324cfec8bd7c976 Mon Sep 17 00:00:00 2001 From: Matarener Date: Wed, 21 Oct 2026 21:34:38 +0300 Subject: [PATCH] Audit #4100 BIT.md local-ring algebra, row stock, depth, recurrence and wrapper (no new kappa) Co-Authored-By: Claude Opus 5.5 (1M context) --- research/bit-stock-recurrence-audit/README.md | 106 +++++ .../bit-stock-recurrence-audit/SOURCES.json | 30 ++ research/bit-stock-recurrence-audit/TASK.json | 74 ++++ .../bit-stock-recurrence-audit/VERIFY.json | 30 ++ .../bit-stock-recurrence-audit/test_verify.py | 66 ++++ research/bit-stock-recurrence-audit/verify.py | 365 ++++++++++++++++++ 6 files changed, 671 insertions(+) create mode 100644 research/bit-stock-recurrence-audit/README.md create mode 100644 research/bit-stock-recurrence-audit/SOURCES.json create mode 100644 research/bit-stock-recurrence-audit/TASK.json create mode 100644 research/bit-stock-recurrence-audit/VERIFY.json create mode 100644 research/bit-stock-recurrence-audit/test_verify.py create mode 100644 research/bit-stock-recurrence-audit/verify.py diff --git a/research/bit-stock-recurrence-audit/README.md b/research/bit-stock-recurrence-audit/README.md new file mode 100644 index 000000000..2b24598d6 --- /dev/null +++ b/research/bit-stock-recurrence-audit/README.md @@ -0,0 +1,106 @@ +# Audit of #4100 BIT.md: local-ring algebra, row stock, depth, recurrence and wrapper + +**Audit with a positive scoped result; no new κ.** The selected κ is unchanged. + +Pinned #4100 head: `337d04939571a637aeab54a0f58f54fde811663a` +(`research/unconditional-batched-crt/BIT.md`). + +The #4311 coverage map lists BIT.md's "row stock, whole high packets and fixed +tapes" and "weighted leaves and fast header routes" with no #4100-specific +audit. This package checks every finite or algebraic statement in BIT.md +§1–§4. The moments `mu_1`, `mu_(1-c)` and the body word are #4330's and are +used, not re-certified. + +## Results + +| Statement in BIT.md | Outcome | +|---|---| +| §1: `q = 2^89-1` prime; `2·64^3/q < 10^-16` | Holds. Primality is confirmed by Lucas–Lehmer, independently of #4100's Pratt certificate. The bound is 8.5·10^-22. | +| §1: `N(P) = [[I-P,PC],[CP,I-CPC]]` gives `N(P)^2 = I`; disjoint commuting projectors add; nested differences multiply; conjugation by `diag(A,I)` keeps these; split rank over `Z/q^w` | Holds when C is an involution, as NORMAL-FORM.md's coordinate reversal is (proof below). Checked exactly on 150 instances over `Z/3`, `Z/27`, `Z/25`, `Z/49` and `Z/961`, with C the reversal or a random involution. The rank of P mod q equals its split rank. | +| §2: `K(w) = 269241600 N q^(4096(w-1))`, and splitting family, low matrix and high matrix leaves `U = 9475 t K^(D-1)` | Holds exactly, with `269241600 = W·9475` and `N q^(4096(w-1)) = |GL_64(Z/q^w)|`. The order formula is proved below and brute-forced on small rings. Exact big-integer division (N has 364,544 bits) gives U, and every child's stock divides it, for w = 1, 2 and D = 1..3. | +| §2: D monotone; every child `r floor(n/64)`, r ≤ 42, drops D; two coarse levels halve; `h ≤ 2 log2 e + O(1)` | Holds (proof below). Scanned for n ≤ 3·10^5 (3·10^6 in the development run); `max(D − 2 log2 n) = −11`. | +| §3: header routes come from the low-digit carry plus an ordered high-field addition; doubling and halving are a two-piece split/interleave | Holds: the identity `xA = (x_lo A_lo mod q) + q(carry + x_lo A_hi + x_hi A)` is checked over `Z/q^w` on 800 cases, and doubling maps `[0,(M+1)/2)` onto the evens and the rest onto the odds. | +| §3: `F = 32·64^2·ΣH = 118606528512`; the fallback's at most 49,911 children fit `32·64^2` | Holds. | +| §3: the recurrence gives `T ≤ C/(1-mu_1) n + C/(1-mu_t) n^t w^(1-a)` | The induction is correct (below). The equality recursion stays below the envelope for a ∈ {0, a1}, w ∈ {1, 10, 10^6} and n up to 10^9, and it is monotone, which the `max_(j≤f)` step needs. | +| §4: `k = O_q(w log(e+1)) = o(e)`; elementary fallback only at finitely many widths | Holds. With the literal constants, `2k < e` and `n ≥ 64` first hold at `e = 2^19` (k ≈ 1.6·10^5) and keep holding up to `2^400`. | +| §4: stage cost `O(V e^(1-c+c^2-ca))`, `a_next = c(1-c)+ca`, reaching `a2` | Holds exactly. Every cost term (`n`, `n^(1-c) w^(1-a)`, `k`, `w`) is dominated because `a ≤ c`. | + +**Minor wording.** §2 says each child's payload has volume `V/(WN)`. The split +also divides by the high-matrix fibre `q^(4096(w-1))`, so the payload is +`V/(W N q^(4096(w-1)))`. It is smaller, so the geometric-series step holds a +fortiori. + +## Proofs + +**`N(P)` identities.** Let `C^2 = I`. +- For an idempotent P, the top-left block of `N(P)^2` is + `(I-P)^2 + PC·CP = I - P + P = I`. The top-right block is + `(I-P)PC + PC(I-CPC) = PC - P^2C = 0`, and the other two blocks are + symmetric. +- For idempotents with `P1 P2 = P2 P1 = 0`, multiplying `N(P1) N(P2)` + block by block gives `N(P1+P2)`, because every cross term carries + `P1 P2` or `P2 P1`. +- For `Q = PQ = QP`, the same expansion gives `N(P) N(Q) = N(P-Q)`. +- Conjugation by `diag(A,I)` is an algebra automorphism, so it keeps all + three identities. +- Over the local ring `Z/q^w`, an idempotent's image is free with rank equal + to its rank mod q (idempotents lift uniquely). The characteristic + polynomial is `x^(n-r) (x-1)^r`. + +The controls in `test_verify.py` show that `N(P)^2 = I` fails for a +non-involution C, and that additivity fails for overlapping projectors. + +**`|GL_n(Z/q^w)|`.** Reduction mod q is surjective onto `GL_n(F_q)`. Its +kernel is `{I + qM}`, which has `q^(n^2 (w-1))` elements, because every such +matrix is a unit. So `K(w)/269241600 = |GL_64(Z/q^w)|`, and dividing +`R = t K^D` by `W·|GL_64(Z/q^w)|` leaves `9475 t K^(D-1)`. + +**Depth.** `floor(n/64)` is nondecreasing, so D is nondecreasing. A child +`r floor(n/64)` with `r ≤ 42` has `D ≤ D(42 floor(n/64)) = D(n)-1`. Two levels +multiply the width by at most `(42/64)^2 < 1/2`, so +`D(n) ≤ 2 log2 n + O(1)`. + +**Recurrence.** Put `A = C/(1-mu_1)` and `B = C/(1-mu_t)`, where +`mu_x = (Σ_r H_r (r/64)^x + η F (1/64)^x)/W`. Since `rf ≤ rn/64`, the +children contribute at most `mu_1 A n + mu_t B n^t w^(1-a)`. The `max_(j≤f)` +term is bounded by the envelope at f, which is monotone. The local toll +`C(n + w^(1-a))` is absorbed because `A(1-mu_1) = C` and +`B(1-mu_t) n^t ≥ C`. Leaves need `B ≥ C`. + +**Wrapper exponents.** With `w = e^c` and `n ≈ e^(1-c)`, the term +`n^(1-c) w^(1-a)` gives `e^(1-c+c^2-ca)`. The other terms have exponents +`1-c` and `c`, both at most that, because `a ≤ c`. The identity +`(1-c)^2 + c(1-a) = 1 - a_next` gives `a1 = c(1-c)` and +`a2 = c(1-c^2)`. + +## Not covered + +- The PR2809 block construction itself: the encoded continuation stack, + stream serialization, capsule copying and restoration. +- The finite word certificate and owner data, which `evidence/bit-local-ring-current.log` replays. +- The Schwartz–Zippel density step, which #4227/#4228 audit. + +These are prose constructions or are audited elsewhere. + +## Validation + +```sh +python3.12 -B research/bit-stock-recurrence-audit/verify.py +python3.12 -B research/bit-stock-recurrence-audit/test_verify.py +python3.12 scripts/research_package.py --package research/bit-stock-recurrence-audit +``` + +The finite checks support the written proofs; they do not replace them. + +## Credit + +The construction is by Cinna Davis (#4100). The moment certificate is #4330 +(pinterbanget), and the coverage map is #4311 (Ryanson Jonathan). Prepared for +Mark (Matarener) with substantial Claude Code assistance. Apache-2.0. + +## Open next questions + +Does the PR2809 block construction's encoded stack keep the claimed `O_q(V)` +visited space when the continuation is serialized at every recursive cut, at +the new c? First step: bound the parked bytes per ancestor level against +`V/(W N q^(4096(w-1)))^level`. See `TASK.json`. diff --git a/research/bit-stock-recurrence-audit/SOURCES.json b/research/bit-stock-recurrence-audit/SOURCES.json new file mode 100644 index 000000000..5a08d9f46 --- /dev/null +++ b/research/bit-stock-recurrence-audit/SOURCES.json @@ -0,0 +1,30 @@ +{ + "schema_version": 1, + "base_commit": "12fa8c5c844be21bfefb38a19503fbd867a11244", + "files": [ + { + "path": "certificates/selected-result.json", + "usage": "provenance", + "git_blob": "b1fea0c94436b236cc2b7d41a7be032811b94c81", + "sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "path": "docs/research/agent-task-catalog.json", + "usage": "provenance", + "git_blob": "64bfc5af6709ed9811171e0a2108f05e1219cb65", + "sha256": "352a6184713b8843d96fd2c2071eb7aeaa655906b7c51d45ad444912666c8a96" + }, + { + "path": "research/unconditional-batched-crt/BIT.md", + "usage": "input", + "git_blob": "2aa2544fa61259ac48b936b8c906ec9d31d813fa", + "sha256": "83a0cc082e52219a5f142becc8b5b92019e091cb3ede6c5bccbd8508993fc4da" + }, + { + "path": "research/unconditional-batched-crt/NORMAL-FORM.md", + "usage": "input", + "git_blob": "a0e1d5b370fe34174575f9773bf9981188643fc5", + "sha256": "8649127267c0658b699d5bb16a0a5ac3f1bbd7395d2dbd738c1145a0a65ae09c" + } + ] +} diff --git a/research/bit-stock-recurrence-audit/TASK.json b/research/bit-stock-recurrence-audit/TASK.json new file mode 100644 index 000000000..b740c1b44 --- /dev/null +++ b/research/bit-stock-recurrence-audit/TASK.json @@ -0,0 +1,74 @@ +{ + "schema_version": 1, + "task_id": "uniform-transfer-audit", + "priority_revision": 7, + "base_commit": "12fa8c5c844be21bfefb38a19503fbd867a11244", + "original_card": { + "id": "uniform-transfer-audit", + "question": "Does one quantified step of the proposed #4100 all-size construction hold at its actual caller?", + "starting_points": [ + "docs/research/research-family-status.md", + "certificates/selected-result.json", + "research/unconditional-batched-crt/README.md" + ], + "methods": [ + "Prove one missing implication with explicit quantified hypotheses.", + "Find an exact counterexample at the claimed interface boundary." + ], + "acceptance_gate": "Identify the proposition and what would change downstream; do not equate a finite replay with a universal proof.", + "first_step": { + "command": [ + "python3", + "-B", + "scripts/research_baseline.py" + ], + "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY", + "inspect": "Read PROOF.md and AUDIT.md in #4100. Select one distinct obligation: CRT total repair, native/free-algebra lifting, precision/conditioning, or complete fixed-tape cost. This is an author-submitted proof, not an accepted unconditional theorem.", + "scope": "Source-binding orientation only; no scientific replay, theorem acceptance or supplier admission. Budget any heavy replay separately." + }, + "first_experiment": "Check live handoffs; name one proposition, its caller, quantified domain and dependency closure. Independently prove the implication or exhibit a precise gap/counterexample. Prioritize CRT.md total repair or COMPLEX.md free-algebra lift only if unclaimed. Stop at the chosen obligation; do not substitute finite samples or root arithmetic for an all-size proof, or launch the complete package by default.", + "bottleneck_link": "Finite supplier improvements are useful only under the retained transfer obligations.", + "allowed_variations": [ + "One CRT repair, native lift, precision/conditioning or tape-cost proposition, with its exact caller." + ], + "scope_change_triggers": [ + "Replacing a quantified obligation with another finite smoke test.", + "Treating archive intake, an author PASS receipt or a finite composed fixture as discharge of the universal hypotheses." + ], + "overlap_terms": [ + "uniform transfer", + "all-size", + "precision", + "4100", + "unconditional-batched-crt" + ], + "existing_checks": [ + { + "scope": "PR #4100: proposed unconditional construction and finite checks, without independent acceptance.", + "evidence": "research/unconditional-batched-crt/README.md", + "remaining": "Independent proof audit of one caller-bound all-size implication; all other premises remain explicit.", + "status": "author-checked", + "source_commit": "337d04939571a637aeab54a0f58f54fde811663a", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + } + ] + }, + "scope_status": "refined", + "actual_scope": "One #4100 obligation group: every finite or algebraic statement of BIT.md sections 1-4 (local-ring N(P) algebra, row-stock split, depth recursion, header-route arithmetic, recurrence induction, wrapper threshold and stage exponents).", + "scope_change_reason": "The #4311 coverage map lists BIT.md's row stock and header routes with no #4100-specific audit; #4330 covers only the moments, so this package checks the remaining finite and algebraic claims.", + "bottleneck_link": "Side result for the #4100 admission audit: confirms the finite and algebraic layer of the BIT auxiliary theorem; it does not move the selected frontier.", + "budget_minutes": 90, + "stop_rule": "Stop after the finite and algebraic statements of BIT.md sections 1-4 are checked; do not audit the prose PR2809 block construction or re-certify #4330's moments.", + "outcome": "completed", + "evidence": "Written proofs in README.md; verify.py binds sources and checks every finite statement; test_verify.py has rejecting controls.", + "author": "Mark (Matarener)", + "ai_assistance": "Prepared with substantial Claude Code (Claude) assistance.", + "open_questions": [ + { + "question": "Does the PR2809 block construction's serialized continuation stack keep O_q(V) visited space at the new c when it parks at every recursive cut?", + "why_it_matters": "BIT.md's space bound relies on the parked ancestor volumes forming a geometric series including stack encodings.", + "first_step": "Bound the parked bytes per ancestor level against V/(W N q^(4096(w-1)))^level and add the encoding overhead per level.", + "blocked_by": "none" + } + ] +} diff --git a/research/bit-stock-recurrence-audit/VERIFY.json b/research/bit-stock-recurrence-audit/VERIFY.json new file mode 100644 index 000000000..ba4dbe6b9 --- /dev/null +++ b/research/bit-stock-recurrence-audit/VERIFY.json @@ -0,0 +1,30 @@ +{ + "schema_version": 1, + "status": "ready", + "python": "3.12", + "requirements": [], + "inputs": [ + "research/unconditional-batched-crt/BIT.md", + "research/unconditional-batched-crt/NORMAL-FORM.md" + ], + "commands": [ + { + "argv": [ + "{python}", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 120, + "checks": "Bind BIT.md statements; Lucas-Lehmer for 2^89-1 and the constants; N(P) identities over Z/q^w; exact stock split and child divisibility; depth, halving; header-route carry decomposition and doubling; recurrence envelope; wrapper threshold; stage exponents." + }, + { + "argv": [ + "{python}", + "-B", + "{package}/test_verify.py" + ], + "timeout_seconds": 60, + "checks": "Rejecting controls: non-involution C, overlapping projectors, wrong stock factor, one extra rank-1 child, small widths in the wrapper, and an omitted carry all fail." + } + ] +} diff --git a/research/bit-stock-recurrence-audit/test_verify.py b/research/bit-stock-recurrence-audit/test_verify.py new file mode 100644 index 000000000..bf742297e --- /dev/null +++ b/research/bit-stock-recurrence-audit/test_verify.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +"""Rejecting controls for verify.py.""" +import random +import unittest + +import verify as v + + +class Controls(unittest.TestCase): + def test_non_involution_breaks_square(self): + rng = random.Random(5) + p, m, n = 5, 25, 5 + broke = 0 + for _ in range(20): + S, Si = v.unit(n, p, m, rng) + P = v.mm(v.mm(S, [[int(i == j and i < 2) for j in range(n)] for i in range(n)], m), Si, m) + C, _ = v.unit(n, p, m, rng) + if v.mm(C, C, m) == v.eye(n): + continue + NP = v.N(P, C, m) + broke += v.mm(NP, NP, m) != v.eye(2 * n) + self.assertGreater(broke, 10) + + def test_overlapping_projectors_do_not_add(self): + rng = random.Random(6) + p, m, n = 7, 49, 5 + S, Si = v.unit(n, p, m, rng) + proj = lambda idx: v.mm(v.mm(S, [[int(i == j and i in idx) for j in range(n)] + for i in range(n)], m), Si, m) + C = v.reversal(n) + P, Q = proj([0, 1]), proj([1, 2]) + self.assertNotEqual(v.N(v.add(P, Q, m), C, m), v.mm(v.N(P, C, m), v.N(Q, C, m), m)) + + def test_wrong_stock_factor_fails(self): + Nq = v.gl_order(64, v.q) + K = 269241600 * Nq + U, r = divmod(K ** 2, v.W * Nq) + self.assertEqual(r, 0) + self.assertNotEqual(U, 9474 * K) + + def test_moment_with_extra_rank_one_child_fails(self): + saved = dict(v.H) + try: + v.H[1] += 1 + self.assertGreater(v.mu(1 - float(v.c)), 1.0) # #4330's sensitivity remark + finally: + v.H.clear() + v.H.update(saved) + + def test_wrapper_fails_for_small_widths(self): + e, w, k = v.wrapper_k(10) + self.assertFalse(2 * k < e) + + def test_carry_omitted_breaks_route(self): + p, m = 5, 125 + x, A = [7, 99], [[13, 88], [41, 3]] + y = [sum(x[i] * A[i][j] for i in range(2)) % m for j in range(2)] + xl, xh = [t % p for t in x], [t // p for t in x] + Al, Ah = [[t % p for t in r] for r in A], [[t // p for t in r] for r in A] + lo = [sum(xl[i] * Al[i][j] for i in range(2)) % p for j in range(2)] + hi = [sum(xl[i] * Ah[i][j] + xh[i] * A[i][j] for i in range(2)) % 25 for j in range(2)] + self.assertNotEqual(y, [a + p * b for a, b in zip(lo, hi)]) + + +if __name__ == "__main__": + unittest.main() diff --git a/research/bit-stock-recurrence-audit/verify.py b/research/bit-stock-recurrence-audit/verify.py new file mode 100644 index 000000000..d6253d78d --- /dev/null +++ b/research/bit-stock-recurrence-audit/verify.py @@ -0,0 +1,365 @@ +#!/usr/bin/env python3 +"""Audit #4100 BIT.md: local-ring algebra, stock, depth, recurrence and wrapper. + +Propositions checked (BIT.md at #4100 head 337d049): +1. constants: q = 2^89-1 prime (Lucas-Lehmer), 2*64^3/q < 10^-16, + 269241600 = W*9475, F = 32*64^2*sum(H), a1 and a2; +2. section 1: N(P) = [[I-P,PC],[CP,I-CPC]] squares to I, disjoint commuting + projectors add, nested differences multiply, conjugation by diag(A,I) + keeps them, and split rank equals rank mod q, over Z/q^w; +3. section 2: K(w) = 269241600 |GL_64(Z/q^w)|, the split leaves + U = 9475 t K^(D-1) and every child's stock divides U; the depth D is + monotone, every child drops it and two levels halve the width; +4. section 3: the header-route carry decomposition and the two-piece + doubling map; the recurrence envelope; +5. section 4: k = O(w log e) with the literal constants, and the stage + exponents a_next = c(1-c)+ca reaching a2. +Standard library only. +""" +import itertools +import math +import random +from fractions import Fraction as Q +from pathlib import Path + +PKG = Path(__file__).resolve().parent +BIT = PKG.parents[1] / "research" / "unconditional-batched-crt" / "BIT.md" +NF = PKG.parents[1] / "research" / "unconditional-batched-crt" / "NORMAL-FORM.md" + +q = 2 ** 89 - 1 +c = Q(8667, 5_000_000) +W = 28416 +H = {1: 381024, 2: 400416, 3: 22752, 4: 26880, 5: 29376, 6: 17376, 7: 26784, 42: 288} +ETA = Q(1, 10 ** 16) + + +def check_sources(): + bit = BIT.read_text(encoding="utf-8") + for needle in ["q=2^89-1,", "c=8667/5000000,", "2*64^3/q < 10^-16.", + "N(P)=[[I-P,PC],[CP,I-CPC]].", + "Then N(P)^2=I, nested differences multiply as their corresponding residuals,", + "and disjoint commuting projectors add. These identities survive conjugation", + "Let W=28416, N=|GL_64(F_q)| and", + "K(w)=269241600 N q^(4096(w-1)),", + "D(n)=0 (n<64), D(n)=1+D(42 floor(n/64)) (n>=64).", + "low matrix and high matrix fields leaves U=t*9475*K(w)^(D(n)-1) fresh rows.", + "Two coarse levels reduce width by at least a factor two, so h<=2log2(e)+O(1).", + "Each child's payload has volume V/(WN).", + "is obtained from the low-digit carry, followed by an ordered high-field", + "Doubling/halving on an odd modulus is a two-piece split/interleave; it does", + "H={1:381024,2:400416,3:22752,4:26880,5:29376,6:17376,7:26784,42:288}.", + "With F=32*64^2*sum_r(H_r)=118606528512, t=1-c and eta=10^-16, normalized runtime", + "T(n,w) <= C/(1-mu_1) n + C/(1-mu_t) n^t w^(1-a)", + "a_next=c(1-c)+ca. Apply this first to a=0 and then to the completed a1,"]: + assert needle in bit, needle + assert "C the m-coordinate reversal" in NF.read_text(encoding="utf-8") + + +# ------------------------------------------------------------- constants +def gl_order(n, p): + out = 1 + for i in range(n): + out *= p ** n - p ** i + return out + + +def check_constants(): + s = 4 + for _ in range(87): + s = (s * s - 2) % q + assert s == 0, "Lucas-Lehmer" + assert Q(2 * 64 ** 3, q) < ETA + assert 269241600 == W * 9475 + assert 32 * 64 ** 2 * sum(H.values()) == 118606528512 + assert c * (1 - c) == Q(173039532444, 10 ** 14) + assert c * (1 - c * c) == Q(1733394791695384296, 10 ** 21) + assert max(H) == 42 and 49911 <= 32 * 64 ** 2 + for n, p, w in ((1, 3, 3), (2, 3, 2), (2, 2, 3), (2, 5, 2)): + m = p ** w + cnt = 0 + for ent in itertools.product(range(m), repeat=n * n): + det = ent[0] if n == 1 else ent[0] * ent[3] - ent[1] * ent[2] + cnt += det % p != 0 + assert cnt == gl_order(n, p) * p ** (n * n * (w - 1)), (n, p, w) + + +# --------------------------------------------------- section 1 algebra +def mm(A, B, m): + return [[sum(a * b for a, b in zip(r, col)) % m for col in zip(*B)] for r in A] + + +def eye(n): + return [[int(i == j) for j in range(n)] for i in range(n)] + + +def sub(A, B, m): + return [[(a - b) % m for a, b in zip(r, s)] for r, s in zip(A, B)] + + +def add(A, B, m): + return [[(a + b) % m for a, b in zip(r, s)] for r, s in zip(A, B)] + + +def inverse(A, p, m): + k = len(A) + A = [r[:] for r in A] + I = eye(k) + for col in range(k): + piv = next((r for r in range(col, k) if A[r][col] % p), None) + if piv is None: + return None + A[col], A[piv] = A[piv], A[col] + I[col], I[piv] = I[piv], I[col] + iv = pow(A[col][col], -1, m) + A[col] = [x * iv % m for x in A[col]] + I[col] = [x * iv % m for x in I[col]] + for r in range(k): + if r != col and A[r][col]: + f = A[r][col] + A[r] = [(a - f * b) % m for a, b in zip(A[r], A[col])] + I[r] = [(a - f * b) % m for a, b in zip(I[r], I[col])] + return I + + +def rank_mod(A, p): + rows = [[v % p for v in r] for r in A] + n, rk = len(rows[0]), 0 + for col in range(n): + piv = next((r for r in range(rk, len(rows)) if rows[r][col]), None) + if piv is None: + continue + rows[rk], rows[piv] = rows[piv], rows[rk] + iv = pow(rows[rk][col], -1, p) + rows[rk] = [v * iv % p for v in rows[rk]] + for r in range(len(rows)): + if r != rk and rows[r][col]: + f = rows[r][col] + rows[r] = [(a - f * b) % p for a, b in zip(rows[r], rows[rk])] + rk += 1 + return rk + + +def unit(n, p, m, rng): + while True: + S = [[rng.randrange(m) for _ in range(n)] for _ in range(n)] + Si = inverse(S, p, m) + if Si: + return S, Si + + +def block(A, B, C, D): + return [ra + rb for ra, rb in zip(A, B)] + [rc + rd for rc, rd in zip(C, D)] + + +def N(P, C, m): + I = eye(len(P)) + return block(sub(I, P, m), mm(P, C, m), mm(C, P, m), sub(I, mm(mm(C, P, m), C, m), m)) + + +def reversal(n): + return [[int(i + j == n - 1) for j in range(n)] for i in range(n)] + + +def check_local_ring(trials=30, seed=4100): + rng = random.Random(seed) + count = 0 + for p, w, n in [(3, 1, 4), (3, 3, 5), (5, 2, 6), (7, 2, 6), (31, 2, 5)]: + m = p ** w + for t in range(trials): + S, Si = unit(n, p, m, rng) + + def proj(idx): + D = [[int(i == j and i in idx) for j in range(n)] for i in range(n)] + return mm(mm(S, D, m), Si, m) + + if t % 2: + C = reversal(n) + else: + T, Ti = unit(n, p, m, rng) + sg = [[(rng.choice([1, m - 1]) if i == j else 0) for j in range(n)] for i in range(n)] + C = mm(mm(T, sg, m), Ti, m) + assert mm(C, C, m) == eye(n) + k = rng.randrange(1, n) + idx = rng.sample(range(n), k) + P = proj(idx) + I2 = eye(2 * n) + NP = N(P, C, m) + assert mm(NP, NP, m) == I2 + a, b = idx[:max(1, k // 2)], idx[max(1, k // 2):] + P1, P2 = proj(a), proj(b) + assert N(add(P1, P2, m), C, m) == mm(N(P1, C, m), N(P2, C, m), m) + assert mm(NP, N(P1, C, m), m) == N(sub(P, P1, m), C, m) + Ad = [rng.choice([x for x in range(1, m) if x % p]) for _ in range(n)] + Z = [[0] * n for _ in range(n)] + Dg = block([[Ad[i] if i == j else 0 for j in range(n)] for i in range(n)], Z, Z, eye(n)) + Dgi = block([[pow(Ad[i], -1, m) if i == j else 0 for j in range(n)] for i in range(n)], Z, Z, eye(n)) + Mc = mm(mm(Dg, NP, m), Dgi, m) + assert mm(Mc, Mc, m) == I2 + assert mm(Mc, mm(mm(Dg, N(P1, C, m), m), Dgi, m), m) == \ + mm(mm(Dg, N(sub(P, P1, m), C, m), m), Dgi, m) + assert rank_mod(P, p) == k + count += 1 + return count + + +# ------------------------------------------- section 2 stock and depth +def D(n): + d = 0 + while n >= 64: + n, d = 42 * (n // 64), d + 1 + return d + + +def check_stock(): + Nq = gl_order(64, q) + for w in (1, 2): + K = 269241600 * Nq * q ** (4096 * (w - 1)) + split = W * Nq * q ** (4096 * (w - 1)) # family x |GL_64(Z/q^w)| + for d in (1, 2, 3): + for t in (1, 7): + U, r = divmod(t * K ** d, split) + assert r == 0 and U == t * 9475 * K ** (d - 1) + for dc in range(d): + assert U % (9475 * t * K ** dc) == 0 + return Nq.bit_length() + + +def check_depth(limit=300_000): + prev, worst = 0, -1e9 + for n in range(1, limit + 1): + d = D(n) + assert d >= prev + prev = d + if n >= 64: + f = n // 64 + assert 42 * f < n and all(D(r * f) <= d - 1 for r in H) + assert 2 * 42 * ((42 * f) // 64) <= n + worst = max(worst, d - 2 * math.log2(n)) + return worst + + +# --------------------------------------- section 3 routes and recurrence +def check_routes(seed=89): + rng = random.Random(seed) + cases = 0 + for p, w, n in [(3, 2, 4), (5, 3, 4), (7, 3, 3), (101, 2, 4)]: + m, mh = p ** w, p ** (w - 1) + for _ in range(200): + A = [[rng.randrange(m) for _ in range(n)] for _ in range(n)] + x = [rng.randrange(m) for _ in range(n)] + y = [sum(x[i] * A[i][j] for i in range(n)) % m for j in range(n)] + xl, xh = [v % p for v in x], [v // p for v in x] + Al, Ah = [[v % p for v in r] for r in A], [[v // p for v in r] for r in A] + raw = [sum(xl[i] * Al[i][j] for i in range(n)) for j in range(n)] + ylo, carry = [v % p for v in raw], [v // p for v in raw] + yhi = [(carry[j] + sum(xl[i] * Ah[i][j] + xh[i] * A[i][j] for i in range(n))) % mh + for j in range(n)] + assert y == [lo + p * hi for lo, hi in zip(ylo, yhi)] + cases += 1 + for M in (3, 7, 9, 25, 343, 2 ** 13 - 1, 101 ** 2): + img = [2 * x % M for x in range(M)] + half = (M + 1) // 2 + assert img[:half] == list(range(0, M, 2)) and img[half:] == list(range(1, M, 2)) + return cases + + +def mu(x, eta=float(ETA)): + F = 32 * 64 ** 2 * sum(H.values()) + return (sum(h * (r / 64) ** x for r, h in H.items()) + eta * F * (1 / 64) ** x) / W + + +def recurrence_T(n, ww, memo, eta=float(ETA)): + if n < 64: + return ww + if n in memo: + return memo[n] + f = n // 64 + F = 32 * 64 ** 2 * sum(H.values()) + v = sum(h / W * recurrence_T(r * f, ww, memo, eta) for r, h in H.items()) + v += eta * F / W * recurrence_T(f, ww, memo, eta) + n + ww + memo[n] = v + return v + + +def check_recurrence(): + t = 1 - float(c) + m1, mt = mu(1.0), mu(t) + assert m1 < 1 and mt < 1 + A, B = 1 / (1 - m1), 1 / (1 - mt) + worst = 0.0 + for a in (0.0, float(c * (1 - c))): + for w in (1, 10, 10 ** 6): + ww = w ** (1 - a) + memo = {} + for n in (64, 100, 4096, 10 ** 5, 10 ** 7, 10 ** 9): + worst = max(worst, recurrence_T(n, ww, memo) / (A * n + B * n ** t * ww)) + assert worst <= 1 + memo = {} + vals = [recurrence_T(n, 1.0, memo) for n in range(1, 100_000)] + assert all(vals[i] <= vals[i + 1] for i in range(len(vals) - 1)) + return m1, mt, worst + + +# ------------------------------------------------ section 4 wrapper +def wrapper_k(L): + lq = math.log(q) + lqN = 4096 + sum(math.log1p(-(q ** (i - 64))) for i in range(64)) / lq + e = 2.0 ** L + w = math.ceil(e ** float(c)) + d, n = 0, e + while n >= 64: + n, d = 42 * (n // 64), d + 1 + k = math.ceil(d * (math.log(269241600) / lq + lqN + 4096 * (w - 1))) + return e, w, k + + +def check_wrapper(): + first = None + for L in range(1, 400): + e, w, k = wrapper_k(L) + good = 2 * k < e and (e - k) // w >= 64 + if good and first is None: + first = L + if first is not None: + assert good, L + return first + + +def check_exponents(): + for a in (Q(0), c * (1 - c)): + main = 1 - c + c * c - c * a + for term in (1 - c, (1 - c) ** 2 + c * (1 - a), c): + assert term <= main + assert (1 - c) ** 2 + c * (1 - a) == 1 - (c * (1 - c) + c * a) + assert c * (1 - c) + c * (c * (1 - c)) == c * (1 - c * c) + + +def main(): + check_sources() + print("PASS source binding: BIT.md statements and NORMAL-FORM.md's reversal C") + check_constants() + print("PASS constants: 2^89-1 prime (Lucas-Lehmer); 2*64^3/q < 1e-16; 269241600 = W*9475; " + "F = 118606528512; a1, a2 exact; |GL_n(Z/p^w)| lift formula on small cases") + n = check_local_ring() + print(f"PASS local ring: N(P)^2 = I, disjoint add, nested multiply, diag(A,I) conjugation, " + f"split rank = rank mod q on {n} instances over Z/q^w (C reversal or random involution)") + bits = check_stock() + print(f"PASS stock: R = t K^D splits by W x |GL_64(Z/q^w)| into U = 9475 t K^(D-1); " + f"every child stock divides U (exact; |GL_64(F_q)| has {bits} bits)") + worst = check_depth() + print(f"PASS depth: D monotone, children drop D, two levels halve (n <= 3e5); " + f"max D(n) - 2 log2 n = {worst:.2f}") + cases = check_routes() + print(f"PASS header routes: carry decomposition on {cases} cases; doubling is a two-piece interleave") + m1, mt, ratio = check_recurrence() + print(f"PASS recurrence: mu_1 = {m1:.10f}, mu_t = {mt:.10f}; equality recursion " + f"<= envelope (max ratio {ratio:.2e}) and monotone") + first = check_wrapper() + print(f"PASS wrapper: 2k < e and n >= 64 from log2 e = {first} on (checked to 2^400)") + check_exponents() + print("PASS exponents: every stage cost term <= 1-c+c^2-ca; a_next reaches a2 exactly") + print("PASS_BIT_STOCK_RECURRENCE_AUDIT") + + +if __name__ == "__main__": + main()